https://github.com/xieguigang/sciBASIC

根据积分表达式,微分方程的数值解关键在于微分方程的初值及计算微分方程式在tm(上一时刻)与tm+d(下一时刻)与坐标轴围成面积,若这个面积计算得越准确则得到的数值解也就越精确。微分表达式中与坐标轴围成的面积可表示如下,在实施算法的时候可以结合这个图更加直观点:

从上面的示意图可以看出,一段需要进行面积积分的曲线实际上是由多个梯形构成的多边形。那我们实际上只需要将这些梯形的面积都求出来,然后加起来就好了。

这里的梯形分割就是一种欧拉逼近的思想,欧拉逼近的几何意义,就是我们可以使用一段折线来近似的逼近一条曲线。

利用欧拉逼近,我们可以将一个精确的微分方程曲线

近似的使用线段来表示

但是用梯形面积的和来求解数值解,会存在一个精度问题:例如,如果你观察上图的梯形的顶端与曲线的交接的部分,会发现有空隙。这个空隙就导致了使用梯形面积求解的误差的产生。理论上我们只需要将横坐标轴无限分割,这样子梯形顶端与曲线的空隙就会无限接近于零。但是因为目前计算机精度和算力的限制,我们没有办法做到无限分割。

那有更好的方法来求解积分么?

龙格库塔法求解微分方程

虽然目前我们无法做到无限分割,求取得到最精确的解;但是好在我们可以通过应用改进的欧拉逼近4阶龙格库塔法来进行较为精确的常微分方程组的求解。

下面所示的代码为实现4阶龙格库塔法的伪代码

y(tm + 1) = y(tm) + (k1 + 2 k2 + 2 k3 + k4) / 6 * h   # (1)
k1 = f(y(tm), tm)                                     # (2)
k2 = f(y(tm) + k1 / 2 * h, tm + h / 2)                # (3)
k3 = f(y(tm) + k2 / 2 * h, tm + h / 2)                # (4)
k4 = f(y(tm) + k3     * h, tm + h)                    # (5)

在上面的计算公式中,y(tm + 1)表示下一时刻的数值解,y(tm)表示上一时刻数值解,f()表示微分方程表达式

请注意,y在这里是一个编程语言概念上的数组,所以下一个时刻的数值解是数组的下一个元素T+1,其等价于数学分析上的T+d

龙格库塔法的主要思想:在 Tm 点的附近选取一些特定的点,然后把这些点的函数值进行线性组合,使用组合值代替泰勒展开中 Tm 点的导数值。在上面的伪代码公式中:

  • k1 * h 是对阴影面积的一个逼近,欧拉似逼近。
  • k2 h 亦是对阴影面积的一个逼近,k1 / 2 h表示用欧拉法求得的是从当前时刻算起跨步长为h / 2的面积,则y(tm)  + k1 / 2 h表示在 tm + h/2处的数值近似解,则k2代表微分表达式tm + h / 2处的值,则2k2 h表示计算了阴影面积的两倍。
  • k3 h 更是对阴影面积的逼近,而且比k2 h更为准确,因为k2 h基于欧拉法求得原函数在tm + h / 2处的值从而得到微分表达式在tm + h / 2处的值,k3也是微分表达式在tm + h / 2处的值,只是基于k2,更加精确。2k3 h表示2倍的阴影面积。
  • k4 h是对阴影面积的一个逼近,采用y(tm) + h k3得到原函数tm+1时刻的数值解从而求得k4 = f( y(tm) + h k3, tm + h)  为微分表达式在tm + h时刻的解。则k4 h是表示用tm+1时刻的微分表达式值求得的阴影部分面积,也是用矩阵法求得的,属欧拉法。

https://blog.csdn.net/misskissc/article/details/8913941

二阶龙格-库塔法中,泰勒展开到了τ阶,通过与泰勒展开系数进行对比,可以得到含3个未知数的2个方程。依次类推,如果泰勒展开到了 τp-1 阶,可以得到 m 级 p 阶龙格-库塔法。

我们将上面的伪代码公式改写为具体的VisualBasic代码,代码如下:

''' <summary>
''' 四阶龙格库塔法解解微分方程,分块数量为n, 解的区间为[a,b], 解向量为(x,y),方程初值为(x0, y0)
''' 参考http://blog.sina.com.cn/s/blog_698c6a6f0100lp4x.html 和维基百科
''' </summary>
''' <param name="df"></param>
''' <param name="n"></param>
''' <param name="a"></param>
''' <param name="b"></param>
<Extension>
Public Sub RK4(ByRef df As ODE, n As Integer, a As Double, b As Double)
    Dim h As Double = (b - a) / n

    df.x = New Double(n - 1) {}
    df.y = New Double(n - 1) {}
    df.x(Scan0) = a
    df.y(Scan0) = df.y0

    Dim x = df.x, y = df.y
    Dim k1, k2, k3, k4 As Double

    For i As Integer = 1 To n - 1
        x(i) = a + h * i
        k1 = df(x(i - 1), y(i - 1))
        k2 = df(x(i - 1) + 0.5 * h, y(i - 1) + 0.5 * h * k1)
        k3 = df(x(i - 1) + 0.5 * h, y(i - 1) + 0.5 * h * k2)
        k4 = df(x(i - 1) + h, y(i - 1) + h * k3)
        y(i) = y(i - 1) + h / 6 * (k1 + 2 * k2 + 2 * k3 + k4)
    Next
End Sub

利用龙格库塔法求解劳伦兹混沌系统

1963年,Lorenz发现了第一个混沌吸引子——Lorenz系统,从此揭开了混沌研究的序幕,该系统也称为Lorenz混沌系统。Lorenz系统是数值试验中最早发现的呈现混沌运动的耗散系统,其状态方程为:

Lorenz equations

按照上面的劳伦兹方程,我们可以用R#脚本定义如下所示的常微分方程组进行模拟分析:

require(grDevices3D);
require(charts);

const sigma = 10;
const rho   = 28;
const beta  = 8 / 3;

# define the Lorenz system, and then run runge kutta 4
# via deSolve function.
const lorenz = [
    x -> sigma * (y - x),
    y -> x * (rho - z) - y,
    z -> x * y - beta * z
]
|> deSolve(
    y0 = list(x = 1, y = 1, z = 1),
    a  = 0,
    b  = 10
)
;

# export result as table
lorenz
|> as.data.frame
|> write.csv(file = `${dirname(@script)}/lorenz.csv`)
;

# and then do simple data visualization
bitmap(file = `${dirname(@script)}/lorenz.png`) {
    const view = camera(viewAngle = [30, 30, 30]);

    lorenz
    |> plot(
        vector = list(x = "x", y = "y", z = "z" ),
        camera = view,
        color  = "brown"
    )
    ;
}

在上面的示例代码中,我们的混沌动力学系统使用三个R#语言之中的lambda函数来表示。在设定好了对应常数和系统的初始状态后,将所有的对象都送入到deSolve函数中执行。对上面的常微分方程组通过4阶的龙格库塔法进行数值求解,对结果进行简单的作图可视化之后:

Latest posts by xie guigang (see all)

Attachments

One response

发表评论

您的电子邮箱地址不会被公开。 必填项已用*标注