常微分方程数值解的若干Matlab函数文件课件.pptVIP

第八章常微分方程数值解法 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]; 8.7　常微分方程数值解的若干Matlab函数文件 8.7.1 Euler 方法的 Matlab 函数文件 8.7.2 经典Runge-Kutta 法的 Matlab 函数文件 function [x,y] =rk4(f,tspan,y0,n) % 解初值问题：y’ =f(x,y),y(a)=y0。 % 使用 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/6; 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]; function [x,y] =ad3 (f,tspan,y0,n) a=tspan(1);b=tspan(2); h=(b-a)/n;hh=h/12; x=(a+h:h:b); %使用中点方法提供开始值。 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 函数文件 k1=h*z(1);k2=h*feval (f,x(1)+h/2,y(1)+k1/2) ; y(2)=y(1)+k2 ; % 继续使用三阶 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];