欧拉法求解微分方程

Posted conmajia

tags:

篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了欧拉法求解微分方程相关的知识,希望对你有一定的参考价值。

by Conmajia
2014

原文是我在2014年写的,用C# 完成,这里改成javascript了,特基础。当然最方便的还是用数学库,或者Matlab、Mathematics这些数学软件(如果你只求值的话),或者可以换成C、Java、Go、Erlang任何其他的语言实现。欧拉(Euler)和中心差分逼近,是最朴素的想法,可惜代数精度太低了,而龙格库塔的稳定性又是个问题。总之只能用来计算普通的东西,高大上问题一般还是使用广义$alpha$,Wilson-$Theta$之类的,利用系数调节人工阻尼,让低频的部分少受算法影响,高频阻尼增大,增加微分方程的稳定性,又不影响算法的精度。
Math.NET计划是用于.NET的数学运算库,包括了数值计算库Iridium.NET。这是非常古老的的库,最近一次更新已经是10年前了。但是,数学的东西永远不会过时,所以仔细研究一下,还是可以吃到很多干货的。

$ git clone git://github.com/mathnet/mathnet-iridium.git

工程中有时候需要解微分方程,比如这种:

[ y'=y-cfrac{2x}{y} ]

(y(0)=1). 两边积分,得到:

[ egin{align*} int_{x_n}^{x_{n+1}}y'mathrm{d} x&=int_{x_n}^{x_{n+1}}left(y-cfrac{2x}{y} ight)mathrm{d} x y_{n+1}-y_n&=int_{x_n}^{x_{n+1}}left(y-cfrac{2x}{y} ight)mathrm{d}x y&=sqrt{1+2x} end{align*} ]

工程上只想要数值解,一般采用差分近似代替微分. 最简单最朴素的办法是用欧拉公式

[ ag{1} y_{n+1}=y_n+Delta x f(x_n,y_n)quad n=0,1,2,cdots ]

推导很简单,对(x_n),有:

[ y'_n=f(x_n,y_n). ]

(y'(x_n))有:

[ ag{2} y'_napprox cfrac{y_{n+1}-y_n}{Delta x}. ]

式(2)左边叫微商(不是喜提迪丽热巴的那个微商),右侧叫差商(不是淘宝伪劣产品卖家),(Delta x=left|x_{n+1}-x_n ight|).

代回式(2),有:

[ egin{align*} y'_n=f(x_n,y_n) Rightarrowcfrac{y_{n+1}-y_n}{Delta x} &approx f(x_n,y_n) y_{n+1}-y_n &approx Delta x f(x_n,y_n) y_{n+1} &approx y(x_n)+Delta x f(x_n,y_n). end{align*} ]

式(1)的欧拉公式成了:

[ ag{3} y_{n+1}=y_n+Delta xleft(y_n-cfrac{2x_n}{y_n} ight). ]

假设用myFn函数表示余项(y_n-dfrac{2x_n}{y_n})

function myFn(x, y) {
  return y - 2 * x / y;
}

Euler可以这样实现:

function calculate(x0, y0, delta, xn) {
  var yn;
  while(x0 < xn) {
    yn = y0 + delta * myFn(x0, y0);
    y0 = yn;
    x0 = x0 + delta;
  }
  return yn;
}

现在可以开始试验了。

理论上:

[ y=sqrt{1+2x}Rightarrow y(1)=sqrt{3}approx 1.7321 ]

(sqrt{3})是方程的真值,程序的目标是通过计算,得到尽量接近真值的结果.

(x_0=0)(y_0=1)(Delta x=0.1),程序计算结果为:

ans = 1.784771

误差3.04%. 减小(Delta x),比如(Delta x=0.0001),计算结果为:

ans = 1.732112

误差0.0007%十分接近真值了.

虽然这个方法可以求到比较精确的解,但是(Delta x)太小的话,while会执行很多次,效率低下.

引入定积分的梯形公式:

[ int_{x_n}^{x_{n+1}}f(x,y)mathrm{d}xapproxcfrac{Delta x}{2}left[f(x_n,y_n)+f(x_{n+1},y_{n+1}) ight] ]

式(3)可以写成:

[ ag{4} y_{n+1}approx y_n+cfrac{Delta x}{2}left[f(x_n,y_n)+f(x_{n+1},y_{n+1}) ight]. ]

式(3)和式(4)的特点是,前者速度快,精度低;后者速度慢,精度高. 所以对于粗算,可以使用:

[ y'_{n+1}=y_n+Delta x f(x_n,y_n) ]

精算,可以用:

[ y_{n+1}=y_n+cfrac{Delta x}{2}left[f(x_n,y_n)+f(x_{n+1},y'_{n+1}) ight]. ]

结合起来,先通过粗算得到(y'_{n+1})的近似值,再进行精算,得到高精度的最终解. 这样的话,既保证了计算结果的准确度,又没有消耗太多的计算资源,保证了计算效率. 下面是改进后的计算程序:

function calculate(x0, y0, delta, xn) {
  var yp, yc;
  while(x0 < xn) {
    yp = y0 + delta * myFn(x0, y0);
    yc = y0 + delta * myFn(x0 + delta, yp);
    y0 = 1 / 2 * (yp + yc);
    x0 = x0 + delta;
  }
  return y0;
}

(x_0=0)(y_0=1)(Delta x=0.1),计算结果为:

ans = 1.737867

误差0.33%.

和最早版本(误差3.04%)相比,在(Delta x)相同——意味着循环次数相近——的情况下,精度提高接近10倍.

具体问题具体分析,手工求微分方程基本是这么个思路. 当然,对于《数值分析》和《工程数学》这些课程来说,我上面写的东西不过是小儿科了. 更多的时候,做个伸手党其实很不错的,大把的现成玩意儿可以用,我干嘛还要费那劲呢?

The End. (Box)

以上是关于欧拉法求解微分方程的主要内容,如果未能解决你的问题,请参考以下文章

MATLAB-1个实例-欧拉法改进欧拉法ode45求解微分方程

欧拉法求解微分方程

数值分析:用改进欧拉法解微分方程初值问题(vf编程) 100

MATLAB从入门到精通-欧拉法与梯形法求解微分方程(含MATLAB源码)

matlab用龙格库塔法求解变系数常微分方程

matlab代码实现四阶龙格库塔求解微分方程