前言

科学计算中常遇到一阶微分方程的初值问题:

$$ y' = f(x,y), x\in[x_0,b],\\ y(x_0) = y_0. $$

以下介绍几种常用方法

欧拉法

公式

经过推导可以得到公式
$y_{n+1} = y_n+hf(x_n,y_n)$

代码

matlab代码如下
eular.m

syms x;
x=input('x0='); %初始值x0
y=input('y(x0)='); %初始值y0
f(x,y)=input('dy/dx='); %输入y'
h=input('h='); %步长
t=input('x_final=')-x/h; %x的目标值
for i=1:t
    y=y+h*f(x,y);
    x=x+h;
    disp(['x' num2str(i) '=' num2str(roundn(x,-4))]); %保留4位小数
    disp(['y' num2str(i) '=' num2str(roundn(y,-4))]);
end

改进欧拉公式

公式

经过推导可以得到公式

$$ y_p = y_n+hf(x_n,y_n),\\ y_c = y_n+hf(x_{n+1},y_p), \\ y_{n+1} = \frac{1}{2}(y_p+y_c). $$

代码

matlab代码如下
improved_eular.m

syms x;
x=input('x0='); %初始值x0
y=input('y(x0)='); %初始值y0
f(x,y)=input('dy/dx='); %输入y'
h=input('h='); %步长
t=input('x_final=')-x/h; %x的目标值
fp=@(x,y) y+h*f(x,y);
fc=@(x,y) y+h*f(x+h,fp(x,y));
for i=1:t
    y=1/2*(fp(x,y)+fc(x,y));
    x=x+h;
    disp(['x' num2str(i) '=' num2str(roundn(x,-4))]); %保留4位小数
    disp(['y' num2str(i) '=' num2str(roundn(y,-4))]);
end

四阶龙格-库塔法(RK)

公式

$$ y_{n+1} = y_n+\frac{h}{6}(K_1+2K_2+2K_3+K_4),\\ K_1 = f(x_n,y_n),\\ K_2 = f(x_n+\frac{h}{2},y_n+\frac{h}{2}K_1),\\ K_3 = f(x_n+\frac{h}{2},y_n+\frac{h}{2}K_2),\\ K4 = f(x_n+h,y_n+hK_3). $$

代码

RK_4.m

syms x;
x=input('x0='); %初始值x0
y=input('y(x0)='); %初始值y0
f(x,y)=input('dy/dx='); %输入y'
h=input('h='); %步长
t=input('x_final=')-x/h; %x的目标值
K1=@(x,y) f(x,y);
K2=@(x,y) f(x+h/2,y+h/2*K1(x,y));
K3=@(x,y) f(x+h/2,y+h/2*K2(x,y));
K4=@(x,y) f(x+h,y+h*K3(x,y));
for i=1:t
    y=y+h/6*(K1(x,y)+2*K2(x,y)+2*K3(x,y)+K4(x,y));
    x=x+h;
    disp(['x' num2str(i) '=' num2str(roundn(x,-4))]);
    disp(['y*' num2str(i) '=' num2str(roundn(y,-4))]);
    ya=1/3*exp(-50*x)+x^2;
    disp(['y' num2str(i) '=' num2str(roundn(ya,-4))]);
    disp(['delta=' num2str(abs(ya-y))]);
end

项目地址

本文代码已上传至GitHub