王能超 计算方法——算法设计及MATLAB实现课后代码

算例6：对于微分方程解：

function z=f2(x,y) z=x^2-y;

为衡量数值解的精度，我们求出该方程的解析解为y??e?x?x2?2x?2.在此也以文件的形式表示如下： function y=solvef2(x) y=-exp(-x)+x^2-2*x+2;

令f=@f2; a=0; b=1; N=10; ya=1; 运行：E=MendEuler(f,a,b,N,ya)； y=solvef2(a:(b-a)/N:b); m=[E,y’]

其中m内的y’表示准确解.

3.2 四阶Runge-Kutta方法

用四阶Runge-Kutta方法求解常微分方程. MATLAB文件：（文件名：RungKutta4.m） function R=RungKutta4(f,a,b,N,ya) %f是微分方程右端函数句柄

%a,b是自变量的取值区间[a,b]的端点 %N是区间等分的个数 %ya表初值y(a)

%R=[x’,y’]是自变量X和解Y所组成的矩阵 h=(b-a)/N;

x=zeros(1,N+1); y=zeros(1,N+1); x=a:h:b; y(1)=ya; for i=1:N

k1=feval(f,x(i),y(i));

k2=feval(f,x(i)+h/2,y(i)+(h/2)*k1); k3=feval(f,x(i)+h/2,y(i)+(h/2)*k2); k4=feval(f,x(i)+h,y(i)+h*k3);

y(i+1)=y(i)+(h/6)*(k1+2*k2+2*k3+k4); end

R=[x’,y’];

算例7：对于微分方程解：

9

dy?x2?y,y(0)?1,0?x?1，用改进的Euler方法求解. dxdy?x2?y,y(0)?1,0?x?1，用四阶Runge-Kutta方法求解. dxfunction z=f2(x,y) z=x^2-y;

为衡量数值解的精度，我们求出该方程的解析解为y??e?x?x2?2x?2.在此也以文件的形式表示如下： function y=solvef2(x) y=-exp(-x)+x^2-2*x+2;

令f=@f2; a=0; b=1; N=10; ya=1; 运行：R=RungKutta4(f,a,b,N,ya)； y=solvef2(a:(b-a)/N:b); m=[R,y’]

其中m内的y’表示准确解,精度比前者高.

3.3 二阶Adams预报校正系统

用二阶Adams预报校正系统求解常微分方程. MATLAB文件：（文件名：Adams2PC.m） function A=Adams2PC(f,a,b,N,ya) %f是微分方程右端函数句柄

%a,b是自变量的取值区间[a,b]的端点 %N是区间等分的个数 %ya表初值y(a)

%A=[x’,y’]是自变量X和解Y所组成的矩阵 h=(b-a)/N;

x=zeros(1,N+1); y=zeros(1,N+1); x=a:h:b; y(1)=ya; for i=1:N if i==1

y1=y(i)+h*feval(f,x(i),y(i)); y2=y(i)+h*feval(f,x(i+1),y1); y(i+1)=(y1+y2)/2; dy1=feval(f,x(i),y(i));

dy2=feval(f,x(i+1),y(i+1)); else

y(i+1)=y(i)+h*(3*dy2-dy1)/2; P=feval(f,x(i+1),y(i+1));l y(i+1)=y(i)+h*(P+dy2)/2; dy1=dy2;

dy2=feval(f,x(i+1),y(i+1)); end end A[x’,y’];

10

3.4 改进的四阶Adams预报校正系统

用改进的四阶Adams预报校正系统求解常微分方程. MATLAB文件：（文件名：CAdams4PC.m） function A=CAdams4PC(f,a,b,N,ya) %f是微分方程右端函数句柄

%a,b是自变量的取值区间[a,b]的端点 %N是区间等分的个数 %ya表初值y(a)

%A=[x’,y’]是自变量X和解Y所组成的矩阵 if N<4 break; end

h=(b-a)/N;

x=zeros(1,N+1); y=zeros(1,N+1); x=a:h:b; y(1)=ya; F=zeros(1,4); for i=1:N

if i<4 %用四阶Runge-Kutta方法求初始解 k1=feval(f,x(i),y(i));

k2=feval(f,x(i)+h/2,y(i)+(h/2)*k1); k3=feval(f,x(i)+h/2,y(i)+(h/2)*k2); k4=feval(f,x(i)+h,y(i)+h*k3);

y(i+1)=y(i)+(h/6)*(k1+2*k2+2*k3+k4); elseif i==4

F=feval(f,x(i-3:i),y(i-3:i));

py=y(i)+(h/24)*(F*[-9,37,-59,55]’); %预报 p=feval(f,x(i+1),py); F=[F(2) F(3) F(4) p];

y(i+1)=y(i)+(h/24)*(F*[1,-5,19,9]’); %校正 p=py; c=y(i+1); else

F=feval(f,x(i-3:i),y(i-3:i));

py=y(i)+(h/24)*(F*[-9,37,-59,55]’); %预报 my=py-251*(p-c)/270; %改进 m=feval(f,x(i+1),my); F=[F(2) F(3) F(4) m];

cy=y(i)+(h/24)*(F*[1,-5,19,9]’); %校正 y(i+1)=cy+19*(py-cy)/270; %改进 p=py; c=cy; end

11

end A=[x’,y’];

附件(补充)：四阶Adams预报校正系统的程序： MATLAB文件：（文件名：Adams4PC.m） function A=Adams4PC(f,a,b,N,ya) %f是微分方程右端函数句柄

%a,b是自变量的取值区间[a,b]的端点 %N是区间等分的个数 %ya表初值y(a)

%A=[x’,y’]是自变量X和解Y所组成的矩阵 if N<4 break; end

h=(b-a)/N;

x=zeros(1,N+1); y=zeros(1,N+1); x=a:h:b; y(1)=ya; F=zeros(1,4); for i=1:N

if i<4 %用四阶Runge-Kutta方法求初始解 k1=feval(f,x(i),y(i));

k2=feval(f,x(i)+h/2,y(i)+(h/2)*k1); k3=feval(f,x(i)+h/2,y(i)+(h/2)*k2); k4=feval(f,x(i)+h,y(i)+h*k3);

y(i+1)=y(i)+(h/6)*(k1+2*k2+2*k3+k4); else

F=feval(f,x(i-3:i),y(i-3:i));

py=y(i)+(h/24)*(F*[-9,37,-59,55]’); %预报 p=feval(f,x(i+1),py); F=[F(2) F(3) F(4) p];

y(i+1)=y(i)+(h/24)*(F*[1,-5,19,9]’); %校正 end end A=[x’,y’];

算例8：分别用二阶Adams预报校正系统，四阶Adams预报校正系统和改进四阶Adams预报校正系统求解如下微分方程初值问题：解：

function z=f3(x,y) z=-y+x+1;

12

dy??y?x?1,y(0)?1,0?x?1. dx