2016-08-16 72 views
2

我必须执行其具有如下形式的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_ig_i是给定函数和j = -N,..N。这是一个初始值问题,任务是在特定时间T找到解决方案。

如果g_i(t) = h_i(t) = 0,那么可以独立求解不同值j的方程。在这种情况下,借助四阶Runge-Kutta方法获得稳定和精确的解。但是,一旦打开联轴器,就时间网格步骤和函数的明确形式而言,结果变得非常不稳定,其形式为g_i,​​。

我认为尝试使用隐式Runge-Kutta方案是合理的,这种方法在这种情况下可能是稳定的,但如果我这样做,我将不得不评估一个大小为4*N*c的矩阵的逆,其中c取决于该方法的顺序(例如用于高斯 - 勒让德方法的c = 3)。当然,矩阵将大部分包含零并且具有块三对角形式,但它似乎仍然非常耗时。

所以我有两个问题:

  1. 有其工作的稳定明确的方法即使是在连接功能g_i和​​是(非常)大?

  2. 如果一个隐式方法确实是一个好的解决方案,块三对角矩阵反演的最快方法是什么?目前我只是执行简单的高斯方法,避免了由于矩阵的特定结构而产生的冗余操作。

附加信息和细节,可以帮助我们:

  • 我使用的Fortran 95

  • 我目前考虑g_1(t) = h_1(t) = g_2(t) = h_2(t) = -iAF(t)sin(omega*t),其中i是虚数单位,Aomega给出常数和F(t)是一个光滑的包络,它首先从0到1,然后是从1到0,因此F(0) = F(T) = 0

  • 最初u_j = v_j = 0除非j = 0。对于所有t而言,绝对值为j的函数u_jv_j对于所有的t来说都是非常小的,所以初始峰值没有达到“边界”。

回答