20100319-第三講：求導積分與微分方程數(shù)值解（2次課）

內(nèi)容：本講針對高等數(shù)學一元微積分學補充極限、導數(shù)、積分相關(guān)運算；介紹Funtool符號計算器；目的：學習極限/導數(shù)/積分相關(guān)函數(shù)的指令實現(xiàn)，為學習微分方程數(shù)值解作準備；要求：能夠解決高等數(shù)學中的一類極限/導數(shù)/積分求解問題；了解并會使用Funtool符號計算器；掌握極限(左、右極限)函數(shù)limit掌握導數(shù)(1階導、高階導、偏導)函數(shù)diff掌握積分(不定積分、定積分、數(shù)值積分)函數(shù)

int

trapz

quad

quadlquad8第三講極限、導數(shù)、積分(補充)求極限、求導數(shù)與求積分...

極限,導數(shù),積分是我們在高等數(shù)學學習中接觸過的最基本也是最重要的概念.一方面它們是很多數(shù)學工具的基礎(chǔ)(比如微分方程);另一方面它們又是工程計算和科學研究直接面對的問題.微分(導數(shù))運算比較簡單,任何一個由基本初等函數(shù)經(jīng)過四則及復合運算構(gòu)成的函數(shù),都可以用導數(shù)公式和求導法則算出它們的導數(shù).積分運算則相對復雜得多,仍有許多函數(shù)“積不出來”,由于它們的原函數(shù)無法由基本初等函數(shù)經(jīng)過四則及復合運算構(gòu)成,計算這類定積分問題我們也只能采用數(shù)值方法.

借助MATLAB我們得以快速解決這些問題!基本調(diào)用格式：limit(f)功能:計算limit(f,x,a)功能:計算limit(f,x,inf)功能:計算limit(f,x,a,'right')功能:計算limit(f,x,a,'left')功能:計算求極限運算的調(diào)用格式注意：默認x趨于0；在左，右極限不相等，或有一個不存在時，默認為求右極限；求極限運算的應用示例應用示例（熟悉應用類型）：例1求極限syms

x;y=((1+tan(x))/(1+sin(x)))^(1/x^3);limit(y)例2求極限syms

n;y=(1+1/n)^n;limit(y,n,inf)

例3求極限symsx;y=5*x+log(sin(x)+exp(sin(x)));limit(y,x,3,'left')求導數(shù)運算的調(diào)用格式[1]一元函數(shù)求導基本調(diào)用格式：diff(f)功能-求函數(shù)f的一階導數(shù)diff(f,n)功能-求函數(shù)f的n階導數(shù)應用示例：例4求的一階、二階導數(shù)symsabx;y=(a*x+tan(3*x))^(1/2)+sin(x)*cos(b*x);d1y=diff(y),disp('***'),pretty(d1y),disp('***')d2y=diff(y,2),disp('***'),pretty(d2y),求導數(shù)運算的調(diào)用格式[2]多項式擬合求導（表達式未知或不易求導）方法說明：先利用polyfit將函數(shù)擬合成多項式函數(shù)，然后利用多項式函數(shù)求導命令polyder求導或diff求導應用示例：例5用5階多項式擬合函數(shù)并求x=2處的二階導函數(shù)值x=0:.1:8;y=cos(x).*log(3+x.^2+exp(x.^2));p=polyfit(x,y,5),y2=polyval(p,x);plot(x,y,'b',x,y2,'r');legend('y','y2',2);%產(chǎn)生數(shù)據(jù)點，擬合成5階多項式函數(shù)，并作圖比較p1=polyder(p);p2=polyder(p1);ans1=polyval(p2,2),%利用多項式函數(shù)專用求導函數(shù)polyder求導，并代值y2=poly2sym(p,'x'),y2d2=diff(y2,2),ans2=subs(y2d2,2),%利用通用求導函數(shù)diff求導，并代值求導數(shù)運算的調(diào)用格式[3]參數(shù)方程求導方法說明：對參數(shù)方程x=x(t);y=y(t);先求出dy/dt和dx/dt然后代入公式dy/dx=dy/dt/dx/dt

即可應用示例：例6求參數(shù)方程syms

t;x=t*(1-sin(t));y=t*cos(t);ezplot(x,y);gridon;dx=diff(x,t);dy=diff(y,t);dydx=dy/dx;pretty(dydx)%下面在t=4.1處作出參數(shù)方程的切線（導數(shù)）holdon;t=4.1;x=eval(x);y=eval(y);plot(x,y,'ro');k=eval(dydx);line([x,x+1],[y,y+k],'color','r')求導數(shù)運算的調(diào)用格式[4]多元函數(shù)求導方法說明：對指定變量求導，求偏導數(shù)應用示例：例7求對z的偏導數(shù)symsabxyz;u=a*exp(b*x+y+z^2);pretty(diff(u,z))例8對symsxy;z=x^3*y^2+sin(x*y);diff(z,x,3)求導數(shù)運算的應用示例例9以為例驗證羅必塔法則：symsabxf=a^x-b^x;g=x;l1=limit(f/g,x,0)df=diff(f,x);dg=diff(g,x);l2=limit(df/dg,x,0)ifl1==l2disp('羅必塔法則得到驗證！')end求不定積分運算的調(diào)用格式[1]不定積分方法說明：int(f)對默認變量積分；int(f,v)對指定變量積分應用示例：例10計算syms

x;y=1/(sin(x)^2*cos(x)^2);pretty(int(y))例11計算symsax;y=1/(a^2-x^2);pretty(int(y,x))例12計算二重不定積分symsxy;F=int(int(x*exp(-x*y),x),y)求定積分運算的調(diào)用格式[2]定積分-解析解法方法說明：int(f,x,a,b)依據(jù)微積分基本公式計算應用示例：例13計算symsax;f=sqrt(x^2+a);pretty(int(f,x,-2,2))例14對變上限函數(shù)求導symstx;f=sqrt(1-t^2);pretty(diff(int(f,t,0,x^2)))求定積分運算的調(diào)用格式[3]定積分-數(shù)值解法方法說明：當定積分-符號解法失效時，必須用定積分-數(shù)值解法來近似計算定積分的值。矩形公式sum，復合梯形公式trapz，復合辛普森公式quad/quad8的區(qū)別在于替代等距曲邊梯形的方式不同：求定積分運算的應用示例應用示例：sum使用一次用于求向量或矩陣每一列的和，若使用兩次則先按列求和再按行求和(行列總和)例15矩形法計算在x=0與x=10之間所圍面積dx=0.1;x=0:dx:10;y=-x.^2+115;sum(y(1:length(x)-1))*dx(的近似值)求定積分運算的調(diào)用格式trapz(x,y)用復合梯形公式計算定積分，x為積分變量分點向量，y為被積函數(shù)分點函數(shù)值向量quad('fun',a,b,tol,trace)用復合辛普森公式計算定積分，fun為被積函數(shù)表達式字符串或m函數(shù)文件名，a，b是積分下上限，tol表示精度(缺省0.001)，trace=1圖示積分過程(默認=0不顯示)%quadl采用Lobatto算法，精度和速度要優(yōu)于quad

%quad8采用8階NewtonCotes算法，精度優(yōu)于quad求定積分運算的應用示例例16用兩種方法求定積分x=2:.1:5;y=log(x)./(x.^2);tt=trapz(x,y)%復合梯形公式

fun=inline('log(x)./(x.^2)','x');ss=quad(fun,2,5)%復合辛普森公式Funtool符號計算器-界面Funtool符號計算器-功能圖形化符號函數(shù)計算器的使用：f=為圖形窗口1的控制函數(shù)，其缺省值為x；g=為圖形窗口2的控制函數(shù)，其缺省值為1；x=為兩窗口函數(shù)的自變量取值范圍，缺省[-2*pi,2*pi]a=為常數(shù)，缺省值為1/2。df/dx

計算函數(shù)f對x的導法式，并賦給f。intf

計算函數(shù)f的積分函數(shù)，并賦給f。simplef計算函數(shù)f的最簡表達式，并賦給f。(symsx)simplify(cos(x)^2+sin(x)^2);simplify((x^2+5*x+6)/(x+2));

expand(cos(x+y));expand((x-2)*(x-4));

symsxy;factor(x^3-y^3);factor(x^3+3*x^2+3*x+1);numf

取表達式f的分子，并賦給f。denf

取表達式f的分母，并賦給f。1/f

求f的倒數(shù)函數(shù)，并賦給f。finv

求f的反函數(shù)，并賦給f。Funtool符號計算器-功能f±a

計算f(x)±a，并賦給f。f*a

計算f(x)*a，并賦給f。f/a

計算f(x)/a，并賦給f。f^a

計算f^a，并賦給f。f(x+a)

計算f(x+a)，并賦給f。f(x*a)

計算f(ax)，并賦給f。f+a

計算f(x)+a，并賦給f。f±g

計算兩函數(shù)之和/差，并賦給f。f*g

計算兩函數(shù)之積，并賦給f。f/g

計算兩函數(shù)之比，并賦給f。f(g)

計算復合函數(shù)f(g(x))。g=f

將f的函數(shù)值賦給g。swap

交換f與g的函數(shù)表達式。That’sall~3Q!

第三講微分方程數(shù)值解內(nèi)容：本講首先演練單擺微分方程求解的全過程；隨后由例題入手介紹基于MATLAB的微分方程求解函數(shù)，然后重點講解微分方程(組)的圖形圖像解法；最后介紹歐拉方法、改進歐拉方法，R-K方法目的：掌握微分方程數(shù)值解的一般思路和方法。要求：能夠處理應用類型微分方程數(shù)值解問題。掌握單擺微分方程求解的完整過程(課本引例)掌握基于常用微分方程求解函數(shù)dsolveode掌握圖形圖像求解:斜率場/相平面/等值線了解歐拉方法、改進歐拉方法、R-K方法大多數(shù)微分方程無法求解析解?微分方程是研究函數(shù)變化規(guī)律的有力工具，在科技、工程、經(jīng)濟管理、及生態(tài)、環(huán)境、人口、交通等各個領(lǐng)域有著廣泛的應用。

建立微分方程可以依據(jù)物理的、或其他原理和規(guī)律建立的平衡關(guān)系，但是！更重要的問題是如何求解這些微分方程(組)部分微分方程可以求得解析解，但是絕大多數(shù)的非線性、變系數(shù)微分方程或“難以求解”或“求不出解”，所以對于實際問題，研究微分方程的數(shù)值解具有重要意義！我們所熟悉的微分方程？微分方程初值問題的最簡單形式：簡單定義：含有導數(shù)的方程就稱為微分方程一般形式：解析解：求得具體解析式y(tǒng)=f(x)(早先學過的...)數(shù)值解：求得系列散點xi對應的近似值yi(表格法)圖像解：用圖像表示解曲線(有何優(yōu)勢？)表示函數(shù)的三種方法？單擺微分方程求解：建立方程提示：這是一個完整的微分方程建立、求解的過程，通過本例的學習，要求完成p58容器刻度問題全程求解?。ū局v實驗題）

lmg由牛頓第二定律建立方程：我們要找到符合條件的theta與t的函數(shù)關(guān)系，但是此2階非線性微分方程并不容易求得解析解…除非單擺微分方程求解：求近似解除非一簡化方程求其近似解：取x0=10°即x0=

0.1745，在此弧度范圍內(nèi)sin

所以原微分方程可以簡化為：此線性常系數(shù)微分方程可以直接用dsolve函數(shù)求得：dsolve('D2theta+g/l*theta=0','theta(0)=a0','Dtheta(0)=0','t')其解析解為：a0*cos((g/l)^(1/2)*t)這個解可以作為原方程的近似解單擺微分方程求解：求數(shù)值解除非二求其數(shù)值解：也就是部分點xi對應的近似值yi我們采用ode23函數(shù)求解，所以首先改寫方程：由簡化方程求得近似解為：代入g=9.8l=25，得到周期T

10s，下面考察ts=0到tf=10內(nèi)若干點處的近似值i，即所謂的數(shù)值解單擺微分方程求解：求數(shù)值解首先建立被調(diào)函數(shù)danbai.mfunction

xdot=danbai(t,x)g=9.8;l=25;xdot(1)=x(2);xdot(2)=-g/l*sin(x(1));xdot=xdot';然后是主調(diào)指令，也可寫成主調(diào)文件loaddanbai.mwarningoffts=0;tf=10;a0=0.1745;cond0=[a0,0];%初始化變量[t,x]=ode23('danbai',ts,tf,cond0);%調(diào)用ode23函數(shù)求解g=9.8;l=25;w=sqrt(g/l);y=a0*cos(w*t);%近似解[t,x(:,1),y]%輸出t對應的數(shù)值解和近似解subplot(1,2,2);stem(t,x(:,1),'ro');title('數(shù)值解')subplot(1,2,1);holdon;stem(t,y,'bp');plot(t,y,'b-');title('近似解')用dsolve函數(shù)求解微分方程MATLAB求解微分方程解析解的函數(shù)dsolveSymbolicsolutionofordinarydifferentialequations.Syntax~r=dsolve('eq1','eq2',...,'cond1','cond2',...,'v')題例：p49-4.4.1/ex1，ex2dsolve('Dy=1+y^2')dsolve('Dtheta=1+theta^2','theta(0)=1','xi')dsolve('x^2*D2y+x*Dy+(x^2-(1/2)^2)*y=0','y(pi/2)=2','Dy(pi/2)=-2/pi','x')pretty(ans)提示：一些需要注意的細節(jié)…用dsolve函數(shù)求解微分方程組MATLAB求解微分方程組解析解的函數(shù)dsolve題例：p50-4.4.1/ex3，ex4[f,g]=dsolve('Df=3*f+4*g','Dg=-4*f+3*g')[f,g]=dsolve('Df=3*f+4*g','Dg=-4*f+3*g','f(0)=0','g(0)=1','x')下面的指令有否區(qū)別？dsolve('Dy=x*sin(x)/cos(y)')dsolve('Dy=x*sin(x)/cos(y)','x')提示：用dsolve求解存在解析解的微分方程相當方便，在“只要結(jié)果，不求過程”的場合，節(jié)約了大量時間。練習：Malthus人口模型計算Malthus認為單位時間內(nèi)人口凈增長率為常數(shù)：

dsolve('Dx=r*x','x(0)=x0','t')

d=1790:10:1900;x=[3.95.37.29.612.917.123.231.438.650.262.976];t=(d-1790)/10;y=log(x);p=polyfit(t,y,1);y=poly2str(p,'t')r=p(1),a=p(2);x0=exp(a)plot(t,x,'r+');holdon;t=min(t):0.01:max(t);x=x0*exp(r*t);plot(t,x,'b-');[x,y]=ode23('tbp51',[pi/2,pi],[2,-2/pi])plot(x,y(:,1),x,y(:,2));holdon;stem(x,y(:,1))title('x^2*D2y+x*Dy+(x^2-n^2)*y,cond0')loadtbp51.m基于R-K算法的數(shù)值解函數(shù)odeMATLAB求解微分方程數(shù)值解的函數(shù)odeSolveinitialvalueproblemsforordinarydifferentialequations(ODEs).Syntax(ode23):[T,Y]=ode23(odefun,tspan,y0)[T,Y]=ode23(odefun,tspan,y0,options)題例：p51-ex1functiondequ=tbp51(x,y)dequ=[y(2);-y(2)/x+((1/2/x)^2-1)*y(1)];tbp51.m解析法/數(shù)值法的局限...解析法~能找到精確解固然好，但適應面太窄。數(shù)值法~只能得到一些離散點處的近似值，不能較好展示全局和趨勢。圖像法~可能也只有圖像法，才能避免上述缺陷[方法概述]：斜率場法是依據(jù)y'=f(x,y)得到平面上一些點的斜率值，然后過這些點引出自該點出發(fā)的短直線，通過觀察趨勢，了解解曲線的分布和性態(tài)。題例：p45-ex1圖解法-斜率場1用斜率場法求解微分方程：y'=sinx

sinysymsxy;fun=sin(x)*sin(y);hx=16/40;hy=16/40;x0=-8;y0=-8;holdonfori=1:40x=x0+(i-1)*hx;forj=1:40y=y0+(j-1)*hy;k=eval(fun);圖解法-斜率場2ifabs(hx*k)>hy

plot([x,x+hy/k*2/3],[y,y+hy*2/3])elseplot([x,x+hx*2/3],[y,y+hx*k*2/3])endendendtitle('dy/dx=sinx*siny');xlabel('x');ylabel('y');hxhyk圖解法-斜率場3圖解法-相平面軌跡1[方法概述]：相平面軌跡法是依據(jù)不同的初值條件，先用數(shù)值解法求出各自對應的數(shù)值解（x(t),y(t)），最后用plot描點繪圖，通過觀察趨勢，了解解曲線的分布和性態(tài)。題例：p53-ex1先用數(shù)值解法求出若干初值條件下的(x(t),y(t))%tbp53.mfunctiondequ=tbp53(t,x)dequ=[2*x(1)-1*x(1)*x(2);-1*x(2)+1