《解微分方程》PPT课件

1、MATLABMATLABODEODE初值问题的数值解初值问题的数值解PDEPDE问题的数值解问题的数值解问题提出问题提出 倒葫芦外形容器壁上的刻度问题倒葫芦外形容器壁上的刻度问题. .对于如下图圆柱对于如下图圆柱外形容器壁上的容积刻度外形容器壁上的容积刻度, ,可以利用圆柱体体积公式可以利用圆柱体体积公式HDV22其中直径其中直径D D为常数为常数. . 对于几何外形不是规那么的容器对于几何外形不是规那么的容器, ,比如倒葫芦外形比如倒葫芦外形容器壁上如何标出刻度呢容器壁上如何标出刻度呢? ?下表是经过丈量得到部下表是经过丈量得到部分容器高度与直径的关系分容器高度与直径的关系. .x1o根据上

2、表的数据根据上表的数据, ,可以拟合出倒可以拟合出倒葫芦外形容器的图葫芦外形容器的图, ,建立如下图建立如下图的坐标轴的坐标轴, ,问题为如何根据恣意问题为如何根据恣意高度高度x x标出容器体积标出容器体积V V的刻度的刻度, ,由由微元思想分析可知微元思想分析可知H 0 0.2 0.4 0.6 0.8 1.0D 0.04 0.11 0.26 0.56 1.04 1.17dxDdV241其中其中x x表示高度表示高度, ,直径直径D D是高度是高度x x的函数的函数, ,记为记为D(x).D(x).x1o0)0()(412VxDdxdV只需求解上述方程只需求解上述方程, ,就可求出体积就可求出

3、体积V V与高度与高度x x之间的之间的函数关系函数关系, ,从而可标出容器壁上容积的刻度从而可标出容器壁上容积的刻度, ,但问题但问题是函数是函数D(x)D(x)无解析表达式无解析表达式, ,无法求出其解析解无法求出其解析解. . 因此因此, , 得到如下微分方程初值问题得到如下微分方程初值问题包含自变量、未知函数及未知函数的导数或微分的包含自变量、未知函数及未知函数的导数或微分的方程称为微分方程。方程称为微分方程。在微分方程中在微分方程中, , 自变量的个数只需一个自变量的个数只需一个, , 称为常微称为常微分方程。分方程。自变量的个数为两个或两个以上的微分方程叫偏微自变量的个数为两个或两

4、个以上的微分方程叫偏微分方程。分方程。微分方程中出现的未知函数最高阶导数的阶数称为微分方程中出现的未知函数最高阶导数的阶数称为微分方程的阶数。微分方程的阶数。微分方程分类微分方程分类常微分方程：常微分方程：0)(),() 1 (yaybxayxfy )(,)(),()2(0ayyaybxayyxfy nybyyaybxayyxfy)(,)(),()3(02002212210012111)(),()(),()4(yxyyyxfyyxyyyxfy但能求解析解的常微分方程是有限的，大多数的常但能求解析解的常微分方程是有限的，大多数的常微分方程是给不出解析解的微分方程是给不出解析解的. . 22yxy

5、 这个一阶微分方程就不能用初等函数及其积这个一阶微分方程就不能用初等函数及其积分来表达它的解。分来表达它的解。 例例 例例 1)0(yyy,xey 的解的解的值仍需插值方法来计算的值仍需插值方法来计算. .xe222204,xxtdyxydxyee dt例 ：这是微积分的发明者之一Leibniz在1686年曾经让当时数学界人士求解的一阶微分方程式，吸引了许多数学家的注意，大约经过150年的探索，到1838年。刘维尔（Liouville)在理论上证明了这个微分方程不能用初等积分法求解，得借助于数值方法y =1-2xy例5：y(0)=0的解为为计算具体函数值，还需要数值积分方法。如果需要许多点处的

6、y值，则计算工作量较大。可见，用求解析的方法来找微分方程的数值解有时是不适宜的，必须研究微分方程的数值解法来求出其数值解现实上，从实践问题当中笼统出来的微分方程，通现实上，从实践问题当中笼统出来的微分方程，通常主要依托数值解法来处理。常主要依托数值解法来处理。可以证明可以证明: : 假设函数在带形区域假设函数在带形区域 R=axb, R=axb,-y y内延续，且关于内延续，且关于y y满足李普希兹满足李普希兹(Lipschitz)(Lipschitz)条件，即存在常数条件，即存在常数L(L(它与它与x,yx,y无关无关) )使使 2121),(),(yyLyxfyxf对对R R内恣意两个内恣

7、意两个 都成立都成立, ,那么方程的解那么方程的解 在在a, ba, b上存在且独一。上存在且独一。 21, yy)(xyy 00)(),(yxyyxfy在区间在区间a x ba x b上的数值解法。上的数值解法。 主要讨论一阶常微分方程初值问题主要讨论一阶常微分方程初值问题 微分方程数值方法的根本思想微分方程数值方法的根本思想 对常微分方程初值问题的数值解法，就对常微分方程初值问题的数值解法，就是要算出准确解是要算出准确解y(x)y(x)在区间在区间a,ba,b上的一系上的一系列离散节点列离散节点 处的函数值处的函数值bxxxxann110的近似 值)(,),(),(10nxyxyxynyy

8、y,10iixxh1niihxxi, 2 , 1 ,0相邻两个节点的间距相邻两个节点的间距 称为步长，步称为步长，步长可以相等，也可以不等。假定长可以相等，也可以不等。假定h h为定数，称为定步长，为定数，称为定步长，这时节点可表示为这时节点可表示为数值解法需求把延续性的问题加以离散化，从数值解法需求把延续性的问题加以离散化，从而求出离散节点的数值解。而求出离散节点的数值解。离散化离散化 对常微分方程数值解法的根本出发点就是离散化。对常微分方程数值解法的根本出发点就是离散化。其数值解法的根本特点：采用其数值解法的根本特点：采用“步进式，步进式，即求解过程顺着节点陈列的次序一步一步地向前推进，即

9、求解过程顺着节点陈列的次序一步一步地向前推进，描画这类算法，要求给出用知信息描画这类算法，要求给出用知信息 计算计算 的递推公式。的递推公式。建立这类递推公式的根本方法是在这些节点上用数建立这类递推公式的根本方法是在这些节点上用数值积分、数值微分、泰勒展开等离散化方法，对初值积分、数值微分、泰勒展开等离散化方法，对初值问题值问题021,yyyyiii1iy00)(),(yxyyxfyy中的导数中的导数 进展不同的离散化处置。进展不同的离散化处置。数值解和准确解数值解和准确解用数值方法求解初值问题，不是求出它的解析解或其近似解析式，而是给出它的解在某些离散节点上的近似值。用y(x)表示问题的准确

10、解y(x0), y(x1),y(xN) 表示解y(x)在节点x0, x1, xN处的准确值y0,y1,y N表示数值解，即问题的解y(x) 在相应节点处的近似值。单步法和多步法单步法和多步法单步法：在计算yi+1 时只利用y i多步法：在计算yi+1 时不仅利用y i , 还要利用 yi1, yi2,k步法：在计算yi+1 时要用到yi,yi1,yik+1显式格式可写成:yk+1=yk+hf(xk，yk;h隐式格式:yk+1=yk+hfxk，yk,yk+1;h它每步求解yk+1需求解一个隐式方程。欧拉欧拉EulerEuler方法方法在x= x0 处，用差商替代导数：hxyxyxxxyxyxy)

11、()()()()(0101010由由00000)(),()(yxyyxfxy10001),()(yyxhfyxy得得同理，在同理，在x= xn x= xn 处，用差商替代导数：处，用差商替代导数：hxyxyxxxyxyxynnnnnnn)()()()()(111),()()(1nnnnyxhfxyxy),()(nnnyxfxy由由得得假设记假设记11)(,)(nnnnyxyyxy那么上式可记为那么上式可记为),(1nnnnyxhfyy例例: : 用用EulerEuler方法求解常微分方程初值问题方法求解常微分方程初值问题 yyxyxy203002()( ).并将数值解和该问题的解析解比较。并将

12、数值解和该问题的解析解比较。21)(xxxy解析解：解：解：EulerEuler方法的详细格式：方法的详细格式：)2(21nnnnnyxyhyyh=0.2;y(1)=0.2;x=0.2:h:3;h=0.2;y(1)=0.2;x=0.2:h:3;for n=1:14for n=1:14 xn=x(n);yn=y(n); xn=x(n);yn=y(n); y(n+1)=yn+h y(n+1)=yn+h* *(yn/xn-2(yn/xn-2* *ynyn* *yn);yn);endendx0=0.2:h:3;y0=x0./(1+x0.2);x0=0.2:h:3;y0=x0./(1+x0.2);plo

13、t(x0,y0,x,y,x,y,o)plot(x0,y0,x,y,x,y,o)程序实现程序实现xn y(xn) yn yn-y(xn)0.00000.20.19230.20000.00770.40.34480.38400.03920.60.44120.51700.07580.80.48780.58240.09461.00.50000.59240.09241.20.49180.57050.07871.40.47300.53540.0624h = 0 . 2 , x n = n h , ( n = 0 , 1 , 2 , 1 5 ) , h = 0 . 2 , x n = n h , ( n =

14、0 , 1 , 2 , 1 5 ) , f(x,y)=y/x 2y2 f(x,y)=y/x 2y2 计算中取计算中取f(0,0)=1. f(0,0)=1. 计算计算结果如下：结果如下：xn y(xn) yn yn-y(xn)1.60.44940.49720.04781.80.42450.46050.03592.00.40000.42680.02682.20.37670.39660.01992.40.35500.36980.01472.60.33510.34590.01082.80.31670.32460.00793.00.30000.30570.0057由表中数据可以看到，微分方程初值问题的数

15、值解和解由表中数据可以看到，微分方程初值问题的数值解和解析解的误差普通在小数点后第二位或第三位小数上，阐析解的误差普通在小数点后第二位或第三位小数上，阐明明EulerEuler方法的精度是比较差的。方法的精度是比较差的。二阶二阶Runge-KuttaRunge-Kutta方法方法),(),()(12122122111hKyhxfKyxfKKcKchyynnnnnn),(),()(2122111nnnnnnnnhfyhxfcyxfchyxyT即即021021012122221cccc212112122221cccc例例 蛇形曲线的初值问题蛇形曲线的初值问题令令f(x,y)=y/x 2y2, f(

16、x,y)=y/x 2y2, 取取 f(0,0)=1, f(0,0)=1, h=0.2,xn=hn , ( n = 1,2,15)h=0.2,xn=hn , ( n = 1,2,15)2 2阶龙格阶龙格- -库塔公式计算格式：库塔公式计算格式： k1=yn/xn 2yn 2, k1=yn/xn 2yn 2, k2 = (yn+hk1)/(xn+h) k2 = (yn+hk1)/(xn+h) 2(yn+hk1)22(yn+hk1)2 yn+1=yn + 0.5h k1 + k2 yn+1=yn + 0.5h k1 + k222( 03 )(0)0.yyyxxy 00.511.522.530.10.

17、0.6x0=0;y0=0;h=.2;x=.2:h:3;k1=1;k2=(y0+h*k1)/x(1)-2*(y0+h*k1)2;y(1)=y0+.5*h*(k1+k2);for n=1:14 k1=y(n)/x(n)-2*y(n)2; k2=(y(n)+h*k1)/x(n+1)-2*(y(n)+h*k1)2; y(n+1)=y(n)+0.5*h*(k1+k2);endy1=x./(1+x.2);plot(x,y,o,x,y1),()2,2()2,2(),()22(6342312143211hKyhxfKKhyhxfKKhyhxfKyxfKKKKKhyynnnnnnnnnn四阶

18、四阶Runge-KuttaRunge-Kutta方法方法)(txx )()()(txtxty)()()()()(12tytytxtxty function ydot = harmonic(t,y) ydot=y(2); -y(1) y=inline(0 1;-1 0*y,t,y);System of Equations33)(/ )()()(/ )()(trtvtvtrtutu 22)()()(tvtutr)()()()()(tvtutvtuty33)(/ )()(/ )()()()(trtvtrtutvtuty function ydot =twobody(t,y) r=sqrt(y(1)2

19、+y(2)2); ydot=y(3);y(4);-y(1)/r3;-y(2)/r3; Two Body Problem)()(),(),(ccccyyJttytfytf),(),(ccccytyfJyttf),.,(),.,(),.,()()()(1121121nnnnnyytfyytfyytftytytydtdnnnnnnyfyfyfyfyfyfyfyfyfJ212221212111Jyy Linearized Differential EquationsJ的特征值是kkki)(diagk1VVJyVx kkkxx)()()(cttktxetxck0k0k0k解增长解衰减解振荡 基于龙格库塔

20、法，基于龙格库塔法， MATLABMATLAB求常微分方程数求常微分方程数值解的函数，普通调用格式为：值解的函数，普通调用格式为： t,y=ode23(fname,tspan,y0)t,y=ode23(fname,tspan,y0) t,y=ode45(fname,tspan,y0) t,y=ode45(fname,tspan,y0)其中其中fnamefname是定义是定义f(t,y)f(t,y)的函数文件名，该函数文的函数文件名，该函数文件必需前往一个列向量。件必需前往一个列向量。tspantspan方式为方式为t0,tf,t0,tf,表表示求解区间。示求解区间。y0y0是初始形状列向量。是

21、初始形状列向量。t t和和y y分别给分别给出时间向量和相应的形状向量。出时间向量和相应的形状向量。MATLABMATLAB求常微分方程数值解的函数求常微分方程数值解的函数 ode23: Bogacki,Shampine(1989)和Shampine(1994), 23表示用两算法：一个2阶，一个3阶)9865(72),()432(9)43,43()2,2(),(432111143211123121ssssheytfsssshyyhtthsyhtfsshyhtfsytfsnnnnnnnnnnnnnBogacki, P. and LF Shampine, A 3(2) pair of Runge

22、-Kutta formulas, Appl.Math. Letters, Vol. 2, 1989, pp 1-9. BS23 algorithm F=inline(y(2);-y(1),t,y) ode23(F,0 2*pi,1;0) opts=odeset(reltol,1.e-4,abstol,1.e-6,outputfcn)Examples ode23(twobody,0 2*pi,1;0;0;1);01234567-1.5-1-0.500.511.5Examples y0=1;0;0;3; ode23(twobody,0 2*pi,y0); 01234567-202468101214

23、1618 y0=1;0;0;3; t,y=ode23(twobody,0 2*pi,y0); plot(y(:,1),y(:,2); axis equal-1.5-1-0.500.51024681012141618 y0=1;0;0;3; t,y=ode23(twobody,0 2*pi,y0); plot(y(:,1) plot(y(:,2) A problem is stiff if the solution being sought is varying slowly, but there are nearby solutions that very rapidly, so the nu

24、merical method must take small steps to obtain satisfactory results.A model of flame propagation火焰熄灭:/20)0(32tyyyy y是球的半径，y2和y3与球的外表积和体积有关想一下，点燃一根火柴，光球迅速增长，到达一个关键的大小，然后维持它的大小由于进入球内氧气和耗费的氧气平衡Stiff Problem(刚性问题)010203040506070809010000.811.21.4eta=0.02;sym y; F=inline(y2-y3,t,y);ode23(F,0 2/e

25、ta,eta);eta=0.00002;sym y; F=inline(y2-y3,t,y);ode23(F,0 2/eta,eta);012345678910 x 10400.811.21.4ode23seta=0.00002; ode23s(inline(y2-y3,t,y),0 2/eta,eta);例例 蛇形曲线的常微分方程初值问题蛇形曲线的常微分方程初值问题 22211yxy 0)0( yMATLABMATLAB数值求解命令数值求解命令F=inline(1./(1+x.2)-2F=inline(1./(1+x.2)-2* *y.2)y.2)；ode23(F,0,6,

26、0)ode23(F,0,6,0)输出结果为图形输出结果为图形 T,y=ode23(f,0,6,0)T,y=ode23(f,0,6,0)将将得到自变量和函数的离散数得到自变量和函数的离散数据据 T =0 0.0001 0.0005 0.0025 0.0125 0.0496 0.1085 0.1863 0.2837 0.4091 0.5991 0.8513 1.0567 1.2680 1.5110 1.8050 2.1788 2.6842 3.2842 3.8842 4.4842 5.0842 5.6842 6.0000y =0 0.0001 0.0005 0.0025 0.0125 0.0495

27、 0.1073 0.1800 0.2626 0.3505 0.4411 0.4944 0.5000 0.4868 0.4607 0.4242 0.3793 0.3270 0.2783 0.2411 0.2122 0.1891 0.1705 0.1620例例 洛伦兹模型由如下常微分方程组描画洛伦兹模型由如下常微分方程组描画 zyxydtdzzydtdyyzxdtdx )(取取 =8/3=8/3， =10=10，=28=28。初值：初值：x(0)=0 x(0)=0，y(0)=0y(0)=0，z(0)=0.01z(0)=0.01。利用利用MATLABMATLAB求解常微分方程数值解命令计算出求解常微

28、分方程数值解命令计算出t0t0，8080内，三个未知函数的数据值，并绘出相空间在内，三个未知函数的数据值，并绘出相空间在Y-XY-X平面的投影曲线平面的投影曲线 气候学家Lorenz提出一篇论文，名叫一只蝴蝶拍一下翅膀会不会在Taxas州引起龙卷风？论述某系统假设初期条件差一点点，结果会很不稳定，他把这种景象戏称做蝴蝶效应。Lorenz为何要写这篇论文呢？这故事发生在1961年的某个冬天，他如往常普通在办公室操作气候电脑。平常，他只需求将温度、湿度、压力等气候数据输入，电脑就会根据三个内建的微分方程式，计算出下一刻能够的气候数据，因此模拟出气候变化图。这一天，Lorenz想更进一步了解某段纪录

29、的後续变化，他把某时辰的气候数据重新输入电脑，让电脑计算出更多的後续结果。当时，电脑处置数据资料的数度不快，在结果出来之前，足够他喝杯咖啡并和友人闲聊一阵。在一小时後，结果出来了，不过令他目瞪口呆。结果和原资讯两相比较，初期数据还差不多，越到後期，数据差别就越大了，就像是不同的两笔资讯。而问题并不出在电脑，问题是他输入的数据差了0.000127，而这些微的差别却呵斥天壤之别。所以长期的准确预测天气是不能够的。天气预告的准确性：nani.tw/senior/seniorteach/stfastnews/stnfastnews/stnfastnews4/stnfastnews4_2.htmLore

30、nz景象的数学： /news_detail.php?news_id=225分形艺术电子版：/personal/huajie/fractalart/html/notice.htm混沌实际：/yzhao/doc/chaos.html记向量记向量 y1y1，y2y2，y3 = xy3 = x，y y，zz，创建，创建MATLABMATLAB函函数文件如下数文件如下 function z=flo(t,y)function z=flo(t,y)z(1,:)=-8z(1,:)=-8* *y(1)/3+y(2).y(1)/3+y(2).*

31、*y(3);y(3);z(2,:)=-10z(2,:)=-10* *(y(2)-y(3);(y(2)-y(3);z(3,:)=-y(1).z(3,:)=-y(1).* *y(2)+28y(2)+28* *y(2)-y(3);y(2)-y(3);用用MATLABMATLAB命令求解并绘出命令求解并绘出Y-XY-X平面的投影图平面的投影图 y0=0;0;0.01;y0=0;0;0.01;x,y=ode23(flo,0, 80,y0)x,y=ode23(flo,0, 80,y0)plot(y(:,2),y(:,1) plot(y(:,2),y(:,1) -20-15-10-5051015200102

32、0304050plot(y(:,3),y(:,1) plot(y(:,3),y(:,1) plot3(y(:,1),y(:,2),y(:,3) y0=30;0;-40; plot(y(:,i)非刚性系统： ode45(Runge-Kutta45) ode23(Runge-Kuatta23) ode113(Adams-Bashforth-Moulton PECE)多步方法刚性系统： ode15s(Gear方法) ,多步方法 ode23s(二阶modified Rosenbroack formula)，单步 ode23t(trapezoidal rule),solve DAEs ode23tb(T

33、R-BDF2) low order methodMatlabs ODE Solvers2222yxLaplacian 算子： Poisson方程elliptic：fu Laplacian 算子的特征值问题:fuuHeat equation(parabolic):utuWave equation(hyperbolic):utu22PDE Model五点离散22),(),(2),( ),(),(2),(),(hhyxuyxuhyxuhyhxuyxuyhxuyxuh2)(4)()()()()(hPuSuEuWuNuPuh0)(Puh)()(PfPuh)()(PuPukkh Poisson方程离散：

34、特征值问题:Finete Difference Methods),(),(),(yxutyxutyxuh),(),(),(yxutyxutyxuh),(),(),(2),(2yxutyxutyxutyxuh),(),(),(2),(2yxutyxutyxutyxuh热方程：动摇方程：-4 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 -4 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 -4 1 1 0 0 0 0 0 0 0 0 0 0 0 0 1 1 -4 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 -4 1 1 0 0 0 0 0 0 0 0 0 0 0 0 1 1 -4 0 1 0 0 0 0 0 0 0 0 0