现实当中的不少物理问题、工程问题都涉及到微分方程,其中微分方程有常微分方程(Ordinary Differential Equation)和偏微分方程(Partial Differential Equation)之分。一般来说,所谓的常微分方程是指只有一个自变量的方程,如 u' = 2*u+x+2 。
不少工程问题中涉及的微分方程,我们很难求出方程的解析解,或者说根本不存在精确的解析解。此时,我们需要利用电脑,结合数值分析的方法来近似求出微分方程的相关解,并研究其性质。通过求出多个自变量的值,并求出对应的解,那么可以绘制出图形来辅助研究方程的特征。
1 Runge–Kutta算法
根据百度百科的相关介绍,龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法。由于此算法精度高,采取措施对误差进行抑制,所以其实现原理也较复杂。该算法是构建在数学分析的基础之上的。一般的Runge-Kutta算法的形式如下:
这个公式看起来也比较的抽象,在实际计算中,一般采用4阶Runge-Kutta算法进行求值,其计算公式如下:
为了更加直观,这里选择一个示例,假设有如下的一个微分方程需要进行求解:
此时,我们需要计算一下u(0.2)对应的方程解是什么?那么此示例的解题过程如下:
实际上,本微分方程的解析解表达式为:
即解析式的解近似为 1.3472599 ,而4阶Runge-Kutta算法进行求值的结果为1.3472。由此可知,精确度还是可以的。上述示例来自网址 :
http://www.public.asu.edu/~hhuang38/example_Runge-Kutta.pdf
2 F# Runge–Kutta算法实现
下面给出F#的算法实现,示例代码如下:
letodeIntf (a:float) (b:float) x0f0=letn=5leth= (b-a)/float(n) letx= [|foriin0 .. n->a+float(i) *h|] letu= [|foriin0 .. (n+1) ->0.|] u.[0] <-f0fori=0tondoletk1=f(x.[i], u.[i]) letk2=f(x.[i]+h/2. , u.[i]+k1*h/2.) letk3=f(x.[i]+h/2. , u.[i]+k2*h/2.) letk4=f(x.[i]+h , u.[i]+k3*h) u.[i+1] <-u.[i] + (k1+2.*(k2+k3) +k4) *h/6. (x,u)
3 测试
下面给出测试示例,代码如下:
lettestODEInt=letf1(x,u) =-2.*u+x+4.letdx,du=odeIntf10.1.00.1.printfn"%A"dxprintfn"%A"du
执行示例代码,结果如下:
[|0.0; 0.2; 0.4; 0.6; 0.8; 1.0|] [|1.0; 1.3472; 1.61292288; 1.824023499; 1.998505354; 2.148437989; 2.281912828|]
由此可知,在0.2时,u(x)值为1.3472,与示例一致。