MATLAB解微分方程

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

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(xD(x).).x1o40)0()(412VxDdxdV只要求解上述方程只要求解上述方程, ,就可求出

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

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

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

6、果需要许多点处的y值，则计算工作量较大。可见，用求解析的方法来找微分方程的数值解有时是不适宜的，必须研究微分方程的数值解法来求出其数值解9事实上，从事实上，从实际问题实际问题当中抽象出来的微分方程，通当中抽象出来的微分方程，通常主要依靠常主要依靠数值解法数值解法来解决。来解决。可以证明可以证明: : 如果函数在带形区域如果函数在带形区域 R=axb,R=axb,-y y内连续，且关于内连续，且关于y y满足李普希兹满足李普希兹(Lipschitz(Lipschitz) )条件，即存在常数条件，即存在常数L(L(它与它与x,yx,y无关无关) )使使 2121),(),(yyLyxfyxf对对R

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

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

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

10、上的近似值。用y(x)表示问题的准确解y(x0), y(x1),y(xN) 表示解y(x)在节点x0, x1, xN处的准确值y0,y1,y N表示数值解，即问题的解y(x) 在相应节点处的近似值。14单步法和多步法单步法和多步法单步法：在计算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+hf（xk，yk,yk+1;h）它每步求解yk+1需要解一个隐式方程。15欧拉（欧拉（EulerEuler）方法）方法在x

11、= x0 处，用差商代替导数：hxyxyxxxyxyxy)()()()()(0101010由由00000)(),()(yxyyxfxy10001),()(yyxhfyxy得得16同理，在同理，在x= xx= xn n 处，用差商代替导数：处，用差商代替导数：hxyxyxxxyxyxynnnnnnn)()()()()(111),()()(1nnnnyxhfxyxy),()(nnnyxfxy由由得得若记若记11)(,)(nnnnyxyyxy则上式可记为则上式可记为),(1nnnnyxhfyy17例例: : 用用EulerEuler方法求解常微分方程初值问题方法求解常微分方程初值问题 yyxyxy2

12、03002()( ).并将数值解和该问题的解析解比较。并将数值解和该问题的解析解比较。21)(xxxy解析解：解：解：EulerEuler方法的具体格式：方法的具体格式：)2(21nnnnnyxyhyy18h=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 xn=x(n);yn=y(n);=y(n); y(n+1)=yn+h y(n+1)=yn+h* *(yn/xn-2(yn/xn-2* *ynyn* *ynyn););endendx0=0.2:h:3;y0=x0./(1+x0.2)

13、;x0=0.2:h:3;y0=x0./(1+x0.2);plot(x0,y0,x,y,x,y,o)plot(x0,y0,x,y,x,y,o)程序实现程序实现19xn 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 h=0.2, =0.2, x xn n= =nhnh,(,(n n=0,1,2=0,1,2

14、,15), ,15), f f( (x x, ,y y)=)=y y/ /x x 2 2y y2 2 计算中取计算中取f f(0,0)=1. (0,0)=1. 计算结果如下：计算结果如下：20 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方法的精度是比较差的。方法的精度是比较差的。21二阶二阶Runge-KuttaRunge-Kutta方法方法),(),()(12122122111hKyhxfKyxfKKcKchyynnnnnn),(),()(2122111nnnnnnnnhfyhxfcyxfchyxyT22即即021021012122221cccc212112122221cccc23例例 蛇形曲线蛇形曲线的初值问题的初值问题令令f(x,y)=

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

17、 )(0)0.yyyxxy 2400.511.522.5x0=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)25),()2,2()2,2(),()22(6342312143211hKyhxfK

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

19、 function ydot =twobody(t,y) r=sqrt(y(1)2+y(2)2); ydot=y(3);y(4);-y(1)/r3;-y(2)/r3; Two Body Problem28)()(),(),(ccccyyJttytfytf),(),(ccccytyfJyttf),.,(),.,(),.,()()()(1121121nnnnnyytfyytfyytftytytydtdnnnnnnyfyfyfyfyfyfyfyfyfJ212221212111Jyy Linearized Differential Equations29J的特征值是kkki)(diagk1VVJyVx

20、 kkkxx)()()(cttktxetxck0k0k0k解增长解衰减解振荡30 基于龙格库塔法，基于龙格库塔法， MATLABMATLAB求常微分方程数值解求常微分方程数值解的函数的函数，一般调用格式为：一般调用格式为： t,yt,y=ode23(fname,tspan,y0)=ode23(fname,tspan,y0) t,y t,y=ode45(fname,tspan,y0)=ode45(fname,tspan,y0)其中其中fnamefname是定义是定义f(t,yf(t,y) )的函数文件名的函数文件名，该函数文，该函数文件必须返回一个列向量。件必须返回一个列向量。tspantspa

21、n形式为形式为t0,tf,t0,tf,表表示求解区间。示求解区间。y0y0是初始状态列向量。是初始状态列向量。t t和和y y分别给分别给出时间向量和相应的状态向量。出时间向量和相应的状态向量。MATLABMATLAB求常微分方程数值解的函数求常微分方程数值解的函数31 ode23: Bogacki,Shampine(1989)和Shampine(1994), ”23”表示用两算法：一个2阶，一个3阶)9865(72),()432(9)43,43()2,2(),(432111143211123121ssssheytfsssshyyhtthsyhtfsshyhtfsytfsnnnnnnnnnnn

22、nnBogacki, P. and LF Shampine, A 3(2) pair of Runge-Kutta formulas, Appl.Math. Letters, Vol. 2, 1989, pp 1-9. BS23 algorithm32 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)Examples33 ode23(twobody,0 2*pi,1;0;0;1);01234567-1.5-1-0.500.511.5Examples34 y0

23、=1;0;0;3; ode23(twobody,0 2*pi,y0); 01234567-202468101214161835 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.5102468101214161836 y0=1;0;0;3; t,y=ode23(twobody,0 2*pi,y0); plot(y(:,1) plot(y(:,2)37 A problem is stiff if the solution being sought is varying slow

24、ly, but there are nearby solutions that very rapidly, so the numerical method must take small steps to obtain satisfactory results.A model of flame propagation（火焰燃烧）:/20)0(32tyyyy y是球的半径，y2和y3与球的表面积和体积有关想一下，点燃一根火柴，光球迅速增长，到达一个关键的大小，然后维持它的大小（由于进入球内氧气和消耗的氧气平衡）Stiff Problem(刚性问题)380102030405060708090100

25、00.811.21.4eta=0.02;sym y; F=inline(y2-y3,t,y);ode23(F,0 2/eta,eta);39eta=0.00002;sym y; F=inline(y2-y3,t,y);ode23(F,0 2/eta,eta);40012345678910 x 10400.811.21.4ode23seta=0.00002; ode23s(inline(y2-y3,t,y),0 2/eta,eta);41例例 蛇形曲线的常微分方程初值问题蛇形曲线的常微分方程初值问题 22211yxy 0)0( yMATLABMATLAB数值求

26、解命令数值求解命令F=inline(1./(1+x.2)-2F=inline(1./(1+x.2)-2* *y.2)y.2)；ode23(F,0,6,0)ode23(F,0,6,0)输出结果为图形输出结果为图形 42T,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.2

27、842 3.8842 4.4842 5.0842 5.6842 6.0000y =0 0.0001 0.0005 0.0025 0.0125 0.0495 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.162043例例 洛伦兹模型洛伦兹模型由如下常微分方程组描述由如下常微分方程组描述 zyxydtdzzydtdyyzxdtdx )(取取 = =8/38/3， = =1010， = =2828。初值初值

28、：x x(0)=0(0)=0，y y(0)=0(0)=0，z z(0)=0.01(0)=0.01。利用利用MATLABMATLAB求解常微分方程数值解命令计算出求解常微分方程数值解命令计算出t t00，8080内内，三个未知函数的数据值，并绘出相空间在三个未知函数的数据值，并绘出相空间在Y-XY-X平面的投影曲线平面的投影曲线 44气象学家Lorenz提出一篇论文，名叫一只蝴蝶拍一下翅膀会不会在Taxas州引起龙卷风？论述某系统如果初期条件差一点点，结果会很不稳定，他把这种现象戏称做蝴蝶效应蝴蝶效应。Lorenz为何要写这篇论文呢？这故事发生在1961年的某个冬天，他如往常一般在办公室操作气象

29、电脑。平时，他只需要将温度、湿度、压力等气象数据输入，电脑就会依据三个内建的微分方程式，计算出下一刻可能的气象数据，因此模拟出气象变化图。这一天，Lorenz想更进一步了解某段纪录的後续变化，他把某时刻的气象数据重新输入电脑，让电脑计算出更多的後续结果。当时，电脑处理数据资料的数度不快，在结果出来之前，足够他喝杯咖啡并和友人闲聊一阵。在一小时後，结果出来了，不过令他目瞪口呆。结果和原资讯两相比较，初期数据还差不多，越到後期，数据差异就越大了，就像是不同的两笔资讯。而问题并不出在电脑，问题是他输入的数据差了0.000127，而这些微的差异却造成天壤之别。所以长期的准确预测天气是不可能的。45天气

30、预报的准确性：http:/.tw/senior/seniorteach/stfastnews/stnfastnews/stnfastnews4/stnfastnews4_2.htmLorenz现象的数学： /news_detail.php?news_id=225分形艺术电子版：http:/ y y1 1，y y2 2，y y3 3 = = x x，y y，z z ，创建创建MATLABMATLAB函函数文件如下数文件如下 function z=flo(t,yfunction z=flo(t,y) )z(1,:)=-8z(1,:)=-8* *y(1)/

31、3+y(2).y(1)/3+y(2).* *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) 47

32、-20-15-10-50510152001020304050plot(y(:,3),y(:,1) plot(y(:,3),y(:,1) 48plot3(y(:,1),y(:,2),y(:,3)49 y0=30;0;-40; plot(y(:,iy0=30;0;-40; plot(y(:,i)505152非刚性系统： ode45(Runge-Kutta45) ode23(Runge-Kuatta23) ode113(Adams-Bashforth-Moulton PECE)多步方法刚性系统： ode15s(Gear方法) ,多步方法 ode23s(二阶modified Rosenbroack f

33、ormula)，单步 ode23t(trapezoidal rule),solve DAEs ode23tb(TR-BDF2) low order methodMatlabs ODE Solvers532222yxLaplacian 算子： Poisson方程（elliptic）：fu Laplacian 算子的特征值问题:fuuHeat equation(parabolic):utuWave equation(hyperbolic):utu22PDE Model54五点离散22),(),(2),( ),(),(2),(),(hhyxuyxuhyxuhyhxuyxuyhxuyxuh2)(4)(

34、)()()()(hPuSuEuWuNuPuh0)(Puh)()(PfPuh)()(PuPukkh Poisson方程离散： 特征值问题:Finete Difference Methods55),(),(),(yxutyxutyxuh),(),(),(yxutyxutyxuh),(),(),(2),(2yxutyxutyxutyxuh),(),(),(2),(2yxutyxutyxutyxuh热方程：波动方程：56-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 0 0 0 1 0 -4