《matlab第六章常微分方程的數(shù)值解》由會員分享,可在線閱讀,更多相關(guān)《matlab第六章常微分方程的數(shù)值解(9頁珍藏版)》請?jiān)谘b配圖網(wǎng)上搜索。
1、第六章 常微分方程的數(shù)值解 一、解法步驟:將高階方程轉(zhuǎn)化為一階方程組。(兩種情況)建立相應(yīng)的函數(shù)文件。調(diào)用求解函數(shù)。 t,z=odeij(dzdtk,H,z0,tol)剛性:設(shè)有一階常系數(shù)線性微分方程組y=Ay+f,如果它的Jacobian矩陣的特征值相差十分懸殊。 odeij問題類型精度適用對象ode45非剛性中等多數(shù)情況下可優(yōu)先選用,但不能用來解剛性問題ode23非剛性較低可用來解中等剛性問題,或誤差允許范圍較寬的問題ode113非剛性低到高不能解剛性問題,當(dāng)誤差容限要求嚴(yán)格時(shí)效果較ode45好ode15s剛性低到高可用于解剛性問題,當(dāng)采用ode45失敗或效果很差時(shí),可考慮使用ode23s
2、剛性低可用于解剛性問題,當(dāng)誤差容限較寬時(shí)效果比ode15s好ode23t適度剛性低可用于解剛性問題,但要求無數(shù)值衰減 ode23tb剛性低可用于解剛性問題,當(dāng)誤差容限較寬時(shí)效果比ode15s好 二、例題例1 在區(qū)間H=0.1,30上滿足初值條件x=0.1時(shí), 的特解(1)轉(zhuǎn)化原方程為一階方程組。(2)建立描述上面微分方程組的函數(shù)文件,并設(shè)文件名為dzdt01.m。function dz=f(t,z)dz(1)=z(2);dz(2)=z(3);dz(3)=(1/t)*z(3)-(2/t2)*z(2)+(2/t3)*z(1)+10*t3*sin(t);dz=dz(1);dz(2);dz(3) xx
3、yyxyxyx sin1022 623 0,0,0 yyy (3)調(diào)用求解。H=0.1,30;z0=0;0;0;t,z=ode45(dzdt01,H,z0);plot(t,z(:,1),m-) plot(t,z(:,1),m-) 0 5 10 15 20 25 30-2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5x 10 5 例2 H=0,50(1)轉(zhuǎn)化方程。(2)建立文件,設(shè)文件名為dzdt02.m。function dz=f(t,z)c=1;dz(1)=z(2);dz(2)=c*(1-z(1)2)*z(2)-z(1);dz=dz(1);dz(2);(3)調(diào)用求解。H=0,5
4、0;z0=2;1; t,z=ode45(dzdt02,H,z0);plot(t,z(:,1),m-,t,z(:,2),g-.) 0 5 10 15 20 25 30 35 40 45 50-3 -2 -1 0 1 2 3 0)1( 2 yyycy 1,2,0 yyt 例3 H=0,100(1)轉(zhuǎn)化方程。(2)建立文件,設(shè)文件名為dzdt03.m。function dz=f(t,z)a=-0.2;dz(1)=z(2);dz(2)=a*z(2)-sin(t);dz=dz(1);dz(2);(3)調(diào)用求解。H=0,100;z0=0;2; t,z=ode45(dzdt03,H,z0);plot(z(:
5、,1),z(:,2),m-) 0 1 2 3 4 5 6-1 -0.5 0 0.5 1 1.5 2 tyyyx sin2.0, 2,0,0 yxt 例4 H=-30,30(1)轉(zhuǎn)化方程。(2)建立文件,設(shè)文件名為dzdt04.m。function dz=f(t,z) dz(1)=z(2);dz(2)=-2*t*z(4);dz(3)=z(4);dz(4)=2*t*z(2);dz=dz(1);dz(2);dz(3);dz(4);(3)調(diào)用求解。H=-30,30;z0=0;1;0;0; t,z=ode45(dzdt04,H,z0);plot(z(:,1),z(:,3),m-) -0.4 -0.2 0
6、 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6-2 -1.5 -1 -0.5 0 0.5 02 02 xty ytx 0,1,0,0,0 0000 yxyxt 例5 (1)轉(zhuǎn)化方程。(2)建立文件,設(shè)文件名為dzdt05.m。function dz=f(t,z) a=8/3;b=10;c=28;dz(1)=-a*z(1)+z(2)*z(3);dz(2)=-b*(z(2)-z(3);dz(3)=c*z(2)-z(3)-z(1)*z(2);dz=dz(1);dz(2);dz(3);(3)調(diào)用求解。H=0,80;z0=0;0;0.222/1015; t,z=ode23(dzdt05,H,z0);plot3(z(:,1),z(:,2),z(:,3),m-) -20 0 20 40 60-20-1001020-30-20-100102030 xyzcyzzybyyzaxx ),(,