科學(xué)工程計(jì)算與matlab編程5-2

1、5.2 微分方程問題的數(shù)值解法5.2.1 算法概述0.( , ), ( )dxf t xatbdtx ax 微分方程求解的誤差與步長問題：000121( )1,2,.,:nkkkkatx tattttbyknhtth時(shí)刻系統(tǒng)狀態(tài)向量的值為數(shù)值解法：即求微分方程初值問題的解在若干點(diǎn) 的近似值 （）的方法。步長（一般取等步長 ）應(yīng)用差商近似導(dǎo)數(shù)舍入誤差10000( ,)xxhf txR0( )x t提高數(shù)值解精度的方法之一： 減小步長h的值； function outx,outy=MyEuler(fun,x0,xt,y0,PointNum) % fun 表示f(x,y); x0,xt:自變量的初值

2、和終值; %y0:函數(shù)在x0處的值,其可以為向量形式; %PointNum表示自變量在x0,xt上取的點(diǎn)數(shù)if nargin=4 | PointNumf1=(x,y)sin(x)+y; 2*pi/0.4 % 計(jì)算 xt=2pi 時(shí)應(yīng)取得點(diǎn)數(shù)ans = 15.7080 x1,y1=MyEuler(f1,0,2*pi,1,16); %h=0.4 歐拉法所的的解x11,y11=MyEuler(f1,0,2*pi,1,32); % h=0.2 歐拉法所的的解h=0.2和h=0.4的數(shù)值解。y=dsolve(Dy=y+sin(t),y(0)=1); %該常微分方程的解析解for k=1:33t(k)=x

3、11(k);y2(k)=subs(y,t(k); %求其對(duì)應(yīng)點(diǎn)的離散解endplot(x1,y1,+b,x11,y11,og,x11,y2,*r)legend(h=0.4的歐拉法解,h=0.2的歐拉解,符號(hào)解)誤差判斷：改進(jìn)的Euler算法 ( , )dxf t xdt1 kkttdt11()( )kkkkx tx txx1 kkttdt11( ,)(,)2kkkkhf txf tx1111( ,)( ,)(,)2kkkkkkkkkkxxhf txhxxf txf tx改進(jìn)Euler迭代公式Euler公式梯形公式function Xout,Yout=MyEulerPro(fun,x0,xt,y

4、0,PointNumber) %MyEulerPro 用改進(jìn)的歐拉法解微分方程if nargin=4 | PointNumber5.2.2 四階定步長Runge-Kutta算法及Matlab實(shí)現(xiàn)Euler算法：一階Taylor多項(xiàng)式近似函數(shù)基本思想：2()O h誤差：RK一般算法：推廣：高級(jí)Taylor多項(xiàng)式近似函數(shù)提高數(shù)值解精度11111( ,)(,) (2,3,4,.)nniininijjpjnnjjkf txkf ta h xb kixxhkp為展開的項(xiàng)數(shù)例如：p4，即Matlab 實(shí)現(xiàn)：實(shí)現(xiàn)： function tout,yout=rk_4(odefile,tspan,y0) y0初值

5、列向量 t0=tspan(1); th=tspan(2); if length(tspan) comet3(x(:,1),x(:,2),x(:,3) 描述微分方程是常微分方程初值問題數(shù)值求解的關(guān)鍵。 f1=inline(-8/3*x(1)+x(2)*x(3); -10*x(2)+10*x(3);,. -x(1)*x(2)+28*x(2)-x(3),t,x); t_final=100; x0=0;0;1e-10; t,x=ode45(f1,0,t_final,x0); plot(t,x), figure; plot3(x(:,1),x(:,2),x(:,3); axis(10 42 -20 20

6、 -20 25);得出完全一致的結(jié)果。5.2.3.3 MATLAB 下帶有附加參數(shù)的微分方程求解 例： 編寫函數(shù)function xdot=lorenz1(t,x,flag,beta,rho,sigma) % flag變量是不能省略的 xdot=-beta*x(1)+x(2)*x(3); -rho*x(2)+rho*x(3); -x(1)*x(2)+sigma*x(2)-x(3);3/8,10,28 t_final=100; x0=0;0;1e-10; b2=3/8; r2=10; s2=28;%變量名不必一定為 等 t2,x2=ode45(lorenz1,0,t_final,x0,b2,r2

7、,s2); plot(t2,x2), %options位置為，表示不需修改控制選項(xiàng) figure; plot3(x2(:,1),x2(:,2),x2(:,3); axis(10 42 -20 20 -20 25);例：求解時(shí)的Lorenz微分方程 t_final=100; x0=0;0;1e-10; b2=2; r2=5; s2=20;%變量名不必一定為 等 t2,x2=ode45(lorenz1,0,t_final,x0,b2,r2,s2); plot(t2,x2), %options位置為，表示不需修改控制選項(xiàng) figure; plot3(x2(:,1),x2(:,2),x2(:,3); axis(0 72 -20 22 -35 40);f2=inline(-beta*x(1)+x(2)*x(3); -rho*x(2)+rho*x(3);,. -x(1)*x(2)+sigma*x(2)-x(3), t,x,flag,beta,rho,sigma); flag變量是不能省略t_final=100; x0=0;0;1e-10; b2=2; r2=5; s2=20;%變量名不必一定為 等 t2,x2=ode45(f2,0,t_final,x0,