一阶常微分方程的初值问题求解
一阶常微分方程的初值问题
⎨
⎧y '=f (x , y ), ⎩y (a ) =y 0, x ∈[a , b ],
1. 欧拉方法
欧拉公式:y n +1=y n +hf (x n , y n )
算法(欧拉方法)
(1)给定整数N ;
(2)h =b -a
N ;
(3)y (1) =y 0;
(4)对k=2,3,…N+1做:
x =a +(k -1) h ;
y (k ) =y (k -1) +hf (x , y (k -1));
(5)返回y .
基于欧拉方法的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. 改进的欧拉公式:
⎧y n +1=y n +hf (x n , y n ), ⎪⎨h y =y +(f (x n , y n ) +f (x n +1, y n +1)), ⎪n +1n 2⎩
算法(改进欧拉方法):
(1) 给定整数N ;
(2)
(3) h =b -a N ; y (1) =y 0;
(4) 对k=2,3,…N+1做:
x =a +(k -1) h ;
y p =y (k -1) +hf (x , y (k -1));
y c =y (k -1) +hf (x +h , y p );
y (k ) =0. 5(y c +y p );
(5) 返回y .
基于改进欧拉方法的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⎧y =y +(k 1+2k 2+2k 3+k 4) n ⎪n +16⎪
⎪k 1=hf (x n , y n )
⎪11⎪ k =hf (x +h , y +k 1) ⎨2n n 22⎪
11⎪k =hf (x +h , y +k 2) n ⎪322⎪⎪⎩k 4=hf (x n +h , y n +k 3)
算法(标准四阶四段龙格-库塔方法)
(1)给定整数N ;
(2)
(3)h =b -a N ; y (1) =y 0;
(4)对k=2,3,…N+1做:
x =a +(k -1) h ;
k 1=hf (x , y (k -1));
k 2=hf (x +0. 5h , y (k -1) +0. 5k 1);
k 3=hf (x +0. 5h , y (k -1) +0. 5k 2);
k 4=hf (x +h , y (k -1) +k 3);
y (k ) =y (k -1) +(k 1+2k 2+2k 3+k 4) /6;
(5)返回y .
基于标准四阶四段龙格-库塔方法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. 标准四阶四段龙格-库塔方法应用实例
⎧y '=f (x , y ), x ∈[a , b ], ⎨y (a ) =y , 0⎩
已知f (x , y ) =y -2x
y , x ∈[0, 1],y (0) =1, 并取h =0. 1
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.4832 1.5492 1.6125
Columns 10 through 11
1.6733 1.7321 1.1832 1.2649 1.3416 1.4142