一阶常微分方程的初值问题求解.doc

一阶常微分方程的初值问题

1.欧拉方法
欧拉公式：
算法（欧拉方法）
（1）给定整数N；
（2）；
（3）
（4）对k=2,3,…N+1做：

（5）返回
基于欧拉方法的Matlab程序如下：
function[x,y]=odeEuler(f,y0,a,b,h)
n=(b-a)/h;
x=a:h:b;
y(1)=y0;
fori=1:n
y(i+1)=y(i)+h*feval(f,x(i),y(i));
end

2.改进的欧拉公式：

算法（改进欧拉方法）：
给定整数N；
；

对k=2,3,…N+1做：

返回
基于改进欧拉方法的Matlab程序如下：
function[x,y]=odeIEuler(f,a,b,y0,h)
n=(b-a)/h;
x=a:h:b;
y(1)=y0;
fori=1:n
k1=h*feval(f,x(i),y(i));
k2=h*feval(f,x(i+1),y(i)+k1);
y(i+1)=y(i)+0.5*(k1+k2);
end

3.一类应用广泛的高精度的显式单步法-----龙格-库塔（Runge-Kutta）方法，简称R-K方法
标准四阶四段龙格-库塔公式

算法（标准四阶四段龙格-库塔方法）
（1）给定整数N；
（2）；
（3）
（4）对k=2,3,…N+1做：

（5）返回
基于标准四阶四段龙格-库塔方法Matlab程序如下：
function[x,y]=oderk4(f,a,b,h,y0)
n=(b-a)/h;
x=a:h:b;
y(1)=y0;
fori=1:n;
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,y(i)+k3);
y(i+1)=y(i)+1/6*(k1+2*k2+2*k3+k4);
end
其中：f为常微分方程的右端项，a,求解区间的左右端点，h为自变量的步长，y0微分方程的初值，x,y分别为计算完成时的自变量取值和对应点上的函数值。

4.标准四阶四段龙格-库塔方法应用实例

functionz=f(x,y)
z=y-2x/y;

function[x,y]=oderk4(f,a,b,h,y0)
n=(b-a)/h;
x=a:h:b;
y(1)=y0;
fori=1:n;
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,y(i)+k3);
y(i+1)=y(i)+1/6*(k1+2*k2+2*k3+k4);
end

运行结果：
>>[x,y]=oderk4(‘f’,0,1,0.1,1)
x=

Columns1through9

00.10000.20000.30000.40000.50000.60000.70000.8000

Columns10through11

0.90001.0000

y=

Columns1through9

1.00001.09541.18321.26491.34161.41421.48321.54921.6125

Columns10through11

1.67331.7321