微分方程(1)---微分方程的数值解法

常微分方程

常微分方程是由函数在某点导数值和一个与x,y相关函数组成的方程

y′=f(x,y)y'=f(x,y)

式子 y′=x2+y2y'=x^2+y^2 就是一种典型的常微分方程，并且这个方程很难用分离变量的方法解出。matlab程序表达为

function f=f_xy(x,y)
f=x^2-y^2;
end

1.欧拉法

利用欧拉法解微分方程的思路是利用线性定理对下一步的函数值进行估计，欧拉给出的公式为

y′=f(x,y)y'=f(x,y)

yn+1=yn+h∗y′y_{n+1}=y_n+h*y'

xn+1=xn+hx_{n+1}=x_n+h

欧拉法的误差为一阶，故称为一阶方法

其代码为

x0=0;
y0=1;
x=[x0];
y=[y0];
a=.01;
for h=0:a:1
y0=y0+a*f_xy(x0,y0);
x0=x0+a;
x=[x x0];
y=[y y0];
end
plot(x,y,'r+')

由欧拉法所解得的图像为

二阶龙哥库塔法

龙哥库塔法的出现时为了减小欧拉法带来的误差，其思想在于利用估计斜率的组合来获得真实斜率值得估计。

Bn=(An+An′)2Bn=\frac{(An+An')}2

An′=f(xn+1,yˆn+1)An'=f(x_{n+1},\widehat y_{n+1})

yˆn+1=xn+y′∗h\widehat y_{n+1}=x_n+y'*h

An=y′An=y'

BnBn 为对下一步斜率的最佳估计。即 yn+1=xn+Bn∗hy_{n+1}=x_n+Bn*h

x0=0;
y0=1;
x=[x0];
y=[y0];
a=.01;
for h=0:a:1
y1=y0;
x1=x0+a;
y1=y0+a*f_xy(x0,y0);
Hypf_xy=f_xy(x1,y1);
y0=y0+a*(f_xy(x0,y0)+Hypf_xy)/2;
x0=x0+a;
x=[x x0];
y=[y y0];
end
plot(x,y,'g+')
hold on

结果如图



绿色的线是二阶龙哥库塔法的结果，显然要比红线更加贴近真实情况（对于凹曲线而言，微分方程数值解会小于真实值）