Численного устойчивость системы ОДЫ

Question

Я должен выполнить численное решение системы ОДЫ, которая имеет следующий вид:

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 для метода Гаусса-Лежандра) на каждом шаге. Конечно, матрица будет в основном содержать нули и иметь трехдиагональную форму блока, но она по-прежнему кажется очень трудоемкой.

Поэтому у меня есть два вопроса:

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

К 1) Не будет явного метода стабилизации, если ваши функции очень большие. Это связано с тем, что область устойчивости явных (Рунге-Кутта) методов компактна.

2) Если ваши матрицы больше, чем 100x100, вы можете использовать этот метод: Inverses of Block Tridiagonal Matrices and Rounding Errors.