ЖУРНАЛ ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ И МАТЕМАТИЧЕСКОЙ ФИЗИКИ, 2007, том 47, < 7, с. 1187-1191
УДК 519.63
НОВЫЕ МЕТОДЫ РАСЩЕПЛЕНИЯ ДЛЯ ДВУМЕРНЫХ ЭВОЛЮЦИОННЫХ УРАВНЕНИЙ
© 2007 г. Н. В. Широбоков
(454080 Челябинск, пр-т Ленина, 76, Южно-Уральский гос. ун-т) e-mail: snv@susu.ac.ru Поступила в редакцию 15.12.2006 г.
Предлагаются новые методы расщепления второго и третьего порядков для дифференциальных уравнений в частных производных эволюционного типа в двумерном пространстве. Вывод методов основан на виде диагонально-неявных методов, применяемых для численного решения жестких обыкновенных дифференциальных уравнений. Найденные методы являются абсолютно и безусловно устойчивыми. Приведены тестовые расчеты. Библ. 4. Табл. 1.
Ключевые слова: двумерные эволюционные уравнения, методы Рунге-Кутты, диагонально-неявные методы, жесткие задачи, методы расщепления второго и третьего порядков.
1. ОСНОВНЫЕ ОПРЕДЕЛЕНИЯ Пусть в эволюционном уравнении
2
ut = Lu + f(t, x), u = u(t, x), x e U , (1)
дифференциальный оператор L является суммой операторов L1 и L2. Для численного решения уравнения (1) будем применять метод прямых, поэтому в дальнейшем будем считать все перечисленные операторы дискретными аппроксимациями по пространственным переменным. Оператор L1 зависит от производной по x, а L2 - по у. Ниже речь идет только о порядках аппроксимации по временной переменной t.
Предлагаются новые методы решения уравнения (1), основанные на идее построения многостадийных диагонально-неявных методов типа Рунге-Кутты, применяемых для решения жестких обыкновенных дифференциальных уравнений
s
Un + 1 = Un + т£bi (LYi + fi), (2)
i = 1
i - 1
Yi = Un + T £ aj (LYj + f) + Ta^Yi + f), i = 1, 2, ..., s, (3)
j = i
где
Un ~ u(tn, x), т = tn + 1 - tn, fi = f(tn + CiT, x), Лi = L - TaiLL.
Коэффициенты методов aij, ci, bi подбираются так, чтобы порядок аппроксимации метода равнялся s. Ниже находятся условия достижимости указанного порядка аппроксимации и коэффициенты предложенных методов.
Отметим, что уравнения (3) при любом i имеют расщепимый вид
i-i
(E - TaiL)(E - TaL)Yi = Un + т £ a4 (LYj + f) + Taf
j = i
который можно представить через две легко программируемые стадии:
г -1
(Е - тщ^Ж = ип + т ^ аг] (LYj + £) + та,/-, (Е - та^)^ =
1 = 1
Коэффициенты диагонально-неявного ^-стадийного метода Рунге-Кутты будем представлять таблицей Бутчера
|А c1 a11 • • a1,
ьт c, a,1 • • ass
b1 • ■ bs
В дальнейшем предполагается, что |b| < 1, 0 < ci < 1, i = 1, 2, ..., 5.
2. УСЛОВИЯ АППРОКСИМАЦИИ ВТОРОГО И ТРЕТЬЕГО ПОРЯДКОВ
Пусть 5 = 2, ип = и(п х). Порядок аппроксимации метода (2), (3) равен 2, если ип + 1 - и^п + 1, х) = = 0(т3). Применяя к данной разности формулу Тейлора, получаем условия для второго порядка аппроксимации:
bi + b2 = 1, bxan + b2a2i + b2a22 = 1 , biCi + b2C2 = 2-.
(4)
Коэффициенты Ь1, Ь2, а21 уравнений (4) выражаем через а11, а22, с1, с2. В результате получаем общий вид двухстадийных методов расщепления второго порядка аппроксимации:
a11 a22 C1 + С2 + 2 C1 a22 2c2al1
1-2 c1
0
22
2 C2-1
1-2 c1
2(с2-С1) 2(с2- С1)
При 5 = 3 условиями аппроксимации третьего порядка являются уравнения
Ь1 + Ь2 + Ь3 = 1,
Ь1С1 + Ь2С2 + ЬзСз = 1/2,
Ь1 е\ + Ь2 с2 + Ь3 с3 = 1/3,
Ь1ап + Ь2(а21 + а22) + Ьз(аз1 + аз2 + азз) = 1/2, Ь1апС1 + Ь2(а21С1 + а22С2) + Ьз(аз1С1 + аз2С2 + аззСз) = 1/6, Ь1 а2п + Ь2 а222 + Ь3 а233 = 0, Ь а21 + Ь2(а21(ап + а22) + а22) + Ьз(аз1(ап + азз) + аз2(а21 + а22 + азз) + а3з) = 1/6.
(5)
(6)
(7)
(8) (9)
(10) (11)
Заметим, что уравнения (5)-(9) совпадают с обычными условиями аппроксимации неявных методов Рунге-Кутты (см. [1]). Заметным отличием от стандартных условий аппроксимации методов Рунге-Кутты является уравнение (10), которое показывает, что не существует методов (2), (3) третьего порядка аппроксимации с положительными весами Ь. Исключаем последовательно коэффициенты а21, а31, а32, Ь1, Ь2, Ь3, с2 из уравнений (9), (11), (8), (5), (6), (7), (10). Следовательно, коэффициенты метода будут выражены через 5 параметров а11 = а, а22 = в, а33 = у, с1 = с, с3 = d:
c2 =
2,^ 2 2 2 . ^ 2 i2 Y +3у c - а + 3 a d
±4W
3 (Y2 - а2 + 2a2d - 2Y2c)
b1 =
2-3d + 6dc2 - 3c2 6(c - d)(c - c2)
c
a
c
2
a21 =
a32 =
2-3 d + 6 dc -3c 2-3 c + 6 cc2-3 c2
b = -, bi = -,
2 6(С - C2)(d - C2) ' 3 6(C2- d)(С - d)
[(6 c2 - 3 )Y - 3c2 + 1 ]a + [(3-6 c )Y + 3 c - 1 ]ß + 3 (c - c 2 )y - c + c2
1-3 c -3 y + 6 c y '
- c){[(6d- 3)c - 3d + 2]ß + [(3 - 6c2)c + 3c2 - 2]y + 3(c2 - d)c + d- c2}
a31
(2-3c - 3c2 + 6cc2)(c2 - c) _V__
[ a32 b 3 + (ß - Y) b 2 ] b 3 c ,
где
Ж = (1-6 d + 15 d2- 18 d3 + 9 et )a4 + (1-6 c +15 c2- 18 c3 + 9c4 )y 4 + + 3 [ 6 (1-2 d )c3- (7-6 d - 12 d2) c2-2 (1+2 d -6 d2) c -2 d + 3 d2 ] 2 + + { 3 [ 3 (1-4 d + 4 d2) c2 - 2 (1-2 d -3 d2 + 6d3) c + 2 d -7 d 2 + 6 d3 ]ß2 + + [ 3 (- 1 + 6 d -6 d2) c2 + 6 (1-4 d + 4 d2) c -2 + 6 d -3 d2 ]2 }a2,
F = - (Yd + a32c2)a32b^ + (cß - ac2 - ßc2)ßb\ + 6(a + ß - c - 6aßcbj)b2 + + j[- (Yd + a32c2)a - (Yd - a32c + 2a32c2)ß + (y + a32)yc]b2 + ^- acbjja32 Jb3.
3. ВЫБОР КОЭФФИЦИЕНТОВ
Имеется большой произвол в выборе коэффициентов методов (2), (3), поскольку уравнений, определяющих порядок методов, меньше, чем неизвестных коэффициентов. Мы предлагаем выбирать коэффициенты по методике, изложенной в [2], [3].
Рассмотрим модельное уравнение du/dt = Xu, X е С. Считая, что Lu = Xu, L1u = 2 u, L2u = 2 u,
применим к этому уравнению метод расщепления (2), (3). В результате Un + 1 будет выражено через Un с помощью некоторой функции R: Un + 1 = R(Xl)Un. Функция комплексного переменного R(z) называется функцией устойчивости метода. Областью устойчивости называется множество тех z, при которых |R(z)| < 1. Метод называется А(0)-устойчивым, если его область устойчивости содержит интервал 0). Зафиксируем некоторое положительное число Z.
Определение. 5-Стадийный численный метод (2), (3) называется Z-оптималъным, если среди всех А(0)-устойчивых методов вида (2), (3), имеющих порядок аппроксимации 5, он доставляет минимум функционалу
J = max \R(z) - exp (z)|.
z е [-Z; 0]
При 5 = 2 функция устойчивости R(z) является функцией четырех переменных a11, a22, с1, c2, а при 5 = 3 - пяти переменных a, ß, Y, c, d. Таким образом, для нахождения оптимальных методов требуется найти условный экстремум недифференцируемой функции многих переменных.
При 5 = 2 и Z = 1 численное решение задачи минимизации функционала J приводит к методу, который назовем метод 2:
3 100 102 125
173
£И 0
250
14933 249
117500 500 158 235 393 393
При s = 3 и Z = 1 аналогично получаем метод 3:
0.233 0.170444235920 0.705 0.497 0 0.038623264130 0.287 0 -0.312350451890 0.728616156642 0.301
-0.534248810957 0.855224972756 0.679023838201
4. ТЕСТОВЫЕ РАСЧЕТЫ В качестве тестовых примеров рассмотрим две задачи, связанные с уравнением
= vf^-U + ^-U) + ft, x, у), v = const > 0, (x, у) e V = [0; 1] x [0; 1], 0 < t < 1.
dt x2 Эу2)
В этих задачах начальные, граничные условия и функция f(t, x, у) подобраны так, чтобы точным решением первой задачи являлась функция
u(t, x, у) =
а второй - функция
, , I 2 x + 2y-t
1 + exP I-4f—
и(г, х, у) = ехр[-10(х - уг)2]. Численные расчеты проводились на равномерной пространственной сетке с шагом Ах = Лу = 1/М
и с постоянным временнЫм шагом т = 1/N. Производные второго порядка L1u = V —-, L2u = V —-
д х ду
были аппроксимированы обычным образом на трехточечном шаблоне. Условия на границе S квадрата V для оптимальных методов вычислялись по формулам Y¡|S = и(г, х, у) |(ху)е = ( + с т,
= и(, х, у ) |(х,у)е , = + с,т .
Указанные выше методы сравнивались с методом стабилизации (см. [4, формула (3.22)], метод ST) и методом предиктор-корректора (см. [4, формула (3.42)], метод РК). Эти методы имеют второй порядок аппроксимации по времени, и граничные условия выбирались так же, как в [3].
Результаты расчетов приведены в таблице. В ней квадратичная погрешность численных расчетов определяется известной формулой
e2 =
ЛхЛу£ £ [и(1, Xi, yj) - UNj ]2,
i = 1 j = 1
где и(1, х, у) и и^ - соответственно, точное и приближенное значения решения в момент времени г = 1 в расчетных точках (х, у).
Таблица
Численные методы Задача 1 Задача 2
v = 0.01 v = 0.001 v = 0.01 v = 0.001
M = 20, N = 20 M = 40, N = 40 M = 20, N = 20 M = 40, N = 40
Метод 2 0.0128951 0.0782973 0.0006506 0.0000968
Метод 3 0.0050691 0.0907972 0.0005738 0.0000323
Метод ST 0.0137404 0.1601674 0.0007617 0.0001241
Метод PK 0.0156731 0.1595295 0.0067292 0.0015924
Как в случае жесткой задачи (V = 0.01), так и в случае нежесткой задачи (V = 0.001) оптимальные методы справляются с решением задач лучше, чем методы стабилизации и предиктор-корректора. Метод 3 имеет более высокий порядок аппроксимации, чем другие методы, и может применяться с большим временным шагом.
Отметим также, что метод 2 и метод 3 являются не только А(0)-устойчивыми, но и безусловно устойчивыми (см. [4]) в смысле теории устойчивости разностных схем уравнений математической физики.
СПИСОК ЛИТЕРАТУРЫ
1. Деккер К., Вервер Я. Устойчивость методов Рунге-Кутты для жестких нелинейных дифференциальных уравнений. М.: Мир, 1988.
2. Широбоков Н.В. Диагонально-неявные схемы Рунге-Кутты // Ж. вычисл. матем. и матем. физ. 2002. Т. 42. № 7. С. 1012-1017.
3. Широбоков Н.В. Расщепление четвертого порядка эволюционных уравнений в двумерном пространстве на основе диагонально-неявных методов // Ж. вычисл. матем. и матем. физ. 2004. Т. 44. № 10. С. 1824-1828.
4. Марчук ГИ. Основы вычислительной математики. М.: Наука, 1980.
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.