一阶常微分方程的初值问题求解.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;
for i=1:n
y(i+1)=y(i)+h*feval(f,x(i),y(i));
end

2.改进的欧拉公式：

算法（改进欧拉方法）：
(1) 给定整数N；
(2) ；
(3)
(4) 对k=2,3,…N+1做：

(5) 返回
基于改进欧拉方法的Matlab程序如下：
function [x,y]=odeIEuler(f,a,b,y0,h)
n=(b-a)/h;
x=a:h:b;
y(1)=y0;
for i=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;
for i=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.标准四阶四段龙格-库塔方法应用实例

function z= 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;
for i=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 =

Columns 1 through 9

0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000 0.8000

Columns 10 through 11

0.9000 1.0000

y =

Columns 1 through 9

1.0000 1.0954 1.1832 1.2649 1.3416 1.4142 1.4832 1.5492 1.6125

Columns 10 through 11

1.6733 1.7321