c - 使用 C 中的 Runge Kutta 求解二阶偏微分方程组

Question

我在使用龙格库塔算法求解微分方程组时遇到问题。到目前为止，我已将二阶偏微分方程重写为一组两个耦合方程，其中

    f(L1,L2) = L2
    g(L1,L2) = A*(B*L1-C*L2-D)

是两个方程，A、B、C 和 D 是常数。为了获得下一步的值，我对每个时间步 dt 进行如下操作:

    k1 = f(L1,L2)
    l1 = g(L1,L2)

    k2 = f(L1 + 0.5 * dt * k1,L2 + 0.5 * dt * l1 )
    l2 = g(L1 + 0.5 * dt * k1,L2 + 0.5 * dt * l1 )

    k3 = f(L1 + 0.5 * dt * k2,L2 + 0.5 * dt * l2 )
    l3 = g(L1 + 0.5 * dt * k2, L2 + 0.5 * dt * l2 )

    k4 = f(L1 + dt * k1,L2 +  dt * l1 )
    l4 = g(L1 + dt * k1,L2 + dt * l1 ) 

我使用当前时间步长的 L1 和 L2 值并迭代计算系数。

结果，我通过最后对系数进行求和和加权得到了 L1 和 L2。我的问题是，整个算法在 4 个时间步长后变得不稳定。

有谁知道这种实现在技术上是否正确？谢谢!

Answer 1

只是猜测，因为您没有说出您使用的 dt 值:使其尽可能小，因为

local truncation error is on the order of O(h^5), while the total accumulated error is order O(h^4).

(引自 this 维基百科文章，dt 扮演 h 角色)。