二階常微分方程的數(shù)值求解優(yōu)秀課件

1、1二階常微分方程的數(shù)值求解二階常微分方程的數(shù)值求解一一. 教學(xué)要求教學(xué)要求 掌握利用降階把二階常微分方程轉(zhuǎn)化為一階微分方程組，再利用Euler方法數(shù)值求解,并能利用MATLAB軟件進(jìn)行數(shù)值計(jì)算和符號(hào)運(yùn)算。二二. 教學(xué)過(guò)程教學(xué)過(guò)程2q 考慮如下的二階微分方程考慮如下的二階微分方程初值問(wèn)題初值問(wèn)題2012( , ,) , ( ) , ( ), , d yf x y yy ayy ayxa bdx 100( )( ), , ( )( , ( ), ( ), , ( ), ( )yxz xxa bzxf x y x z xxa bz ayzy ay 3q利用利用Euler方法求解上述方程組可得如下數(shù)方

2、法求解上述方程組可得如下數(shù)值格式值格式001111 2( ),( ),(,),.kkkkkkkkkky ayy azyyhzzzhf xyzkxxh 4q利用四階利用四階R-K方法求解上述方程組可得如下方法求解上述方程組可得如下數(shù)值格式數(shù)值格式11234112341121211323224343322622622222222(),(),(,),(,),(,),(,).kkkkkkkkkkkkkkkkkkkkhyyKKKKhzzLLLLKzLf xyzhhhhKzLLf xyK zLhhhhKzLLf xyKzLkKzhLLf xh yhKzhL 1 2, 5例例1：用用 Euler 法求解如下

3、初值問(wèn)題法求解如下初值問(wèn)題220 20101 , ( ),( )d yyxdxyy 當(dāng)當(dāng) h=0.1，即，即 n=20 時(shí)，時(shí)，Matlab 源程序見(jiàn)源程序見(jiàn) Euler_sys1.m解：解：00 20 20101,( )( ), , ( )( ), , ( ), ( ).令則該初值問(wèn)題可以轉(zhuǎn)化為zyyxz xxzxy xxzzy 6clc;clear;h=0.1;a=0;b=2;x=a:h:b;y(1)=1;z(1)=-1;for i=1:length(x)-1 y(i+1)=y(i)+h*z(i); z(i+1)=z(i)+h*y(i);endplot(x,y,r+,x,exp(-x),k

4、-);xlabel(Variable x); ylabel(Variable y); Euler_sys1.m7數(shù)值解與真解如下圖數(shù)值解與真解如下圖8例例2：利用利用4階階R-K方法求解例方法求解例1，并與，并與Euler方法方法進(jìn)行比較。進(jìn)行比較。解解 當(dāng)當(dāng) h=0.1，即，即 n=20 時(shí)，時(shí)，R-K方法的方法的Matlab 源程序見(jiàn)源程序見(jiàn) RK_sys1.m，數(shù)值結(jié)果見(jiàn)下圖，數(shù)值結(jié)果見(jiàn)下圖9function w=rightf_sys1(x,y,z) w=y;clc;clear;h=0.1;a=0;b=2;x=a:h:b;Euler_y(1)=1; Euler_z(1)=-1; %初值R

5、K_y(1)=1; RK_z(1)=-1; %初值for i=1:length(x)-1 %* Euler Method *% Euler_y(i+1)=Euler_y(i)+h*Euler_z(i); Euler_z(i+1)=Euler_z(i)+h*Euler_y(i); %* R-K4 Method*% K1=RK_z(i); L1=rightf_sys1(x(i),RK_y(i),RK_z(i); % K1 and L1 K2=RK_z(i)+0.5*h*L1; rightf_sys1.mRK_sys1.m10L2=rightf_sys1(x(i)+0.5*h,RK_y(i)+0.5

6、*h*K1,RK_z(i)+0.5*h*L1); % K2 and L2 K3=RK_z(i)+0.5*h*L2; L3=rightf_sys1(x(i)+0.5*h,RK_y(i)+0.5*h*K2,RK_z(i)+0.5*h*L2); % K3 and L3 K4=RK_z(i)+h*L3; L4=rightf_sys1(x(i)+h,RK_y(i)+h*K3,RK_z(i)+h*L3); % K4 and L4 RK_y(i+1)=RK_y(i)+1/6*h*(K1+2*K2+2*K3+K4); RK_z(i+1)=RK_z(i)+1/6*h*(L1+2*L2+2*L3+L4); end

7、plot(x,Euler_y,r+,x,exp(-x),k-,x,RK_y,b*);xlabel(Variable x); ylabel(Variable y); 11例例3：分別用分別用 Euler 法和法和R-K4求解如下初值問(wèn)題求解如下初值問(wèn)題22210 20101 () , ( ),( )xd yyexxdxyy 解：解：00 2210 20101,( )( ), , ( )( )(), , ( ), ( ).令則該初值問(wèn)題可以轉(zhuǎn)化為xzyyxz xxzxy xexxzzy 12當(dāng)當(dāng) h=0.1，即，即 n=20 時(shí)，時(shí)，Matlab 源程序見(jiàn)源程序見(jiàn) RK_sys2.m, 數(shù)值結(jié)數(shù)值

8、結(jié)果如下圖果如下圖13function w=rightf_sys2(x,y,z)w=-y+2*exp(-x)*(x-1);clc;clear;h=0.1;a=0;b=2;x=a:h:b;Euler_y(1)=1; Euler_z(1)=1; RK_y(1)=1; RK_z(1)=1; for i=1:length(x)-1%* Euler Method *% Euler_y(i+1)=Euler_y(i)+h*Euler_z(i); Euler_z(i+1)=Euler_z(i)+h*rightf_sys2(x(i),Euler_y(i),Euler_z(i); %* R-K4 Method*

9、% K1=RK_z(i); L1=rightf_sys2(x(i),RK_y(i),RK_z(i); % K1 and L1 rightf_sys1.mRK_sys2.m14K2=RK_z(i)+0.5*h*L1; L2=rightf_sys2(x(i)+0.5*h,RK_y(i)+0.5*h*K1,RK_z(i)+0.5*h*L1); % K2 and L2 K3=RK_z(i)+0.5*h*L2; L3=rightf_sys2(x(i)+0.5*h,RK_y(i)+0.5*h*K2,RK_z(i)+0.5*h*L2); % K3 and L3 K4=RK_z(i)+h*L3; L4=rig

10、htf_sys2(x(i)+h,RK_y(i)+h*K3,RK_z(i)+h*L3); % K4 and L4 RK_y(i+1)=RK_y(i)+1/6*h*(K1+2*K2+2*K3+K4); RK_z(i+1)=RK_z(i)+1/6*h*(L1+2*L2+2*L3+L4); endplot(x,Euler_y,r+,x,cos(x)+x.*exp(-x),k-,x,RK_y,b*);xlabel(Variable x); ylabel(Variable y); 15q dsolve 的調(diào)用格式的調(diào)用格式y(tǒng)=dsolve(eq1,eq2, . ,cond1,cond2, . ,v)其中其

11、中 y 為輸出，為輸出， eq1、eq2、.為微分方程，為微分方程，cond1、cond2、.為初值條件，為初值條件，v 為自變量，如果不指定為自變量，如果不指定v作為自變作為自變量，則默認(rèn)量，則默認(rèn)t為自變量。為自變量。例例 4：求微分方程求微分方程 的通解，并驗(yàn)證。的通解，并驗(yàn)證。22xdyxyxedx y=dsolve(Dy+2*x*y=x*exp(-x2),x) syms x; diff(y)+2*x*y - x*exp(-x2)q利用利用dsolvedsolve 函數(shù)求微分方程解析解函數(shù)求微分方程解析解16q 幾點(diǎn)說(shuō)明幾點(diǎn)說(shuō)明l 如果省略初值條件，則表示求通解；如果省略初值條件，則表

12、示求通解；l 如果省略自變量，則默認(rèn)自變量為如果省略自變量，則默認(rèn)自變量為 t dsolve(Dy=2*x,x); dy/dx = 2xdsolve(Dy=2*x); dy/dt = 2xl 若找不到解析解，則返回其積分形式。若找不到解析解，則返回其積分形式。l 微分方程中用微分方程中用 D 表示對(duì)表示對(duì) 自變量自變量 的導(dǎo)數(shù)，如：的導(dǎo)數(shù)，如：Dy y； D2y y； D3y y17例例 5：求微分方程求微分方程 在初值條件在初值條件 下的特解，并畫(huà)出解函數(shù)的圖形。下的特解，并畫(huà)出解函數(shù)的圖形。0 xxyye y=dsolve(x*Dy+y-exp(x)=0,y(1)=2*exp(1),x) ezplot(y);12( )ye 182220100cos()( ),( )求二階常微分方程的通解d yxydxyy 例例 6在在MatlabMatlab中的命令窗口中輸入下面的命令中的命令窗口中輸入下面的命令 syms x y S=dsolve(D2y=cos(2*x)-y,y(0)=1,Dy(0)=0,x)則可以得到如下的結(jié)果則可以得到如下的結(jié)果