我必须执行一个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
值时的方程可以独立求解。在这种情况下,我借助四阶龙格-库塔方法得到了稳定准确的解。然而,一旦我打开耦合,结果就变得非常不稳定,相对于时间网格步长和函数g_i
, h_i
的显式形式。
我想尝试使用隐式龙格-库塔方案是合理的,在这种情况下它可能是稳定的,但如果我这样做,我将不得不评估大小为4*N*c
的巨大矩阵的逆,其中c
取决于方法的顺序(例如高斯-勒让德方法的c = 3
)在每一步。当然,该矩阵将主要包含零并具有块三对角线形式,但它似乎仍然非常耗时。
我有两个问题:
是否存在稳定的显式方法,即使耦合函数
g_i
和h_i
(非常)大也能工作?如果隐式方法确实是一个很好的解决方案,那么块三对角矩阵的最快反演方法是什么?目前,我只是执行一个简单的高斯方法,以避免由于矩阵的特定结构而产生的冗余操作。
可能对我们有帮助的其他信息和细节:
我目前考虑
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
。j
绝对值较大的u_j
和v_j
函数对于所有t
来说都很小,因此初始峰值没有达到"边界"。
To 1)如果您的函数非常大,将没有稳定的显式方法。这是由于显式(龙格-库塔)方法的稳定区域是紧凑的。
To 2)如果你的矩阵大于100x100,你可以使用这个方法:块三对角矩阵的逆和舍入误差。