常微分方程數(shù)值解的若干Matlab函數(shù)文件實用教案

1、 function x,y=euler(f,tspan,y0,n) % 解初值問題： =f(x,y),y(a)=y0。 % 使用 n 步的 Euler 法，步長 h=(b-a)/h。 a=tspan(1);b=tspan(2);h=(b-a)/h; x=(a+h:h:b); y(1)=y0+h*feval(f,a,y0); for i=2:n y(i)=y(i-1)+h*fevol(f,x(i-1),y(i-1); end x=a x; y=y0 y;y 8.7常微分方程數(shù)值(shz)解的若干Matlab函數(shù)文件8.7.1 Euler 方法(fngf)的 Matlab 函數(shù)文件第1頁/共6頁第

2、一頁，共7頁。8.7.2 經(jīng)典(jngdin)Runge-Kutta 法的 Matlab 函數(shù)文件function x,y =rk4(f,tspan,y0,n) % 解初值問題：y =f(x,y),y(a)=y0。 % 使用(shyng) n 步 4 階 R-K 法。 a=tspan(1);b=tspan(2);h=(b-a)/n; k1=h*feval (f,a,y0); k2=h*feval (f,a+h/2,y0+k1/2); k3=h*feval (f,a+h/2,y0+k2/2); k4=h*feval (f,a+h,y0+k3); y(1)=y0+k1/6+k2/3+k3/3+k4

3、/6; 第2頁/共6頁第二頁，共7頁。for i=1:n-1 k1=h*feval(f,x(i),y(i); k2=h*feval(f,x(i)+h/2,y(i)+k1/2); k3=h*feval(f,x(i)+h/2,y(i)+k2/2); k4=h*feval(f,x(i)+h/2,y(i)+k3/2); y(i+1)=y(i)+k1/6+k2/3+k3/3+k4/6 ； end x=a x; y=y0 y;第3頁/共6頁第三頁，共7頁。function x,y =ad3 (f,tspan,y0,n) a=tspan(1);b=tspan(2); h=(b-a)/n;hh=h/12; x

4、=(a+h:h:b); %使用中點方法(fngf)提供開始值。 z0=feval (f,a,y0); k1=h*z0;k2=h*feval (f,a+h/2,y0+k1/2) ; y(1)=y0+k2; z(1)=feval (f,x(1)+h/2,y(1)+k1/2); 8.7.3 三階 Adams方法的 Matlab 函數(shù)(hnsh)文件第4頁/共6頁第四頁，共7頁。k1=h*z(1);k2=h*feval (f,x(1)+h/2,y(1)+k1/2) ; y(2)=y(1)+k2 ; % 繼續(xù)使用(shyng)三階 Adams 方法。 z(2)=feval (f,x(2),y(2); y(3)=y(2)+hh*(23*z(2)-16*z(1)+5*feval (f,a,y0) ; for i=3:n-1 z(i) =feval (f,x(2),y(2) ; y(i+1)=y(i)+hh(23*z(i)-16*z(i-1)+5z(i-2) : end x=a x ; y=y0 y; 第5頁/共6頁第五頁，共7頁。感謝您的觀看(gunkn)！第6頁/共6頁第六頁，共7頁。NoImage內(nèi)容(nirng)總結(jié)function x,y=euler(f,tspan,y0,n)。% 解初值問題：y =f(x,y),y(a)=y0。第3頁/共6頁。y(2)=y(1)+k2。