ODE系统的数值稳定性

Question

我必须对一个ODE系统进行数值求解，其形式如下：

du_j/dt = f_1(u_j, v_j, t) + g_1(t)v_(j-1) + h_1(t)v_(j+1),

dv_j/dt = f_2(u_j, v_j, t) + g_2(t)u_(j-1) + h_2(t)u_(j+1),

其中u_j(t)和v_j(t)是时间的复值标量函数，t、f_i和g_i是给定的函数，j = -N,..N是给定的。这是一个初值问题，任务是在一定时间T找到解决方案。

如果是g_i(t) = h_i(t) = 0，则可以独立求解不同j值的方程。在这种情况下，借助于四阶Runge方法，我得到了一个稳定而精确的解。然而，一旦我打开联轴器，结果就变得非常不稳定的时间网格步骤和明确形式的函数g_i，h_i。

我想尝试采用隐式Runge格式是合理的，在这种情况下可能是稳定的，但如果这样做，我将不得不计算一个巨大的4*N*c矩阵的逆，其中c取决于方法的顺序(例如，高斯-Legendre方法的c = 3 )。当然，矩阵将主要包含零，并有一个块三对角形式，但它似乎仍然很费时。

所以我有两个问题：

1. 是否有一种稳定的显式方法，即使耦合函数g_i和h_i (很大)也能工作？
2. 如果隐式方法确实是一个很好的解，那么块三对角矩阵反演最快的方法是什么？目前，我只是执行一个简单的高斯方法，避免了由于矩阵的特定结构而产生的冗余操作。

可能对我们有帮助的其他信息和细节：

• 我用Fortran 95。
• 我现在考虑的是g_1(t) = h_1(t) = g_2(t) = h_2(t) = -iAF(t)sin(omega*t)，其中i是虚单位，A和omega是给定的常数，F(t)是一个平滑的包络，首先从0到1，然后从1到0，所以F(0) = F(T) = 0。
• 最初是u_j = v_j = 0，除非是j = 0。u_j和v_j的绝对值很大的函数j对于所有的t都是非常小的，因此最初的峰值不达到“边界”。

Answer 1

如果你的函数很大，就不会有稳定的显式方法。这是因为显式(Runge)方法的稳定性范围很小。

如果矩阵大于100x100，则可以使用以下方法：块三对角矩阵的逆与舍入误差。