文章阅读目录大纲
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时刻的微分表达式值求得的阴影部分面积，也是用矩阵法求得的，属欧拉法。
二阶龙格-库塔法中，泰勒展开到了τ阶，通过与泰勒展开系数进行对比，可以得到含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系统是数值试验中最早发现的呈现混沌运动的耗散系统，其状态方程为：
按照上面的劳伦兹方程，我们可以用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阶的龙格库塔法进行数值求解，对结果进行简单的作图可视化之后：
- 【MZKit】简单自动化组织分区 - 2023年11月5日
- 【MZKit教程】质谱成像原始数据文件查看 - 2023年6月29日
- 生物序列图嵌入算法 - 2023年6月29日
One response
[…] 关于使用R#脚本进行动力学系统的模拟计算的原理，可以阅读之前的一篇文章：《R#语言使用龙格库塔法求解劳伦兹混沌系统》 […]