CC BY
137
Вычислительные технологии
Том 5, № 4, 2000
ЧИСЛЕННОЕ ИНТЕГРИРОВАНИЕ ДВУМЕРНЫХ КРАЕВЫХ ЗАДАЧ С БОЛЬШИМИ ГРАДИЕНТАМИ РЕШЕНИЯ *
Ю. В. Немировский, А. П. Янковский Институт теоретической и прикладной механики СО РАН
Новосибирск, Россия e-mail: shulgin@itam.nsc.ru
The idea of Runge — Kutta methods generalization to the two-dimensional case of approximated integration of initial-boundary problems corresponding to the partial differential equations has been put forward and tested. It is shown that some classical finite-difference schemes of integrating the equations of transfer and non-stationary unidimensional heat conductivity can be obtained as a corollary of such a generalization. New schemes with a high order of accuracy have been obtained for various problems of mathematical physics. The stability of these schemes has been proved and the calculation results for the problems with large solution gradients have been presented. There are given specific examples, which demonstrate that the classical lower-order schemes inadequately describe the solution of such problems, whereas the high-order schemes constructed by means of the suggested generalized Runge — Kutta methods yield good approximations to the accurate solutions.
Исследование многих задач механики и физики приводит к необходимости интегрирования дифференциальных уравнений в частных производных высоких порядков, причем зачастую порядки высших производных по разным переменным могут сильно отличаться. Кроме того, решения краевых задач, соответствующих этим уравнениям, обладают большими и быстроменяющимися градиентами. Так, например, осесимметричные колебания тонкой изотропной цилиндрической оболочки описываются уравнением [1]
где — прогиб; д — распределенная поперечная нагрузка; , , — некоторые постоянные, определяемые геометрией оболочки и материалом, из которого она изготовлена. Это уравнение имеет 4-й порядок по пространственной переменной и лишь 2-й по времени ¿. Если на кромке х = х0 оболочка защемлена (/ш(Ь,хо) = 0, /шх(1,х0) = 0), то, как нетрудно получить с помощью метода разделения переменных, функция в малой окрестности
* Материалы статьи докладывались на Международной школе-семинаре "Информационные технологии в задачах математического моделирования"(Новосибирск, 13 - 18 сентября 1998 г.). Работа выполнена при финансовой поддержке Президиума Сибирского отделения Российской академии наук (Приложение 1 к Постановлению Президиума СО РАН №473 от 18.12.97).
© Ю.В. Немировский, А. П. Янковский, 2000.
Введение
w,xxxx
+ A(w - Bq + Cwttt) = 0,
(1)
кромки x0 будет иметь по модулю большие производные по (краевой эффект), которые могут быстро изменяться по t. Следовательно, в малой окрестности заделки решение начально-краевой задачи, соответствующей уравнению (1), будет иметь большие градиенты. Другими уравнениями, решение которых обладает большими градиентами, являются уравнения краевого эффекта в теории оболочек [1]
AAAw - L(F) - q = 0, AAF + BL(w) = 0, L(w) = RäVii + R-lw22, (2)
где F — функция усилий; Ri — главные радиусы кривизны оболочки; A — оператор Лапласа; остальные функции и величины имеют прежний смысл. Величина области краевого эффекта имеет порядок А = R\Jh/R, где h, R — толщина и характерный радиус кривизны оболочки, и для тонких оболочек (h/R ~ 0.01) является малой величиной по сравнению с их характерными размерами. При этом в области краевого эффекта функции w, F и их производные могут существенно возрастать и быстро изменяться по модулю с приближением к кромке. Если рассматривать многослойные композитные конструкции, то дифференциальные уравнения, аналогичные (1), (2), будут иметь еще более высокие порядки, причем с увеличением слоев эти порядки возрастают, достигая 12 - 16-го и т.д. Кроме того, в многослойных оболочках краевые эффекты возникают в каждом слое, и в силу существенного отличия механических характеристик слоев эти эффекты в них будут существенно количественно и качественно отличаться [2]. Помимо этого, большие градиенты решения могут возникать в задачах с сингулярными возмущениями, например, при исследовании равнонапряженного армирования конструкций [3] и во многих других случаях.
Необходимость получения наиболее достоверной информации о поведении неизвестных функций в таких задачах определяется тем, что, как правило, именно в подобластях, где решения имеют большие градиенты, происходит их резкое возрастание по модулю, причем значения неизвестных функций в точках экстремумов в указанных подобластях могут по модулю на порядки превосходить соответствующие значения этих функций за пределами областей с большими градиентами решения [2, 3], что может вызывать негативные последствия. Так, например, резкая концентрация напряжений в некоторых точках конструкции может привести к разрушению ее материала именно в этих точках и, как следствие, к общей потере несущей способности конструкции; резкое изменение температурного поля может приводить к резкому и нежелательному изменению механических и теплофизических свойств изделия и т. д. Поэтому актуальной становится проблема разработки надежных устойчивых методов численного интегрирования таких краевых задач. Эти методы, кроме того, должны обладать высокими порядками точности, чтобы достоверно определять поведение неизвестных функций в областях с большими и быстроменяющимися градиентами решения.
В одномерном случае хорошо зарекомендовали себя неклассические методы Рунге — Кутты (МРК) [4] для численного интегрирования краевых и начальных задач с жесткими дифференциальными уравнениями, решения которых обладают указанными выше особенностями. Поэтому, естественно, возникает вопрос об обобщении МРК на многомерные случаи. Разработке этой идеи и ее апробации и посвящено настоящее исследование.
1. Обобщение методов Рунге — Кутты
Суть обобщения классического МРК изложим на примере одного уравнения в частных производных первого порядка. Пусть для некоторой функции y(t, x) поставлена на прямой
£ = ¿о задача Коши
У,г = / у(£о,х) = уо(х) (3)
где /, у0 — известные функции, имеющие такую гладкость по всем аргументам, какая потребуется в процессе рассуждений. Найдем приближенное решение начальной задачи (3) на прямой £ = ¿о + т, где т — шаг интегрирования. Согласно идее классического МРК [5], рассмотрим функцию
8
и(т,х) = Уо(х) + т> 6^(т,ж), 5 > 1, (4)
г=1
где
кг = / Но + йт, х, Уо + т ^ аук, (Уо + т ^ -ук^ П ; -у = ау = 0 при ^ > г; (5)
V з=1 з=1 /
сг, агз, -у — постоянные, подлежащие определению. Функцию и(т, х) будем рассматривать в качестве приближенного решения задачи Коши (3) и потребуем, чтобы она имела по т как можно больший порядок касания [5] (при произвольных / и т) к функции V(т, х) = у(£о + т,ж), т.е. функция <^(т, х) = V(т,х) — и(т,х) должна при т = 0 для наибольшего г удовлетворять условиям
р(0,х) = <^т (0, х) = <р'т (0, х) = ••• = <^;)(0,х) = 0, <^;+1)(0,х) = 0. (6)
При этом отклонение и(т, х) от V(т, х) будет равно [5]
Р(х) = тг+У;+%,х)/(г + 1)!, 0 < ?(х) < т.
В [6,7] показано, что система равенств (6) при условии произвольности / определяет систему уравнений, задающих постоянные сг, ау, -у. Из этой системы следует, что -у = ау (1 < г,< 5), а уравнения, определяющие величины сг, ау, совпадают с уравнениями для коэффициентов МРК в одномерном случае [5]. При этом функция и(т, х) на прямой £ = ¿о имеет по т 5-й порядок касания к V (т, х) (или к у (¿о + т, х)).
Известно [4, 5], что в одномерном случае все схемы численного решения задачи Коши по МРК автоматически переносятся на системы уравнений первого порядка. В связи с полученными выше результатами этот вывод естественно переносится на двумерный случай, если система разрешена относительно частных производных по одной из переменных, т. е. начальная задача имеет вид (3), где под у, /, уо нужно понимать вектор-функции. Более того, аналогичные результаты можно получить и в случае, когда функция / в (3) зависит не только от производной первого порядка у,х, но и от производных по более высоких порядков, например, у,хх и т. д., а также если функция у зависит не от двух, а большего количества переменных, например, у(£,х1,х2) и т.д.
и -к
Нетрудно видеть, что вычисление полных производных —— в правых частях равенств
-х
(5) сильно усложняется с увеличением ]. Однако обобщенный МРК (4), (5) позволяет
-к-
интегрировать задачу Коши (3) только приближенно, поэтому и производные в (5)
-х
целесообразно вычислять приближенно, используя формулы конечно-разностного дифференцирования. Естественно, можно использовать различные шаблоны для аппроксимации оператора дифференцирования по переменной . Так, при большем количестве узлов, входящих в шаблон [8], будем получать больший порядок аппроксимации оператора дифференцирования, что положительно скажется на точности решения задачи Коши (3). С
другой стороны, условия устойчивости изложенного метода или особенности конкретной задачи могут потребовать использования центральных или скошенных разностей для аппроксимации производной по . Например, в [9] приводятся схемы высокого порядка интегрирования по переменной ¿, полученные на базе классических МРК, причем операторы конечно-разностного дифференцирования по на разных стадиях (при разных г в (5)) предложенных схем изменяются, обеспечивая свойства монотонности решений с ограниченной полной вариацией (так называемые ТУБ схемы). Если на всех стадиях МРК операторы конечно-разностного дифференцирования по одинаковы, то предложенный выше метод редуцируется в метод прямых [4].
Известно [4], что явные МРК, к которым фактически относится и схема (4), (5), либо неустойчивы, либо накладывают неоправданно жесткие ограничения на шаги интегрирования. В одномерных случаях эти особенности устраняются при помощи использования неклассических МРК [4], которые являются неявными. Очевидно, что в связи со всеми полученными выше результатами естественно было бы обобщить неклассические МРК на двумерный случай по схеме, аналогичной (4), (5). Единственное отличие, возникающее при таком обобщении, будет заключаться лишь в том, что в (5) йгз- = агз- = 0 при ] > г (именно эта особенность приводит к неявным МРК). Для удобства дальнейших исследований обобщенные неклассические МРК целесообразно представить в форме, отличной от
(4), (5):
уп+1(х) = уп(х) + т ^ Ьг/Цп + сгт, х, Уг(х), у'(х)), (7)
г=1
Уг (х) = уп(х) + / (1п + с3 т,х,у (х),у- (х)), г = 1, 2,..., в, (8)
3=1
где уп(х) — приближенное решение задачи Коши (3) на линии ¿п = ¿0 + тп; п — номер слоя (п = 1, 2,...). (Если функция / в (3) не зависит от последнего аргумента у,х, то схемы (4), (5) и (7), (8) редуцируются в обычные МРК, при этом выступает в качестве параметра. Это обстоятельство и запись обобщенного МРК в форме (7), (8) наглядно подтверждают выводы, полученные выше на основе строгого математического анализа при обобщении классических МРК, о том, что коэффициенты Ьг, сг, агз- в (4), (5) и (7), (8) должны совпадать с коэффициентами обычных МРК [4], применяемых в одномерном случае.)
Равенства (8) определяют замкнутую систему дифференциальных уравнений относительно неизвестных функций у(х), зависящих только от одной переменной, т.е. применение обобщенных МРК позволяет как бы понизить размерность исходной задачи на единицу. В случае линейных или квазилинейных уравнений (8) они могут быть легко разрешены относительно производных у (х) и затем проинтегрированы по при определенных начальных или краевых условиях, которые естественно вытекают из краевых условий или условий периодичности решения исходной задачи (3). В некоторых случаях система (8) может быть проинтегрирована аналитически, в общем же случае ее можно проинтегрировать приближенно, применяя, например, обычные МРК. После решения задачи, соответствующей системе (8), с помощью (7) легко строится приближенное решение исходной задачи Коши (3) на (п + 1)-м слое в точках прямой Ьп+1 = ¿п + т и т. д.
Предложенный подход позволяет легко получать численные схемы начальной задачи (3) (и других начально-краевых задач) высоких порядков точности по , и, как будет показано ниже на конкретных примерах, эти схемы обладают устойчивостью и позволяют
получать надежные результаты при интегрировании задач с большими градиентами решения.
Следует подчеркнуть принципиальное отличие метода (7), (8) от широко известного метода прямых [4]. При использовании метода прямых начально-краевая задача сначала дискретизируется по пространственной переменной , по которой нужно как бы решить краевую (или граничную) часть задачи, и затем уже интегрируется каким-либо методом (возможно, аналитически) по временной переменой по которой ставятся начальные условия. При использовании же метода (7), (8), наоборот, задача сначала как бы дискретизируется по переменной £ на слои (расстояния между слоями — т — могут меняться), а лишь затем каким-либо способом (возможно, аналитически) интегрируется по краевая (или граничная) задача для системы (8). При этом в методе прямых интегрирование по £ производится одновременно по всем прямым Хк = х0 + кк (к = 0, ±1, ±2, ...), параллельных оси а, а на каждом этапе применения метода (7), (8) рассматриваются только два ближайших слоя £га+1, и количество фиктивных прямых ^ = + с^ т, параллельных оси , для в-стадийного МРК (7), (8) равно в, что в случае неявных методов приводит к существенному уменьшению количества неизвестных величин, определяемых при переходе от одного слоя к другому.
2. Решение некоторых задач математической физики
Исследуем теперь применение обобщенного МРК (7), (8) к решению некоторых задач математической физики. Сначала рассмотрим решение задачи Коши для уравнения переноса
y,t = аУ,х, y(0,x) = yo (x) (a = ±1), (9)
предполагая, что y0 — периодическая функция с периодом , т. е. y0(x + X) = y0(x).
Для интегрирования начальной задачи (9) по схеме (7), (8) воспользуемся коэффициентами одностадийного метода Радо IIA (далее будем просто говорить: "воспользуемся таким-то методом", используя названия, принятые в [4]) — неявного метода Эйлера, — тогда получим уравнения
yn+i(x) = yn (x) + raY' i(x), Yi(x) = yn(x) + raY' i(x),
откуда
yn+1(x) = Yi(x), Yl(x) = a(Yi(x) - yn(x))/r. (10)
Второе уравнение (10) может быть проинтегрировано аналитически [10] с учетом периодичности решения по , однако нас интересует численное решение задачи (9), поэтому для интегрирования (10) вновь используем неявный метод Эйлера, тогда окончательно получим
ym+i = (1 - aH)-i(ym+i - aHym+i), H = h/т, (11)
где m — номер узла в слое; h = xm+i — xm > 0 — расстояние между двумя соседними узлами в слое (шаг интегрирования по ) — может быть переменной по m величиной. Если для интегрирования уравнения (10) использовать одностадийный метод Радо IA (явный метод Эйлера), то вместо (11) будем иметь схему
ym+i = (1 + aH )-i(ymft + aHym).
(12)
Значения коэффициентов в правой части (11) при а = —1 ив (12) при а = 1 положительны и меньше единицы, что обеспечивает устойчивость этих схем при любых Н > 0. (В действительности схемы (11), (12) являются хорошо известными в теории разностных схем неявными, абсолютно устойчивыми схемами интегрирования уравнения переноса, построенными на трехточечных шаблонах [8].) Методы (11), (12) имеют первый порядок точности по обеим переменным.
Если для интегрирования задачи Коши (9) использовать по £ явный метод Эйлера, а затем по — либо явный, либо неявный метод Эйлера, в конечном итоге получим две явные схемы первого порядка, подобные (11), (12) и построенные на трехточечных шаблонах (п,т), (п,т + 1), (п + 1,т +1) или (п,т), (п,т + 1), (п + 1,т), но устойчивые лишь при 0 < Н-1 < 1 [8]. (Этот результат наглядно подтверждает тезис о том, что требование устойчивости явных МРК накладывает жесткие ограничения на шаги интегрирования, в частности, в рассматриваемом примере должно быть: 0 < т < к.)
Исследуем схемы более высоких порядков точности. Проинтегрируем задачу (9), используя одностадийный метод Гаусса — Лежандра (метод средней точки), тогда получим уравнения: уп+1(х) = уп(х) + таУ^ (х), У1 = уп + 0.5таУ1', откуда
Уп+1 (х) = 2У/(х) — уп(х), У/ (х) = 2а (У/(х) — уп(х))/т. (13)
Для интегрирования второго уравнения (13) применим двустадийный метод Лобатто ША (метод трапеций), тогда из (13) следует
УГ+1 = [(1 + аН )УГ — аН (ут + у^Л/а — аН), уД = 2УГ+1 — у^, (14) УГ = [(1 — аН )УГ+1 + аН (ут + ут+х)]/(1 + аН), у^1 = 2У? — у1. (15)
Если из вторых равенств (14), (15) исключить за счет первых равенств величины У1т+1, У,™ соответственно и во вновь полученных уравнениях заменить У™, У1т+1 на У™ = 0.5(у1+1 + у1), У1т+1 = 0.5(ут+\ + у'т+1), то в конечном итоге получим классические схемы второго порядка по обеим переменным, построенные на 4-точечном шаблоне [8]. При этом схема (14) будет устойчива, по крайней мере в энергетической норме [8], при а = —1, а схема (15) при а =1.
Построим теперь схему 4-го порядка по обеим переменным. Для этого по £ используем двустадийный метод Гаусса — Лежандра:
уп+1(х) = уп(х) + та(б1У1' (х) + &2У2 (х)), У (х) = уп(х) + та^У^ (х) + ^У*' (х)), где Ь = 0.5, аи = 0.25 (г =1, 2), «12 = 0.25 — v/3/6, «21 = 0.25 + ^/6. Отсюда:
У/(х) = 12а[а,,(Уг(х) — уп(х)) — аг,(У,(х) — уп(х))]/т, 3 = 3 — г, г =1, 2; (16)
уп+1(х) = уп(х) + л/3(У2(х) — У1(х)). (17)
Для интегрирования системы (16) применим трехстадийный метод Лобатто ША, который в векторной форме имеет вид
(У}т+1 = (У}т + Н[/(хт, (У}т) + 4/(хт + Н/2, (£}) + /(хт+1, (У}т+0]/6,
(18)
(^} = (У}т + Н[5/(хт, (У}т) + 8/(хт + Н/2, (£}) - /(хт+1, (У}т+1)]/24,
где (У}т = (У™, У2т}Т; (^} = (^1, ^2}т — вспомогательный вектор; / — вектор-функция правых частей уравнений (16), явно зависящая от в силу того, что правые части в (16) зависят от уп(х), которая на п-м слое предполагается уже известной. В развернутом виде система (18) громоздка, поэтому не будем ее приводить. Уравнения (18) образуют замкнутую систему линейных алгебраических уравнений (СЛАУ) относительно У1т+1, У2т+1, при а = —1 и У™, У2т, ^1, при а = 1 (такая зависимость выбора неизвестных в (18) от значений а определяется условиями устойчивости получающейся в конечном итоге схемы). Так как уравнения (16) являются линейными уравнениями с постоянными коэффициентами, то при любых т коэффициенты СЛАУ (18) будут одни и те же, т. е. требуется лишь один раз произвести обращение матрицы 4-го порядка системы (18), причем значения ^1, в дальнейшем нигде не потребуются и их можно не вычислять (они формально потребовались лишь для составления матрицы (18)). Если система (18) разрешена, то из (17) окончательно получим схемы
= Ут+1 + л/ЭДГ1 — УГ+1) (а = —1), С+1 = ут + ^2™ — УГ) (а = 1). (19)
Отметим, что при определении правых частей системы (18) необходимо вычислять значения функции уп(х) в промежуточных точках хт + Н/2. Эти значения можно приближенно определить либо с помощью интерполяционного полинома Лагранжа: уп(хт + Н/2) и
+ 9ут+1 — ут-1 — ут+2)/16 (но в этом случае шаг Н должен быть постоянным, иначе такая замена не будет иметь 4-й порядок точности), либо с помощью интерполяционного полинома Эрмита: уга(хт + Н/2) и + у£+1)/2 — — (у^У]/8 (здесь Н уже может быть переменным по т; штрих означает производную по в узловых точках хт, хт+1; при переходе от слоя к слою эти производные определяются с помощью (17): (Ут+1) = (Ут) + "\/3(У2 (хт) — У1 (хт)), где У вычисляются по формулам (16); следовательно, при использовании полинома Эрмита для интегрирования задачи Коши (9) с 4-м порядком точности по помимо начального условия у(0,х) = у0(х) необходимо задать дополнительное условие у,х(0,х) = у0(х)).
Авторам пока не удалось строго доказать устойчивость методов (18), (19), но серии численных экспериментов, проведенные ими, показали, что схемы эти абсолютно устойчивы и имеют 4-й порядок по обеим переменным.
Для сравнения решений, полученных с помощью схем (11), (14), (18), (19), рассмотрим уравнение переноса (9) при а = —1 и начальном условии
где функция у0(х) является частичной суммой ряда Фурье для периодической функции
к
у0(х) = ^8вт(Акх)/Ак, Ак = 2п(2к + 1), К = 0,1, 2,...
(20)
к=0
г(х)
— 1, —0.5 + I < х < /; + 1, I < х < 0.5 + /, I = 0, ±1, ±2, ..
(21)
Рис. 1. Начальные условия и решения для уравнения теплопроводности.
Рис. 2. Решения начально-краевой задачи для уравнения переноса.
с периодом = 1 (ломаная 1 на рис. 1). Точное решение задачи Коши (9), (20) имеет вид
к
у(Ь,х) = ^8вт(Лк(х — г))/\к.
(22)
к=0
Откуда следует, что на прямых х — Ь = 0.51 частные производные от по модулю равны:
Ы = |у,х I =8
к
^сов(п (2к +1)1)
к=0
= 8(К +1),
(23)
а модуль градиента решения в этих точках равен + 1). Следовательно, с увеличени-
ем модуль градиента возрастает и уже при = 4 его можно считать большим по сравнению со значениями модуля функции (22).
В качестве конкретного примера рассмотрим случай = 9 (график функции (20), соответствующей этому значению , изображен на рис. 1 кривой 2). При этом максимальные значения модуля градиента будут большими и превосходят число 112. На рис. 2 приведены графики точного и приближенных решений этой задачи на прямой х = —0.5
при 0 < t < 4. Кривая 1 соответствует точному решению; кривая 2 построена по схеме (11) при разбивке прямоугольной области G : 0 < t < 4, -0.5 < x < 0.5 на 201 слой по t (tn = пт, т = 4/N, п = 0,1, ... , N, N = 200) и каждого слоя на 101 узел по (xm = mh, h = 1/M, m = 0,1, ... , M, M = 100), величина наибольшего отклонения g от точного решения при этом составила по модулю 1.09; кривая 3 соответствует той же схеме (11) при N = 2000, M = 1000, причем здесь max |g| = 0.897. График 4 получен по схеме (14) при прежних значениях N, M; в этом случае max |g| = 0.39. А решение, построенное по схеме (18), (19) при тех же значениях N, M, отличается от точного менее, чем на 5.2 • 10-4, и поэтому соответствующий график визуально не отличается от кривой 1 на рис. 2.
Кривые 2, 3 на рис. 2 показывают, что схема (11) "отслеживает" только общую (в некотором усредненном смысле) тенденцию поведения решения задачи (9), (20), но не отражает локальные особенности решения (22), в частности, не выделяет локальных экстремумов, причем с увеличением t приближенное решение постепенно затухает, приближаясь к нулю. (Последняя особенность фактически заложена в схеме (11), где коэффициенты в правой части при а = —1 положительны и меньше единицы.) Решение, построенное по схеме (14) (кривая 4), уже отслеживает и локальное поведение решения (22). Однако с увеличением t это приближенное решение приводит к увеличению по модулю "пиковых" значений (экстремумов) в окрестностях точек t = 0.5/ справа и сглаживанию пиков в окрестностях тех же точек — слева. В [8] доказано, что схема (14) абсолютно устойчива в энергетической норме; приведенный же пример наглядно показывает, что в норме пространства , которая наиболее естественна при исследовании задач с большими градиентами решения, следует говорить, по-видимому, о р-устойчивости схемы (14). (Отметим, что во всех рассмотренных примерах было принято: H = h/т = 0.5; если же = 2, то решение, полученное по схеме (14), будет увеличивать пиковые значения в окрестностях точек t = 0.5/ слева и сглаживать пики справа. При = 1 все предложенные схемы приводят к точному решению задачи (9), (20).)
Таким образом, исследованные схемы показывают, что МРК позволяет легко строить устойчивые схемы высоких порядков точности, по крайней мере, для начальных и краевых задач с дифференциальными уравнениями в частных производных первого порядка, решения которых имеют большие градиенты. (Так, решение, построенное по схеме 4-го порядка (18), (19), при =35 в (20) и N = 2000, M = 1000 — градиенты по модулю превосходят значения 400, — отклоняется от точного менее, чем на 5 %.)
Обратимся теперь к другим задачам математической физики, которым соответствуют дифференциальные уравнения в частных производных более высоких порядков. Сначала рассмотрим начально-краевую задачу нестационарной одномерной теплопроводности
y,t = y,xx, y(0, x) = yo(x), y(t, ±0.5) = p±(t) (—0.5 < x < 0.5, t > 0). (24)
Интегрировать эту задачу по t можно по схеме, аналогичной (7), (8), с той лишь разницей, что функция f должна зависеть еще от одного аргумента Y .
Для решения задачи (24) используем одностадийный метод Радо IIA (неявный метод Эйлера), тогда получим: yn+1(x) = yn(x) + тУ1 (x), Y1 = yn + тУ1 . Откуда следует:
yn+1(x) = Yi (x), ту" (x) — Yi(x) = —yn(x). (25)
Исследуем устойчивость схемы (25), используя идею метода разделения переменных [8]. А именно, предположим:
yn(x) = an sin Afcx (Afc = 2n(2k + 1), k = 0,1, 2,...), (26)
где sin Afcx — собственная функция задачи теплопроводности (24). Тогда решение второго уравнения (25) при граничных условиях Yi(±0.5) = yn+1(±0.5) = 0 имеет вид: yn+1(x) = аП+1 sin Akx, где постоянные , &к связаны соотношением
аП+1 = R(rAk К, R(rAk) = (1 + tA2 )-1. (27)
Функцию R по аналогии с терминологией, принятой в [4], целесообразно назвать функцией устойчивости МРК для задачи теплопроводности. Так как при любых т > 0, Ak функция R удовлетворяет неравенствам 0 < R < 1, то любое возмущение коэффициента ак затухает с удалением от n-го слоя, а значит, метод (25) абсолютно устойчив по t. Граничную задачу для конечно-разностного аналога (КРА) второго уравнения (25)
т(ym- - 2ym+1+ym++\)/h2 - ym+1 = -ym (ym+1 = ym) (28)
можно проинтегрировать с помощью метода прогонки [8]. Прогонка для уравнения (28) устойчива, поэтому метод (25), (28) в целом устойчив. (Отметим, что уравнение (28) соответствует известной в теории разностных схем абсолютно устойчивой схеме с опережением [8], имеющей первый порядок точности по т и второй порядок по h.)
Используем для интегрирования задачи (24) одностадийный метод Гаусса — Лежандра, тогда получим: yn+1(x) = yn(x) + тУ/'(x), Y1 = yn + 0.5тУ1", откуда
yn+1(x) = 2Y(x) - yn(x), тУ/'(x) - 2Y(x) = — 2yn(x). (29)
Если из второго уравнения (29) за счет первого равенства исключить функцию У1, а затем вторые производные заменить их КРА, построенными на трехточечном шаблоне, то получим в итоге классическую симметричную шеститочечную схему [8], которая абсолютно устойчива и имеет второй порядок точности по т и h.
После анализа уравнений (29) по схеме (26), (27) получим, что функция устойчивости R, соответствующая им, имеет вид
R^Ak) = (2 - тAk)(2 + тAk)-1. (30)
Отсюда следует, что |R| < 1 для любых тA2 > 0, т. е. метод (29) абсолютно устойчив по t, но при некоторых тA2 > 0 может оказаться
-1 < R^Ak) < 0. (31)
Например, R^A2) = -0.33 при т = 0.1 и R^A2) = -0.28 при т = 0.01.
Из теории уравнений математической физики [11] известно, что решение однородного уравнения теплопроводности (24) при нулевых граничных условиях <^±(x) = 0 монотонно убывает по t. В случае же выполнения неравенств (31) коэффициенты ак+1, ак имеют разные знаки, т. е. последовательность ак немонотонна, но с увеличением n убывает по модулю. Следовательно, в этом случае приближенное решение не удовлетворяет условию монотонного убывания по t и искусственно порождает эффект фиктивных тепловых колебаний, не имеющих места в реальности. Очевидно, что эта особенность решения уравнений (29) переносится и на решение их КРА — симметричную шеститочечную схему. Насколько известно авторам, эта особенность симметричной шеститочечной схемы в классической теории разностных схем не была подмечена и объяснена. (Нетрудно видеть, что при использовании схемы с опережением (25), (28) при любых тAk > 0 функция устойчивости
удовлетворяет неравенствам 0 < Я < 1, т. е. схема (28) удовлетворяет требованию монотонности решения по ¿.)
Авторами были исследованы и другие МРК применительно к интегрированию задачи (24). Так, в силу линейности этой задачи двустадийные методы Лобатто ША, ШБ приводят в конечном итоге к симметричной шеститочечной схеме, а из всех МРК второго порядка по т только двустадийный метод Лобатто ШС удовлетворяет требованию монотонности решения по ¿. Уравнения этого метода имеют вид:
уп+1(х) = уп(х) + 0.5т(К + К,"), у = + 0.5т(К - У2"), 1 = + 0.5т(К + К,"), (32)
а функция устойчивости определяется равенством Я = (1 + + )—1, = тА|/2. Очевидно, что 0 < Я < 1 при любых тА| > 0.
Следует отметить, что схемы (29), (32) имеют 2-й порядок по т, но для реализации метода (29) нужно проинтегрировать на каждом слое только одно уравнение относительно 1 (второе уравнение (29)), а для реализации схемы (32) нужно проинтегрировать уже систему двух уравнений относительно У1,12 (второе и третье уравнения (32)). Такое усложнение метода (32) окупается выполнением требования монотонности решения по ¿. Матричная прогонка для КРА двух последних уравнений (32), как легко получить, будет устойчива, поэтому в целом этот метод абсолютно устойчив и имеет второй порядок по т и Л.
Из всех МРК 4-го порядка самым "удобным" для интегрирования задачи (24) является трехстадийный метод Лобатто ШБ:
уп+1(х) = уп(х) + т[1"(х) + 41"(х) + 1"(х)]/6, 11 (х) = Уп(х) + тЦ"(х) - 1"(х)]/6,
(33)
У(х) = уп(х) + т 1" (х) + 212 (х)]/6, 1з(х) = уп(х) + т 1" (х) + 512 (х)]/6, так как его функция устойчивости имеет вид
Я(тАк) = (1 - 37к + 372)(1 + 37к + 372)-1, = тА2/6 (34)
(т. е. метод абсолютно устойчив и удовлетворяет требованию монотонности решения по Ь : 0 < Я < 1) и при реализации схемы (33) для каждого п необходимо интегрировать граничную задачу для системы двух уравнений относительно 11 12 (2-е и 3-е уравнения (33)). Следовательно, реализация этого метода требует приблизительно тех же вычислительных затрат, что и метода (32), но схема (33) имеет точность на два порядка выше.
Трехстадийный метод Лобатто 111А имеет функцию устойчивости (34); для его реализации на каждом слое также необходимо интегрировать систему 2-х уравнений, но в отличие от (33) в правых частях этих уравнений предварительно нужно дважды продифференцировать функцию уп(х). Остальные МРК 4-го порядка применительно к задаче (24) не удовлетворяют требованию монотонности решения по ¿, т.е. при определенных значениях тА| > 0 их функции устойчивости удовлетворяют неравенствам (31).
Рассмотрим теперь конкретный пример приближенного интегрирования задачи теплопроводности (24) с большими градиентами решения. А именно, пусть начальные условия имеют вид (20), а граничные условия равны нулю. Точное решение такой задачи теплопроводности известно [11]:
к
у(Ь,х) = ^8exp(-АkЬ) 81П(айх)/Ак, к=0
max |y,t(0,x)| = max
— г
> max |8AK sin(AKx)| = 8AK. (36)
и оценка для производной y,t в начальный момент времени имеет вид
к
-8 У^ Ak sin(Afcx) k=0
Следовательно, градиенты решения такой задачи будут превосходить по модулю значения 8Ak = 16n(2K +1), т. е. будут большими по сравнению со значениями функции (35). Так, например, max |y,t(0,x)| > 16п > 50 при K = 0 и max |y,t(0,x)| > 80п при K = 2.
На рис. 1 приведены графики, характеризующие приближенное решение этой задачи, полученное по схеме (33) (4-й порядок по т и второй порядок по h) при K = 2 и т = 0.01, h = 0.001. Кривая 3 соответствует начальному условию (20), кривые 4-7 соответствуют значениям времени t = 0.01 — 0.04 с шагом т, причем графики с номерами 4-7, изображенные сплошными линиями, соответствуют приближенному решению, а штриховыми — точному. Сопоставление точного и приближенного решений показывает, что наибольшее отклонение возникает на первом слое t = 0.01 (кривые 4), где градиенты наиболее велики. Максимальное отклонение на этом слое равно по модулю 0.091. С увеличением t градиенты решения быстро убывают, и приближенное решение удовлетворительно описывает поведение точного уже на всех последующих слоях (кривые 5-7). Если интегрировать уравнения (33) точно, то не получим принципиального улучшения приближенного решения задачи (24). Улучшить приближенное решение можно лишь за счет использования МРК более высоких порядков либо за счет уменьшения т. Так, при t = 0.001 кривые, характеризующие точное и приближенное решения, визуально вообще не отличаются.
Кривая 8 на рис. 1 соответствует приближенному решению задачи теплопроводности при t = т = 0.01, полученному по классической симметричной шеститочечной схеме (29). На форме этой кривой отчетливо проявились немонотонность схемы (29) и малый (второй) порядок точности по т, в силу чего график 8 ни качественно, ни количественно (наибольшее отклонение от точного решения составляет по модулю 0.28) не отражает поведения точного решения (штриховая линия 4). Следовательно, можно заключить, что классические схемы решения задачи теплопроводности малопригодны в задачах с большими градиентами решения.
Исследуем теперь вопрос интегрирования при помощи МРК задачи колебания струны:
y,tt = V,xx, y(0,x) = yo(x), y,t(0, x) = zo(x), y(t, ±0.5) = 0 (—0.5 < x < 0.5). (37)
В классической теории разностных схем для интегрирования такой задачи используют трехслойные схемы [8], получая точность 2-го порядка по т. Покажем, что с помощью МРК можно получит двухслойные схемы различных порядков по т. Для этого перепишем уравнение (37) в виде системы
Z,t = y,xx, y,t = Z. (38)
Для интегрирования системы (38) применим сначала одностадийный метод Радо IIA (неявный метод Эйлера): zn+l(x) = zn(x) + тУ1 (x), yn+1(x) = yn(x) + тZl(x), Yi(x) = yn(x) + тZ1(x), Z1(x) = zn(x) + тУ1 (x), откуда
yn+1 (x) = Y1(x), zn+1 (x) = Z1 (x), т2У1' (x) — Y1(x) = — yn(x) — тzn(x). (39)
Чтобы исследовать устойчивость схемы (39), вновь воспользуемся идеей метода разделения переменных. Предположим, что
yn(x) = ank sin Akx, zn(x) = bnk sin Akx, Ak = 2n(2k + 1),k = 0,1, 2,..., (40)
тогда из (39) с учетом граничных условий Y1 (±0.5) = 0 получим
yn+1(x) = a£+1 sin Afcx, zn+1(x) = b£+1 sin Afc x, (41)
где
аП+1 = К + Tbn)/(1 + T2Ak), b£+1 = (bn - Akran)/(1 + T2Ak), (42)
Используя выражения (42), можно записать
an-1T2 A2 an-2T2 A2 a0t2 A2 b0 т
an+1 = afc___afc T Afc___afc T Afc __afcT Afc |__bfcT
k _1 + t 2 A2 (1 + r2A2 )2 (1 + r2A2 )3 "' (1 + t 2 A2 )n+1 + (1 + t 2A2)n+1'
Отсюда видно, что возмущение коэффициентов а«, Ы с удалением от п-го слоя затухает, т. е. метод устойчив по ¿. Прогонка для КРА уравнения (39) устойчива [8], поэтому в целом метод устойчив и имеет первый порядок по т и второй порядок по Л.
Применим для интегрирования задачи (37), (38) одностадийный метод Гаусса — Ле-жандра: гга+1(х) = гга(х) + тУ/'(х), уга+1(х) = уп(х) + т^1(х), У1(х) = уп(х) + 0.5т^1(х), ^1(х) = гга(х) + 0.5тУ1 (х), откуда
уп+1 = 2У1 — у«, = 4(У1 — уп)/т — гга, т2 У," — 4У1 = —4уп — 2тгга (44)
Проведя анализ уравнений (44) по схеме (40) - (42), получим
„n+l _ (4 - т Afc)afc + rbfc 4(2bfc - TA2afc) m
% = 4 + t 2A| ' = 4 + т 2A| • (45)
Используя эти выражения, можно записать разложение, аналогичное (43) (в силу громоздкости не будем его приводить), из которого следует, что схема устойчива по t. Так как прогонка для КРА последнего уравнения (44) устойчива, то в целом метод устойчив и имеет второй порядок по т и h.
Для построения схемы 4-го порядка по т можно воспользоваться, например, трехста-дийным методом Лобатто IIIB:
zn+l = zn + т (Y" + 4Y2" + Y3" )/6, yn+l = yn + т (zi + 4Z2 + Za)/6,
Zi = zn + т(Y" - Y")/6, Z2 = zn + т(Y" + 2Y")/6, Z3 = zn + т(Y" + 5Y")/6,
Y = yn + т(Zi - Z2)/6, Y = yn + т(Zi + 2Z2)/6, Y3 = yn + т(Zi + 5Z2)/6,
откуда
т2Y" - 12Yi - 12Y2 = — 24yn - 6TZn, т2Y2" + 12Yi = 12yn, (46)
Zi = 2(2Yi + Y2 - 3yn)/t, Z2 = 2(Y> - Yi)/т, Z3 = 2(-4Yi + Y2 + 3уп)/т,
Y3 = -Yi + 2Y2, yn+i = yn + 2(Y2 - Yi), zn+i = zn + 12(yn - y)/t. (47)
Проанализировав решение системы (46), (47) по схеме (40) - (42), получим
an+1 _ + ! 2 (12 - t2 A2)(4an + Tbn) - 2(24 + t2 A2R % % + 12 144 + t 2 A2 (12 + t 2 A2 )
bn+1 = bn +
t
n 12
n
n_ 144ап + 6t2a| (тьп + 4ап) afc 144 + t2a2 (12 + t2a2 )
Используя эти выражения, можно также записать разложение, подобное (43), из которого следует, что схема устойчива по ¿. Матричная прогонка для КРА уравнений (46) устойчива [8], поэтому метод в целом устойчив и имеет 4-й порядок точности по т.
Численные эксперименты, проведенные авторами, подтвердили устойчивость всех рассмотренных выше методов и их порядки точности.
Отметим, что можно интегрировать задачу (37), используя и другой подход. Для этого представим уравнение (37) в виде системы уравнений первого порядка, записанных в инвариантах Римана: 2,^ + 2,х = 0, — у,х = г. Эти уравнения являются уравнениями переноса, интегрировать которые можно последовательно и для которых применение МРК было исследовано выше.
Обратим внимание на то, что при исследовании вопросов устойчивости МРК применительно к решению задач теплопроводности и колебания струны были использованы представления (26), (40), где 8т(А^х) — собственные функции этих задач. Те же результаты можно получить, если представить уп(х), 2п(х) в (26), (40) в виде других собственных функций: есв^х), ^ = /2.
Во всех приведенных выше методах, касающихся решения задач теплопроводности и колебания струны, исследовалась фактически устойчивость по начальным данным. Если решать при помощи МРК эти задачи с неоднородными уравнениями, то при переходе с п-го слоя на (п + 1)-й в уравнениях МРК появятся дополнительные слагаемые — известные функции от — следствия неоднородности исходных уравнений. Эти функции можно разложить в ряды Фурье по собственным функциям задачи (вт(Акх), есв(^х)) и применить уже описанные выше алгоритмы исследования устойчивости. После чего получим устойчивость МРК по правым частям.
В заключение отметим, что по схеме, подобной (37), (38), при помощи МРК можно численно интегрировать граничные задачи и для эллиптических уравнений. Так, в [7] при помощи явного метода 4-го порядка, аналогичного (4), (5), была решена следующая граничная задача:
(Лцу,4 + Л12У,х Ь + (Л21У,4 + Л22У,х ),х = 0, (49)
у(^0,х)= Ую(х), У (¿1, х) = уц (х) , у(£,х0)= Уж0 (х) , у(г,х!)= Уж1 (х), где Лу — переменные коэффициенты положительно определенной квадратичной формы. (Такая граничная задача возникает при решении задачи стационарной двумерной теплопроводности для неоднородной среды.)
Граничная задача (49) была сведена к начально-краевой задаче
У,4 = 2, = —Л-11[22,хЛ12 + У,ххЛ22 + ¿(Лц^ + Л^) + У,х(Л2М + Л22,х)],
(50)
У (¿0, х) = Ую(х), ¿(¿0,х) = 20 (х), у(£,х0)= Уж0 (х) , у(^,х1)= Ух1 (х) и проинтегрирована по схеме (4), (5). Функция 20 в (50) изначально неизвестна, а в начально-краевой задаче (50) осталось незадействованным граничное условие
У(£1,х) = У*1(х). (51)
При помощи метода пристрелки узловые значения функции 20 подбирались так, чтобы в узловых точках на линии £ = ¿1 выполнялось граничное условие (51).
При разбивке прямоугольной области С : ¿0 < £ < ¿1, х0 < х < х1 (¿1 — ¿0 = 1, х1 — х0 = 2п) равномерной сеткой на 75 слоев по £ и при введении на каждом слое 90 узлов по х точность интегрирования граничной задачи составила 10-11 значащих цифр; при
разбивке G на 20 слоев и при введении на каждом слое 40 узлов точность интегрирования составила 5-6 значащих цифр.
Так как при помощи метода пристрелки, совмещенного с МРК (4), (5), удалось проинтегрировать граничную задачу (49) с высокой точностью, то можно утверждать, что обобщенный МРК применительно к интегрированию начально-краевой задачи (50) является устойчивым.
Заключение
Проведенные в настоящей работе исследования позволяют сделать вывод о том, что обобщение метода Рунге — Кутты на двумерный случай открывает широкие возможности для быстрого и удобного построения устойчивых численных схем высоких порядков точности для дифференциальных уравнений в частных производных различных типов. Особенно важно, что это обобщение позволяет строить надежные численные схемы для задач с большими градиентами решения.
Список литературы
[1] РАБотнов Ю.В. Механика деформируемого твердого тела. Наука, М., 1979.
[2] Кобелев В.Н., Коварский Л. М., Тимофеев С. И. Расчет трехслойных конструкций: Справочник. Машиностроение, М., 1984.
[3] НЕмировский Ю. В., Янковский А. П. О некоторых свойствах решений плоских термоупругих задач рационального армирования композитных конструкций. Прикладная математика и механика, 61, В. 2, 1997, 312-321.
[4] Деккер К., Вервер Я. Устойчивость методов Рунге — Кутты для жестких нелинейных дифференциальных уравнений. Мир, М., 1988.
[5] Березин И. С., Жидков Н. П. Методы вычислений. Т. 2. Физматгиз, М., 1959.
[6] НЕмировский Ю.В., Янковский А. П. Численный метод решения плоских задач композитных конструкций с равнонапряженной арматурой. Тр. XIII Межресп. конф. по численным методам решения задач теории упругости и пластичности. Новосибирск, 1995, 130-139.
[7] Янковский А. П. Рациональное армирование пластин при продольном и поперечном нагружении: Дис. ... канд. физ.-мат. наук. Новосибирск, 1996.
[8] Самарский А. А. Теория разностных схем. Наука, М., 1989.
[9] Shu C.-W., Osher S. Efficient implementation of essentially non-oscillatory shock-capturing schemes. J. Comp. Phys., Vol. 77, No. 2, 1988, 439-471.
[10] Петровский И. Г. Лекции по теории обыкновенных дифференциальных уравнений. Гос. изд-во технико-теоретической литературы, М.-Л., 1948.
[11] Тихонов А. Н., Самарский А. А. Уравнения математической физики. Наука, М., 1977.
Поступила в редакцию 12 марта 1999 г.
