西京学院数学软件实验任务书
实验二十四实验报告
一、实验名称：欧拉数值算法（显式，隐式，欧拉预估-校正法），Runge-Kutta数值算法。
二、实验目的：进一步熟悉欧拉数值算法（显式，隐式，欧拉预估-校正法），Runge-Kutta数值算法。
三、实验要求：运用Matlab/C/C++/Java/Maple/Mathematica等其中一种语言完成程序设计。
四、实验原理：
1．欧拉数值算法（显式）：
微分方程里最简单的方程形式莫过于一阶常微分方程的初值问题，即：
其中
，
为常数。
因为其简单但又是求解其他方程的基础，所以发展了许多典型的解法。所有算法中的
就是代表上式中
，而
表示
，
表示
。
简单欧拉法是一种单步递推算法。简单欧拉法的公式如下所示：
简单欧拉法的算法过程介绍如下:
给出自变量
的定义域
，初始值
及步长
。
对
，计算
2．欧拉数值算法（隐式）：
隐式欧拉法也叫退欧拉法，隐式欧拉法的公式如下所示：
隐式欧拉法是一阶精度的方法，比它精度高的公式是：
隐式欧拉的算法过程介绍如下。
给出自变量的定义域，初始值及步长。
对，用牛顿法或其他方法求解方
得出
。
3．欧拉预估-校正法：
改进欧拉法是一种二阶显式求解法，其计算公式如下所示：
4．Runge-Kutta数值算法：
二阶龙格-库塔法有多种形式，除了改进的欧拉法外，还有中点法。中点法计算公式为：
五、实验内容：
%欧拉数值算法（显式）
function [h,k,X,Y,P,REn]=Qeuler1(funfcn,x0,y0,b,n,tol)
x=x0; h=(b-x)/n; X=zeros(n,1); y=y0;
Y=zeros(n,1); k=1; X(k)=x; Y(k)=y';
for k=2:n+1
fxy=feval(funfcn,x,y);
delta=norm(h*fxy,'inf');
wucha=tol*max(norm(y,'inf'),1.0);
if delta>=wucha
x=x+h; y=y+h*fxy; X(k)=x;Y(k)=y';
end
plot(X,Y,'rp')
grid
label('自变量 X'), ylabel('因变量 Y')
title('用向前欧拉（Euler）公式计算dy/dx=f(x,y)，y(x0)=y0在[x0,b]上的数值解')
end
P=[X,Y];
syms dy2,
REn=0.5*dy2*h^2;
%欧拉数值算法（隐式）
function [X,Y,n,P]=Heuler1(funfcn,x0,b,y0,h,tol)
n=fix((b-x0)/h); X=zeros(n+1,1); Y=zeros(n+1,1);
k=1; X(k)=x0; Y(k,:)=y0; Y1(k,:)=y0;
clc,x0,h,y0
for i=2:n+1
X(i)=x0+h; Y(i,:)=y0+h*feval(funfcn,x0,y0);
Y1(i,:)=y0+h*feval(funfcn,X(i),Y(i,:));
Wu=abs(Y1(i,:)-Y(i,:));
while Wu>tol
p=Y1(i,:);
Y1(i,:)=y0+h*feval(funfcn,X(i),p);
Y(i,:)=p;
end
x0=x0+h; y0=Y1(i,:);
Y(i,:)=y0; plot(X,Y,'ro')
grid on
xlabel('自变量 X'), ylabel('因变量 Y')
title('用向后欧拉公式计算dy/dx=f(x,y)，y(x0)=y0在[x0,b]上的数值解')
end
X=X(1:n+1); Y=Y(1:n+1,:); n=1:n+1; P=[n',X,Y]
%欧拉预估-校正法
function [H,X,Y,k,h,P]=gaiEuler(funfcn,x0,b,y0,tol)
pow=1/3;
if nargin < 5 | isempty(tol),
tol = 1.e-6;
end;
if nargin < 6 | isempty(trace),
trace = 0;
end;
x=x0; h=0.0078125*(b-x);
y=y0(:);p=128; n=fix((b-x0)/h);
H=zeros(p,1); X=zeros(p,1);
Y=zeros(p,length(y)); k=1;
X(k)=x; Y(k,:)=y';
while (x<b)&(x+h>x)
if x+h>b
h=b-x;
end
fxy=feval(funfcn,x,y); fxy=fxy(:);
delta=norm(h*fxy,'inf'); wucha=tol*max(norm(y,'inf'),1.0);
if delta<=wucha
x=x+h; y1=y+h*fxy; fxy1=feval(funfcn,x,y1);
fxy=fxy(:);y2=(fxy+fxy1)/2; y=y+h*y2; k=k+1;
if k>length(X)
X=[X;zeros(p,1)]; Y=[Y;zeros(p,length(y))];
H=[H;zeros(p,1)];
end
H(k)=h;X(k)=x;Y(k,:)=y'; plot(X,Y,'mh'), grid
xlabel('自变量 X'), ylabel('因变量 Y')
title('用改进的欧拉公式计算dy/dx=f(x,y)，y(x0)=y0在[x0,b]上的数值解')
end
if delta~=0.0
h=min(h*8,0.9*h*(wucha/delta)^pow);
end
end
if (x<b)
disp('Singularity likely.'), x
end
H=H(1:k); X=X(1:k); Y=Y(1:k,:); n=1:k; P=[n',H,X,Y]
%Runge-Kutta数值算法
function [k,X,Y,fxy,wch,wucha,P]=RK3(funfcn,fun,x0,b,C,y0,h)
x=x0; y=y0;p=128;
n=fix((b-x0)/h);
fxy=zeros(p,1);
wucha=zeros(p,1);
wch=zeros(p,1);
X=zeros(p,1);
Y=zeros(p,length(y));
k=1; X(k)=x; Y(k,:)=y';
clc,x,h,y
for k=2:n+1
x=x+h;a2=C(4);
a3=C(5);b21=C(6);
b31=C(7);
b32=C(8);c1=C(1);
c2=C(2);
c3=C(3);
x1=x+a2*h;
x2=x+a3*h;
k1=feval(funfcn,x,y);
y1=y+b21*h*k1;
k2=feval(funfcn,x1,y1);
y2=y+b31*h*k1+b32*h*k2;
k3=feval(funfcn,x2,y2);
fxy(k)=feval(fun,x);
y=y+h*(c1*k1+c2*k2+c3*k3);
X(k)=x; Y(k,:)=y;
k=k+1;
plot(X,Y,'rp',X,fxy,'bo'),
grid,xlabel('自变量 X'), ylabel('因变量 Y')
legend('用三阶龙格-库塔方法计算dy/dx=f(x,y)，y(x0)=y0在[x0,b]上的数值解','y/dx=f(x,y)，y(x0)=y0的精确解y=f(x)')
end
for k=2:n+1
wucha(k)=norm(Y(k-1)-Y(k));
wch(k)=norm(fxy(k)-Y(k));
end
X=X(1:k);
Y=Y(1:k,:);
fxy=fxy(1:k,:);
n=1:k;
wucha=wucha(1:k,:);
wch=wch(1:k,:);
P=[n',X,Y,fxy,wch,wucha];
第二篇：数学实验“微分方程组数值算法——四阶Runge-Kutta数值算法”实验报告(内含matlab程序)
西京学院数学软件实验任务书
实验二十六实验报告
一、实验名称：微分方程组数值算法——四阶Runge-Kutta数值算法。
二、实验目的：进一步熟悉微分方程组数值算法——四阶Runge-Kutta数值算法。
三、实验要求：运用Matlab/C/C++/Java/Maple/Mathematica等其中一种语言完成程序设计。
四、实验原理：
四阶Runge-Kutta数值算法：
对于求解一阶微分方程组问题
由初值问题的经典Runge-kutta公式可得一阶常微分方程组初值问题的Runge-kutta公式：
五、实验内容：
%四阶Runge-Kutta数值算法
function y=DELGKT4_lungkuta(f, h,a,b,y0,varvec)
format long;
N=(b-a)/h;
y=zeros(N+1,1);
y(1)=y0;
x=a:h:b;
var=findsym(f);
for i=2:N+1
K1=Funval(f,varvec,[x(i-1) y(i-1)]);
K2=Funval(f,varvec,[x(i-1)+h/2 y(i-1)+K1*h/2]);
K3=Funval(f,varvec,[x(i-1)+h/2 y(i-2)+K2*h/2]);
K4=Funval(f,varvec,[x(i-1)+h y(i-1)+h*K3]);
y(i)=y(i-1)+h*(K1+2*K2+2*K3+K4)/6;
end
format short;