「微分方程求解」常微分方程的数值求解

微分方程求解

常微分方程

首先理解一下什么是常微分方程,简单的说就是只有一个未知数的微分方程,具体定义如下:

凡含有参数，未知函数和未知函数导数 (或微分) 的方程，称为微分方程，有时简称为方程，未知函数是一元函数的微分方程称作常微分方程，未知函数是多元函数的微分方程称作偏微分方程。

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

{y′=f(x,y),y(x0)=y0," role="presentation"> ${\begin{cases} y^{'} = f (x, y), \\ y (x_{0}) = y_{0}, \end{cases}$

其中y=y(x)" role="presentation"> $y = y (x)$ 是未知函数,y(x0)=y0" role="presentation"> $y (x_{0}) = y_{0}$ 是初值条件,而f(x,y)" role="presentation"> $f (x, y)$ 是给定的二元函数

简单的数值方法和欧拉公式

简单的数值方法就是用差商代替导数,公式如下:

yn+1=yn+hf(xn,yn)(n=0,1,2,⋯)" role="presentation"> $y_{n + 1} = y_{n} + h f (x_{n}, y_{n}) (n = 0, 1, 2, \dots)$

其中h" role="presentation"> $h$ 是步长;

改进的欧拉方法:

{y¯n+1=yn+hf(xn,yn)预测yn+1=yn+h2[f(xn,yn)+f(xn+1,y¯n+1)]校正" role="presentation"> ${\begin{cases} {\bar{y}}_{n + 1} = y_{n} + h f (x_{n}, y_{n}) & 预测 \\ y_{n + 1} = y_{n} + \frac{h}{2} [f (x_{n}, y_{n}) + f (x_{n + 1}, {\bar{y}}_{n + 1})] & 校正 \end{cases}$

龙格-库塔(Runge-Kutta)方法

欧拉方法是龙格-库塔方法的一个特例,其局部截断误差为一阶泰勒余项O(h2)" role="presentation"> $O (h^{2})$ ,为了使误差更小,我们可以做更高阶的误差截断,这也就是我们Runge-Kutta方法的基本原理,具体推导可参考《数值分析》的第八章.其公式如下:

{yn+1=yn+h∑i=1rciKi,K1=f(xn,yn),Ki=f(xn+λih,yn+h∑j=1i−1uijKj)(i=2,3,⋯,r)," role="presentation"> ${\begin{cases} y_{n + 1} = y_{n} + h \sum_{i = 1}^{r} c_{i} K_{i}, \\ K_{1} = f (x_{n}, y_{n}), \\ K_{i} = f (x_{n} + λ_{i} h, y_{n} + h \sum_{j = 1}^{i - 1} u_{i} j K_{j}) (i = 2, 3, \dots, r), \end{cases}$

当r=1时就是欧拉方法,当r=2时,就是改进的欧拉方法,这里我们不做具体推导,而是看一下Matlab中封装好的4阶Runge-Kutta方法的函数实现ode45函数.

ode45

先看一个简单的例子:dydx=y+3xx2" role="presentation"> $\frac{d y}{d x} = \frac{y + 3 x}{x^{2}}$ ,初值y(0)=−2" role="presentation"> $y (0) = - 2$ ,求解区间为[14]" role="presentation"> $[14]$ ,代码如下:

odefun=@(x,y) (y+3*x)/x^2;
tspan=[1 4];
y0=-2;
[x y]=ode45(odefun,tspan,y0)
plot(x,y)

这里写图片描述

但是我们再看另外一个例子:

{y′=yx−2y20<x<3y(0)=0" role="presentation"> ${\begin{cases} y^{'} = \frac{y}{x} - 2 y^{2} 0<x<3 \\ y (0) = 0 \end{cases}$

我们编写如下代码:

odefun=@(x,y) y/x-2*y^2;
tspan=[0 3];
y0=0;
[x y]=ode45(odefun,tspan,y0)
plot(x,y)
disp(y)

发现没有数值解

这里写图片描述

于是我们把代码进行了如下修改:

odefun=@(x,y) fun(x,y);
tspan=[0 3];
y0=0;
[x y]=ode45(odefun,tspan,y0);
plot(x,y)

function dy = fun(x,y)
if(x==0)
    dy=1-2*y^2
else
    dy = y/x-2*y^2
end
end

这里写图片描述

这里需要强调一点的是,Runge-Kutta法针对的是连续的函数f" role="presentation"> $f$ ,由于在例子中,函数在x=0" role="presentation"> $x = 0$ 处是不连续的,所以在这个地方是需要单独处理的,在这里,yx" role="presentation"> $\frac{y}{x}$ 在x=0" role="presentation"> $x = 0$ 处的导数为1(洛必达法则,即00=1" role="presentation"> $\frac{0}{0} = 1$ )

常微分方程的数值求解