2. Конечноразностные методы
Одними из самых распространенных методов численного решения различных задач математической физики являются конечноразностные методы (методы конечных разностей). В различных вариантах этих методов в области определения искомых функций вводится сетка и решение ищется на этой сетке.
Для значений искомой сеточной функции (т.е. функции, заданной в узлах сетки) строится система скалярных уравнений, решение которой и служит таблицей приближенных значений решения исходной задачи. Один из способов построения этой системы скалярных уравнений состоит в приближенной замене производных, входящих в решаемое диф-ференциальное уравнение и в краевые условия, разностными соотноше-ниями, чем и объясняется название данного класса методов вычислитель-ной математики.2.1. Метод сеток.
2.1.1. Основные понятия и идеи метода. Изложим основные идеи метода сеток, рассмотрев его сначала применительно к простейшей линейной граничной задаче для обыкновенного дифференциального уравнения. Пусть при а <х<Ь рассматривается граничная задача вида
Au = -u" + q(x)u = f(x), хЄ(а,Ь), и(а) = u(b) = 0, (1)
где q(x) > qo = const > 0. Будем считать, что граничная задача (1) имеет единственное решение, это решение непрерывно на [а, Ь] и имеет непрерывные производные на этом отрезке до четвертого порядка включительно.
Метод сеток для решения граничной задачи (1), равно как и для многих других задач, состоит в следующем.
Область задания дифференциального уравнения (1), т.е. отрезок [а,Ь], заменяется некоторой дискретной (сеточной) областью. Это означает, что на отрезке [а,Ь\ выбирается некоторая система точек. Совокупность этих точек называется сеткой. Если положение каждой точки определяется по правилу Xk = a + kh, к = 0,1,... ,N, h = (b — a)/N, то сетку называют равномерной. Точки Хк есть узлы сетки.
Граничная задача (1) на множестве узлов, принадлежащих сетке, заменяется некоторой сеточной задачей.
Под термином сеточная задача понимают некоторые соотношения между приближенными значениями решения граничной задачи (1) в узлах сетки. В рассматриваемом случае это будет система линейных алгебраических уравнений.Полученная сеточная задача решается по какому-либо численному методу и тем самым находятся приближенные значения решения граничной задачи в узлах сепси. Это и является конечной целью метода сеток.
Можно указать на следующие вопросы как на основные в методе сеток.
а) Как заменить область задания дифференциальных уравнений, а в случае дифференциальных уравнений в частных производных еще и границу области некоторой сеточной областью?
б) Как заменить дифференциальное уравнение и граничные условия некоторыми сеточными соотношениями?
в) Будет ли полученная сеточная задача однозначно разрешимой, будет ли она устойчивой, сходящейся?
Поясним смысл последних двух понятий и дадим ответы на поставленные вопросы в применении к (1).
Построение разностной схемы. Выберем равномерную сетку.
xk = a + kh, k = 0,l,...,N, h = (-b~a\
Дифференциальное уравнение из (1) рассмотрим только во внутренних узлах сетки, т.е. будем рассматривать его в точках хк, к = 1 ,...,N - 1:
Аи\х=ч = -u"(xk) + q(xk)u(xk)=f(xk), k=l,2,...,N- 1.
Выразим производные, входящие в это равенство, через значения и(хк) в узлах сетки, пользуясь для них соответствующими конечноразностными представлениями:
(3)
4-і <4 <хк+й
и, ч u(xk+l) — 2u(xk) + ll(xk-\) (4)
и W = +гк (/г)'
г(4)ш- h2 JY(J4h v, ^ ,(4> ^ Г1 rk W -~Ї2 ^ k Xk-l u'(xN) = 3"Ы-4и(ХД,-,) + и^-2) + ^ 2/i С учетом конечноразностных представлений (применительно к (1) только для н"(хк)) получим Au\x=4 = Ai,u(xk) - Rk(h) = f{xk), где - / ч -u{xk-\) + 2u{xk)-u{xk+\) (4) Ahu(xk) = —г + q(xk)u(xk), Rk(h) = r\ '(h). Если для Rk{h) выполняется условие |К*(Л)| < Mh2, к= 1,2,... Пусть /і достаточно мало; тогда величиной Rk(h) можно пренебречь и получим A,,uk=f(xk), к= 1,2,...,W-1, (2) где при выполнении некоторых условий можно предположить, что и « « м(х,), і = 0,1,2,...,N. Вообще говоря, всегда и,фи(х,), и только при Rk(h) = 0 можно ожидать, что будет и, = и(х,), і = О,1,2,... ,N. Равенство (2) называют разностной схемой, аппроксимирующей уравнение Аи — Отметим еще, что (2) есть система линейных алгебраических уравнений, число таких уравнений N — 1, а матрица этой системы трехдиаго- нальная. Неизвестными являются ио,и\, ¦ • • Число этих неизвестных в системе равно N Обратимся к граничным условиям из (1). Используя граничные усло-вия, получим простейшие (в данной задаче) дополнительные уравнения (3) мо = 0, un — 0. Формулы (2), (3) образуют систему N+ 1 линейных алгебраических уравнений с неизвестными uo,ui,...,uN. Иногда (как в случае рассматривае-мой задачи) некоторые из этих неизвестных можно определить из «гра-ничных уравнений» типа (3), и в результате получаем уже задачу для N — 1 неизвестных мі,иг, ¦. ¦, un-і ¦ В других случаях можно из граничных уравнений выразить UO,UN через U\,U2,... ,UN-\. Подставив эти выражения в уравнения (3) вместо UQ,UN, снова приходим к системе уравнений для мі,м2,---,мдг_і. Однако заметим, что в последнем случае оператор А/, и выражения правых частей в (2) могут видоизменяться. Если теперь решить систему (2), (3) по какому-либо алгоритму, то после ЭТОГО МОЖНО принять u(Xk) ~ и k, k = 0,l,...,N. О разрешимости систем разностных уравнений В методе сеток приближенное решение рассматриваемой задачи вычисляется как решение системы разностных уравнений. Выясним условия разрешимости этих уравнений на примере системы (2), (3). Здесь из «граничных уравнений» имеем и о = О, и^ = 0. Поэтому будем рассматривать (2) для и і,..., uN~ і. называется матрицей с диагональным преобладанием величины 8 > 0, если |а„| > ф \ач\ + 5, і = 1,3,...,п. Для такой матрицы существует обратная А-1, причем ее норма, связанная с нормой = max; \xj\, удовлетворяет оценке ||А_!||<»< 1/8. Таким образом, если q(x) > qo = const > 0, то матрица системы (2) является матрицей с диагональным преобладанием; при любых {/(х*)} система (2) имеет единственное решение {у*}, причем тахЫ < — тах|/(х*)|. qo к Данное неравенство одновременно доказывает, что схема (2) устойчива к возможным ошибкам в значениях {/(**)} Решение системы (2) может быть найдено известным методом исключения Гаусса, который в применении к системам уравнений с трехдиаго- нальными матрицами называют также методом прогонки или методом факторизации Оценка погрешности и сходимость метода сеток Рассмотрим эти вопросы на примере задачи (1) Пусть tk = u(xk) - uk, e(h) - max |є*| 0 Ahek = Rk{h), k = 1,2, ,N — 1, Єо = Єуу=0 Считая, что точное решение (1) имеет ограниченные четвертные производные, получаем |/?*(Л)| < Mh2 А поскольку матрица системы для ошибок {є*} является с диагональным преобладанием, заключаем, что справедлива оценка 1 , тах|є*| < —А/Л О, k qo которая гарантирует также сходимость метода сеток в рассматриваемом случае задачи со скоростью 0{h2) Исследование основных вопросов обоснования метода сеток в более сложных случаях задачи осуществляется с привлечением различных специальных подходов и утверждений, развитых и полученных в теории данного метода (принцип максимума, теоремы сравнения и др) 2 1 2 Общие определения метода сеток Теорема сходимости Приведем основные определения и понятия теории метода сеток в более общих формулировках Пусть рассматривается некоторая краевая задача, которую запишем символическим равенством Аи = / в D, (4) где_исходные данные и точное решение и(х) определены в областях D и D = D U dD соответственно (В рассмотренной в п 2 1 1 задаче для обыкновенного дифференциального уравнения имеем D={a Пусть для приближенного вычисления решения дифференциальной краевой задачи (4), т.е. А/|М(Л>=/М. (5) Решение і/'1' = .. системы (5) определено на той же сет ке D/,, что и искомая сеточная функция [и]/,. Его значения Ug'\ uf\... С') /сч ...,uN' в точках XQ,X\,... ,XN последовательно вычисляются из (5) при и = 0,1,... ,ЛЛ Для краткости в уравнении (5) значок Л при не пишется. Способам построения и исследования сходящихся разностных схем (5) посвящено целое направление вычислительной математики, носящее название теория разностных схем. Придадим точный смысл понятию «сходящейся схемы» и требованию сходимости и^ —> [и]/,. Для этого рассмотрим линейное нормированное пространство функций, определенных на сетке D/,. Норма сеточной функции и^ Є ?//, есть неотрица тельное число, которое принимается за меру отклонения функции и^ от тождественного нуля. Норма может быть определена различными способами. Можно, например, принять за норму функции точную верхнюю грань модуля ее значений в точках сетки, положив ||«W||„A=maxK|. Если D С R", то часто используют также норму, определенную равенством ll«(A)lk = (ft"Z«*)l/2, 4 k=0 где N + 1 — число узлов сетки и неизвестных uo,ui,...,uN, определяемых разностной схемой. Эта норма аналогична норме ||м|| = (^J для функций и(х), интегрируемых на D с квадратом. После того как введено нормированное пространство U/,, приобретает смысл понятие отклонения одной функции от другой. Если а"'), произвольные сеточные функции ИЗ U/,, то мерой их отклонения друг от друга считается норма их разности, т.е. число ЦА^ — Дадим теперь строгое определение сходящейся разностной схемы. Система (5), как видим, зависит от А и должна быть выписана для всех тех А, для которых рассматриваются сетка D/, и сеточная функция []/,. Таким образом, разностная краевая задача (5) это не одна система, а семейство систем, зависящее от параметра А. задачи (5), принадлежащее пространству U/,. Говорят, что решение иС1' разностной краевой задачи (5) при измельчении сетки сходится к решению дифференциальной краевой задачи (4), если ||[«],, — u^\\uh —> О, h —> 0. Если сверх того выполнено неравенство ||[м]Л — и^\\цк < chk, где с > 0, к > 0 — некоторые постоянные, не зависящие от h, то говорят, что имеет место сходимость порядка hk, или что разностная схема имеет k-й порядок точности. Обладание свойством сходимости является фундаментальным требованием, которое предъявляется к разностной схеме (5) для численного решения дифференциальной краевой задачи (4). Если оно имеет место, то с помощью разностной схемы (5) можно вычислить решение [и]/, с любой наперед заданной точностью, выбирая для этого h достаточно малым. Придадим теперь точный смысл понятию «аппроксимация задачи (4) на решении к(х)» разностной схемой (5). Условимся считать, что правые части уравнений, записанных символи-ческим равенством (5), являются компонентами вектора fW из некоторого линейного нормированного пространства F/,. Тогда на А/, можно смотреть как на оператор, ставящий в соответствие каждой сеточной функции иМ из Uи некоторый элемент /С'' из F/,. В таком случае имеет смысл выражение А/, [«]/,, возникающее в результате применения оператора А к сеточной функции [и]/, из Uh и являющееся элементом пространства /•},. Невязка 5/М =Л/,[м]/, - /(А' принадлежит пространству Fh как разность двух элементов этого пространства. Под величиной невязки следует пони- мать \\8fW\\FH. Определение 1. Будем говорить, что разностная схема А/,и^ = = /М аппроксимирует задачу Аи = / на решении и, если Нб/^Ц^ —> 0 при h —> 0. Если сверх того имеет место неравенство Цб/^'^Ц^ < chk, где с > 0, к > 0 — некоторые постоянные, то будем говорить, что имеет место аппроксимация порядка hk, или порядка к относительно величины h. Обратимся к общему определению устойчивости схемы (5), которая аппроксимирует задачу (4) на решении и(х) с некоторым порядком hk. Определение 2. Будем называть разностную схему (5) устойчивой, если существуют числа ho > 0, 8 > 0 такие, что при любом /г < /г0 и любом є« Є Fh, ||є(/,)||^ < 8 разностная задача АйгЮ =/<*>+?<*>, полученная из задачи (5) добавлением к правой части возмущения сМ имеет одно и только одно решение г''1), причем это решение отклоняется от решения иМ невозмущенной задачи (5) на сеточную функцию г^ — — «(А), удовлетворяющую оценке где с\ — некоторая постоянная, не зависящая от h. В частности, последнее неравенство означает, что малое возмущение правой части разностной схемы (5) вызывает равномерно относительно h малое возмущение z^ — и^ решения. Пусть оператор А/,, отображающий Ui, в F/,, линейный. Тогда приве-денное выше определение устойчивости равносильно следующему опре-делению. Определение 3. Будем называть разностную схему (5) с линейным оператором Ah устойчивой, если при любом /М Є F/, уравнение Ани^ = = /С1' имеет единственное решение и^> Є U/„ причем ll^'lk^ll/^lk, (6) где с j — некоторая постоянная, не зависящая от h. Можно показать, что определения 2 и 3 эквивалентны. Теперь докажем один из фундаментальных результатов теории разностных схем, а именно докажем, что из аппроксимации и устойчивости следует сходимость. Теорема сходимости. Пусть разностная схема AhuW = /М ап-проксимирует задачу Аи = f на решении с порядком hk и устойчива. Тогда решение «М разностной задачи А/,и^ = сходится к [и]/„ причем имеет место оценка \\[u)n-uW\\uh Доказательство. Положим є^ = 5, [и]/, = z^. Тогда из определения 2 имеем Учитывая оценку из определения 1, сразу получаем доказываемое неравенство. ? В заключение подчеркнем, что схема доказательства сходимости ре-шения задачи (5) к решению задачи (4) путем проверки аппроксимации и устойчивости носит общий характер. При этом под Аи = f можно понимать любое функциональное уравнение, на основе которого строится «разностная задача» (5). 2.1.3. Метод сеток для дифференциальных уравнений в частных производных. Дифференциальные уравнения в частных производных имеют широкие приложения в математической физике, гидродинамике, акустике и других областях знаний". В большинстве своем такие уравнения в явном виде не решаются. Поэтому широкое распространение получили методы приближенного их решения, в частности метод сеток. Построение различных схем метода сеток в случае уравнений в частных производных зависит от типа уравнений и вида граничных (начальных) условий, присоединяемых к уравнениям. Частными примерами данных уравнений являются следующие: уравнение Пуассона (эллиптическое уравнение) /д 2и д 2и\ уравнение теплопроводности (параболическое уравнение) ди д2и . волновое уравнение (гиперболическое уравнение) д2и д2и - а? =/('.')• Изложим некоторые подходы к построению разностных схем применительно к этим уравнениям. Разностные схемы для уравнения параболического типа. Рассмотрим задачу Коши для уравнения теплопроводности ди д2и . . - = ^ + „(,,0, — <*<+-, ,>0) (7) и(х, 0) = i|/(x), < X < Будем считать, что задача (7) имеет единственное решение u(x,t), непре- Э'и дк и рывное вместе СО СВОИМИ производными гг—, і = 1,2, и т, к — 1,2,3,4. at' ох* Запишем задачу (7) в виде Au=f. Для этого достаточно положить Г ди д2и АН=\ а" а-<'<+-• '>0> [ И(Х, 0), —°° < j: < T = 0; ~ [ V(JT), < х < t = 0. Будем далее считать, что t изменяется в пределах 0 < t < Т < В рассматриваемом случае D = {—°° < х < 0 < / < Г < Г — объединение прямых t = 0 и t = Т. Выберем прямоугольную сетку и заменим D = D U Г сеточной областью ?)/,. К области ?>/, отнесем совокупность точек (узлов) {xm,t„), координаты которых определяются по правилу х,„ = mh, /н = 0,±1,±2,..., /г > 0, t„ = пт, и = 0,1,...,ЛГ, т > 0, NX<{N+ 1)т. Заменим задачу Аи — / разностной схемой вида Аи^ = /W. Обозначим через и(х,„,(„) точное значение решения задачи Аи ='/ в узле (xm,t„),? через и','п — соответствующее приближенное сеточное значение. Имеем А"\(Wn) = Э' j 'л ) (-Wn) Um.'n) d« -— , m = 0,±l,..., n=l,...,N, (*т,'пі ф(*.0 vW , m = 0, ±1,±2,..., R= 1,2 ЛЛ, э2« I (-Wn) Для замены выражений (дгт,гя) Эх2 Эг разностными отношениями (•Wn) воспользуемся формулами численного дифференцирования. Имеем и(-*СТ,'я+і)-и(Хст,'«) Т д2и т 2 dt2 (*т,'л) /ц+l) ~ u(xm,t„-l) I2 д3U 6 dt3 (*т,'л) 2Т м(хш+і, /„) - 2u(x„„ t„) + u(xm-i,t„) h2 d 4« /г2 12 Эх4 (¦Wn) _ u(xm,t„)-u(xmft„-i) ^ тЭ2Ц (W,) _ Т 2 Э/2 Эй д t ди dt ди dt д2_и дх2 m > 'и < ^ < '«+1> (•Wn ) (2) (2). > 'и-1 '» < '»> (•WA ) (3) (3) I 'n-1 < <'H+1, Рассмотрим следующую двухслойную разностную аппроксимацию: A"l(;Wn) = ц(х,п,/„+|)-ц(хш,/,,) ц(х,„+|, t„) - 2ц(х„„ /„) + г<(хт-і, /„) _ (Л) г Тт» > = Т /г2 . г<(х„„0), я = 0, где (wi") 12 Эх4 (Л) _ _ X Э2М 2dt2 Введем следующее обозначение: fW _ | _ / и запишем разностную схему для задачи Аи = /: где разностный оператор определяется по правилу АІ'У") = ( ;/"+' _ и" - ?»/' 4- и" ит ит т+1 zum + "m-1 Т А2 { т = 0,±1,±2,..., л = 0,1,...,N-1. Аналогично можно получить разностную схему вида Л<2>м(") =/("), (9) где ит ит _ "m+l ^"w-l , . /г2 ' <-(ft) _ ! ф(Л"ш, t„j, и0 "in і ш = 0, ±1,±2,..., я = 0,1,...,/V — 1, Замечаем, что 1(1 Vi где Г (Л) Г тЭ2" (wi21) 12 5л4 5(')/(") = ] 5(2)//.)= J 2^2" I 10. Выясним порядок аппроксимации разностных схем (8) и (9). В качестве /*), возьмем линейное множество всех пар ограниченных функций _ / <с ~ I р Норму в F/, определим по правилу ||Л = тахК| + тах|рш|. »1,Н Если шахт|„ |aJJ,| или тах,„ |Р„,| не достигается, то будем под нормой понимать величину ||g№|| = supm),|a^| + supjpm|. Пусть т = rhs, r,s — некоторые положительные числа. Предположим, Э2и Э4и ЧТО для -г-у И —7 верны оценки О/- 0X4 э2« д4и <Л/4- шах dt2 Эх4 <Мі, шах М Тогда легко получить ||б(Ч/(Л)1к* = шах|rm,)(Л)| < С-Мг + '^-щ) Л'\ m,i і \2 12 / ||8<2>/(/,)11 = max|r„,j(A)| < (Г-М2 + -h*. in,и \2 12 / В случае схемы (8) здесь можно взять s = 2, в случае схемы (9) можно взять 5 = 1. Из этих формул следует, что разностные схемы (8), (9) аппрок-симируют задачу Аи = f на и(х,у) с погрешностью порядка s (1 < s <2) относительно h. Разностная схема (8) позволяет по значениям решения на нулевом слое, т.е. по значениям и%, т = 0,±1,±2,..., вычислить значения на первом слое ихпо т = 0,±1,±2,... Для этого достаточно в (8) положить п = 0 и произвести вычисления, носящие рекурсивный характер. Потом по значениям и}„ можно аналогично при п — 1 вычислить значения и2 и т. д. В силу таких вычислительных свойств разностную схему (8) называют явной. Разностная схема (9) упомянутыми выше свойствами не обладает. Действительно, если в (9) положить п = 0, то в левой части полученной формулы будет линейная комбинация из значений ujn_j, и)ю в правой части будут значения известной функции ф(х,„,0) и ц(хт). Для вычисления значений на первом слое ..., и1_2, мір "о> м!> м2' •¦¦ в ЭТ0М случае необходимо решать бесконечную систему линейных уравнений. По этой причине разностную схему (9) называют неявной. Изучим вопрос устойчивости схем (9), (10), определив норму в пространстве ?/Л по правилу ||«<%,=тахК|. /п,и Рассмотрим разностную схему (8). Выясним, при каких значениях г, т = = rh2, возможна устойчивость этой схемы. Для доказательства устойчивости надо показать, что разностная схема однозначно разрешима и при любых gw = / . gW є Fln ( Рш, имеет место оценка l|z(A)lk Разностная схема (8), явная, и ее однозначная разрешимость очевидна. Перепишем формулу A^z^ = g^ в виде С1 = г{&+, +4-.) + (1 -2г)? + т<С z° = р„„ т = 0,±1,±2,..., « = 0,1,2, ...,N—1. Из этих уравнений при выполнении ограничения вида г = (т//г2) < 1/2 получаем неравенство шах|^'+1| < max\z!'m\ + ттах|а„| ш т т,н — принцип максимума. Полагая п = 0,1,... ,N — 1 и суммируя получаемые соотношения, получим maxl^l < max|p,„| + /Vxmax|a^| < тах|рш| + Гтах|а;|,| < т т т,п т іп,п < тах(1,Г) (maxК| +тахіні) = М ¦ |U(ft)||fft, \ т,п т / где обозначено М = тах(1, Т); М = 1, если Т < 1, и М — Т, если Т > 1. Следовательно, ІІг(Л)Ік<лні«(й)Ік- Таким образом, схема (8) устойчива при выполнении ограничения г < 0,5. Отметим, что это ограничение налагает жесткие ограничения на выбор т: т < 0,5/і2, и оно приводит к тому, что если необходимо сохранить устойчивость, при вычислениях по схеме (8) шаг по времени приходится выбирать очень малым. Обратимся к разностной схеме (9) и перепишем ее в виде -КС1. + <,-,) + (1+2 /-К+1 = <,+T(p(w„), (10) м° =V(xm), m = 0,±1,±2,..., л = 0,1,..., ЛГ — 1. Данная схема представляет собой бесконечную систему линейных уравнений относительно неизвестных ... к'ір MQ, м}, ... Решение таких систем является сложной и трудоемкой задачей, поэтому разностные схемы (9) неудобны для задач Коши на бесконечных отрезках и применяются редко. Однако если отрезок оси х, на котором рассматривается задача Коши, конечен: а <х <Ь, Ь — а<К, а на прямых х = а и х = b дополнительно заданы некоторые ограничения на решение и(х, t), то разностные схемы вида (9) оказываются весьма эффективными. В частности, такие схемы являются абсолютно устойчивыми, т. е. устойчивыми при любых значениях г = x/h2. Если, например, на отрезках прямых х — а и х — b заданы условия и(а, t) = Yo(')> ') = Yi (')' т0 вид системы (10) при п = 0 существенно изменится: -К",',,+ 1 +"ш_і) + (1 +2r)«i = V(*,„) +T м0 — Yo('i), "м = Yi('i)> т— 1,2,... ,Л/ — 1, li={b-a)/M. Формулы (11) представляют собой систему М+ 1 линейных алгебраических уравнений относительно гг', ..., и[м. Эту систему можно решить по методу прогонки, а затем определить все и"т при других значениях п. Приведенные рассуждения показывают, что реализация неявных разностных схем требует больших вычислительных затрат для вычисления решения на одном временном слое, но таких временных слоев может быть немного из-за того, что в этом случае отсутствуют ограничения на соотношение т//г2. Если пользоваться явной разностной схемой, то вычисление решения на следующем слое осуществляется по рекурсионному правилу и связано с минимальными вычислительными затратами, однако из-за ограничения т/Л2 < 1 /2 число временных слоев в случае явных схем может быть существенно большим по сравнению с числом временных слоев для неявных схем. Отметим следующее утверждение о сходимости разностной схемы (8). Эта схема аппроксимирует задачу (7) на решении н(х,/) с погрешностью порядка 0(т + Л2) и устойчива при г < 1/2. Поэтому схема (8) будет сходящейся; при этом погрешность для приближенного решения будет величиной порядка 0(т + Л). Разностная схема для задачи Дирихле для уравнения Пуассона. Пусть в области D= {0<х<а, 0 < у < ?} рассматривается уравнение Пуассона с условием Дирихле на границе Г области D : м|г=°- (12) Будем считать, что (12) имеет единственное решение и(х, у) в D = DUT, Э4« Э4и И ЭТО решение имеет непрерывные В D производные - И zr—r. ох4 дуг Возьмем прямоугольную сетку, положив хт = mh, т = 0,1,... ,М, А = = а/М; y„=nl, п = 0,1,... ,N, I = b/N. Отнесем к множеству внутренних узлов D® все узлы, лежащие в D, а к множеству граничных узлов Г/, отнесем узлы, лежащие на Г. Пусть узел (//г, п) eoj. Замену дифференциального уравнения из (12) разностным будем осуществлять только во внутренних узлах. Имеем _/Э2_и \Эх2 Э2_и (хт,Уп) ду2 ) =f(xm,y„), (*т,Уп) и{хт+\,уп)-2и(хт,уп) + и(хт-\,уп) /і2 Э4Ц Л2 12 ' Эх4 = f{xm, у„), и(хт, у„+|) - 2ц(х,„, у„) + ц(х„„ у„_1) | /2 Э4и 12 Э/ _>¦<") (т,п) Є Dl, x,„_i < х,',,0 < xm+i, Уп— і < У»0 < Уп-1 • „ , Э4м(х, у) Э4и(х, у) Пусть функции —— и —^— ограничены по абсолютной величине в области D; тогда при достаточно малых Ли/ можно пренебречь членами, содержащими в качестве множителей /г2 и /2, и получим искомое разностное уравнение АьиМ =/">, (13) где Л С') — ~ц»'+1." ~ Um— 1 и —Ц»|,я+1 + 2ц,,, ,! — ит<11-І A,,U = ^ + /2 ' /(">=/(х„„у„), (т, n)eDl "(/|)|гЛ = О-? Здесь через итп обозначается приближенное сеточное значение реше-ния задачи (12) так, что ит„ к и(хт, у„). В силу определений и понятий, введенных в пп. 2.1.1,2.1.2, получаем где h2 Э4м /2 д4и . . - э*« Следовательно, при сделанных выше предположениях относительно т—т ахг и тг—г имеет место оценка дуг / .М-lN-l \\/2\ ii5/(A)i\F„ где С — постоянная, не зависящая от /г. Данная оценка означает, что схема (13) аппроксимирует задачу (12) на решении и(х,у) с погрешностью 0(h2 + l2). Заметим, что оператор Ah является симметричным и положительно определенным в вещественном гильбертовом пространстве F/, с нормой ll lta и скалярным произведением (и{>'\ = Х^Г,' hlumnvm,. Собственные значения этого оператора имеют вид im я) 4 nmh 4 -> плі A '"'"' = —sin2 — + -sin2—, \ Нетрудно видеть, что Если воспользоваться методом разложения по собственным функциям оператора Ah, то для решения задачи (13) несложно получить оценку ll^llr^C-ll/C'llr,, где С = l/X[l'l\ которая гарантирует устойчивость схемы (13). На основе утверждений об устойчивости и аппроксимации для (13) заключаем, что данная схема является сходящейся и справедлива оценка ||и<*> - [uj*||f4 <С-(/г2 + /2), где [«],, есть таблица значений точного решения на а постоянная С не зависит от h и /. Разностная схема для уравнения гиперболического тнпа. Рассмотрим задачу Коши для одного из простейших гиперболи- ческих уравнений — уравнения колебаний струны (волновое одномерное уравнение) А m д2и д2и rf . Au = Du=^-T-T1=/(x,t), -°o at ox (14) и|,=0 = ф(лО, 37f=0 = ^W' где f(x,t), ф(х), g(x) — заданные функции. Решение задачи (14) определяется формулой (решение Даламбера) . . . (*о+'о) ф XQ-?0 +ф XQ+/Q 1 Г /Ч. 1 Г Г /.. . к(*о,'о) = 2 2 J gMdx+2 Ц fdxdt> (дго-fo) О(дг0,г0) где D(xо, to) — треугольник плоскости Ох/, ограниченный осью г = 0 и прямыми («характеристиками») x+t =xo+to, x — t =хо — to, проходящими через (XQ, to). Из этой формулы следует, что u(x,t) определяется значениями ф и g в точках основания треугольника D(xo,to) и значениями / внутри и на контуре треугольника. Рассмотрим прямоугольную сетку xk = kh, tj = jx (k = 0, ± 1,...; j = = 0,1,2,...). Совокупность узлов (k, j) (-°° U/,«/, = ~2 ^ = W 0 = Здесь k = 0,±1,±2,...; j= 1,2,... Из уравнений (15) u„ определяется в нулевом и первом рядах сетки: мдо = Введем обозначения: x/h = r, Dr(xo,to) — треугольник, ограниченный осью t = 0 и прямыми rx + t = rxo+to, rx — t = rxQ — to, проходящими через (хо, to) ¦ Нетрудно видеть, что значение сеточного решения И/, системы (15) в узле (ко, jo) однозначно определяется значениями ф, g в узлах, лежащих на основании треугольника Dr(koh, jox), и значениями / в узлах внутри и на границе Dr(knh, jox). Пусть имеется последовательность сеток таких, что h -4 0, т = rh, г = const > 1, и некоторая точка (*o,'o) ('о > 0) является узлом каждой сетки последовательности. Если uh(xo,to) имеет предел ио(хо,'о) при h —»¦ 0, то последний определяется значениями /, ф, g в Dr(xo,to) и, вообще говоря, отличен от М(А'О,'О) — значения в (хо, to) решения задачи (14), определяемого значениями /, ф, g в треугольнике D(XQ, to), содержащем Dr(xo,to) как часть в случае г > 1. Поэтому в случае г > 1 последовательность и/,, вообще говоря, не сходится к точному решению задачи (14).? Пусть функции f(x, /), ф(х) и достаточно гладкие соответственно в полуплоскости / > 0 и на оси —< х < Тогда справедливо следующее утверждение: если г = const < 1, то при h —»¦ 0 последовательность и/, решений задачи (15) сходится равномерно к точному решению u(x,t) задачи (14) во всякой конечной области полуплоскости. Подходы к построению разностных схем, изложенные выше в применении к простейшим уравнениям в частных производных, могут быть рас-пространены на более сложные уравнения (на уравнения с переменными коэффициентами, многомерные и/или нелинейные уравнения, на системы уравнений). 2.2. Метод прямых. Основная идея метода прямых состоит в том, что производные по одним независимым переменным заменяются их приближенными выражениями через конечные разности, тогда как производные по остальным переменным остаются без изменения. Рассмотрим некоторые варианты данного метода в применении к линейным дифференциальным уравнениям второго порядка. 2.2.1. Метод прямых для уравнений параболического типа. Пусть в полуполосе 0 < х < 1, t > 0 поставлена следующая граничная задача: (16) u,=Au+f(x,t) и(х, 0) = ф(х), M(0,0 = VO(0. И(1,0 = VI(0- Здесь Аи = а(х, t)uxx + b(x, t)ux + c(x, t)u. Предполагаются существование и единственность достаточно гладкого решения рассматриваемой задачи. Будем искать его приближенно методом прямых в прямоугольнике П = {0<х<1, 0< / < Г<°о}. В зависимости от того, аппроксимируются производные по х или по /, получают продольный или поперечный вариант метода. Рассмотрим сначала поперечный его вариант. Обозначим через Пт решеточную область 0 < х < 1, t = t,, = т, п = = 1,2,...,N (Nx = Т), а под ее границей Гт будем понимать множество, состоящее из отрезка 0 <х < 1 прямой t = 0 и точек (0,//у), (І,/,,), п = = 1,2, ...,N. Рассмотрим задачу (16) на Пт + Гт. Производную и,(х, t) при t = t„, п = 1,2,...,N, заменим приближенно левосторонним разностным отно-шением и,(х, /„) Ra (u(x,t„) —u(x,t„-і))/т. Тогда для решеточной функции, задаваемой системой функций и„ = и,,(х), п = 0,1,...,N, приближенно представляющих функцию и(х, t) на отрезках 0 < х < 1 прямых t = t„, л=0,1,... ,N, в Пт + Гт можно поставить следующую дифференциально-разностную граничную задачу: u„(x) - U„-l(x) =a„ u„ + /„ (a) , T moW = Ф0), «„(0) = VO('H), и„(1) = где А»"п = а„(х) и" (х) + Ьп (х)и'„(х)+ сп (х) и„ (х), а„(х) = a(x,t„), bn{x) = b{x,t„), с„(х) = c(x,t„), f„(x) = f(x,t„). Задача (17) распадается на N граничных задач для обыкновенных дифференциальных уравнений второго порядка, так как функция ио(х) известна из начального условия, а неизвестные функции и„(х) % u(x,t„), п = 1,2,...,N, можно находить последовательно, решая при каждом значении п = 1,2,... одно обыкновенное дифференциальное уравнение второго порядка с соответствующими граничными условиями. Построенная вычислительная схема метода прямых достаточно четко поясняет геометрический смысл названия метода. В силу этого в многомерном случае рассматриваемый метод часто называют методом плоскостей или методом гиперплоскостей. Рассмотрим одну из продольных схем метода прямых для задачи (16). Здесь при построении вычислительных схем метода прямых аппроксимируют операцию дифференцирования не по времени t, а по пространственной переменной х. Для этого, подобно случаю поперечных схем, предварительно исходную область приближенно заменим решеточной областью, проведя прямые х = х„ = nh, n = 0,l,...,N (Nh = 1). На каждой из внутренних прямых х = х„, м = 1,2,..., N — 1, производные Uxx и их аппроксимируем через значения функции и на нескольких соседних прямых. Если для этих целей использовать простейшие симметричные выражения, то, подобно случаю поперечных схем метода, можно записать ип = а„ (t) ^ +М0 + + с„(г)и,,(')+/«('), '>0. «o(0 = Vo(0. «w(0=Vi(0> '>0, (18) и„( 0)=\|/(х„), л = 1,2,...,//- 1, где u„(t) ~ и(х„, t), a„(t) = а(х„, t), bn{t) = b(x„, t), c„(t) = c(x„, t), /„(/) = = f(x„,t). Тем самым нахождение приближенных значений u„(t) на прямых х = х,„ п = 1,2,...— 1, сводится к решению задачи Коши для системы N — 1 линейных обыкновенных дифференциальных уравнений первого порядка. Погрешность аппроксимации задачи (16) задачей (18) предопределяется неточностью замены производных Uxx, их соответствующими простейшими симметричными выражениями через значения функции и и будет, очевидно, величиной порядка Л2, если решение и(х, t) обладает достаточной гладкостью. Следует заметить, что в случае продольных схем при аппроксимации дифференциального уравнения на прямых, близких к граничным, обычно возникают трудности, связанные с нарушением однородности процесса приближения. В некоторых частных случаях, правда, эти трудности удается существенно ослабить за счет использования граничных данных. 2.2.2. Метод прямых для уравнений гиперболического типа. Пусть в полуполосе 0 < х < 1, t > 0 поставлена следующая граничная задача: ut,+g{x,t)u, = Аи + f(x,t), М(х,0)=ф(х), к,(х,0)=Ф(х), (19) «(0,0=Vo(0. M(1,0=VI(0- Здесь Аи имеет такой же вид, как и в (16). Будем считать, что а(х, t) > > а > 0. Существование и единственность достаточно гладкого решения задачи (19) предполагаются. Будем искать это решение методом прямых в прямоугольнике П = {0 < х < 1, 0 < г < Г < <»}. Как и в случае параболических уравнений, здесь также можно говорить о продольном и поперечном вариантах метода. Приведем сначала пример поперечной схемы. Прямоугольник П аппроксимируем решеточной областью Пт + Гт; при этом, в отличие от параболического случая (см. п. 2.2.1), к границе Гт рассматриваемой области отнесем также и интервал 0 < х < 1 прямой t —1\ = т. На каждой прямой t —1„+1 для п — 1,2,..., N — 1 производные и, и и„ заменим следующими левосторонними разностными отношениями первого и второго порядков соответственно: и(х, л,+ і) « , т u(x,tn+l)-2u(x,t„) + u(x,t„-\) ««(X, tn+1) « 5 . т Кроме того, на прямой t — 0 для аппроксимации производной и, исполь-зуем правостороннюю замену м,(х, 0) w (м(х, т) — и(х, 0))/т. Это позволит для решеточной функции, задаваемой системой функций и„(х) « и(х, /„), п = 1,2,... ,N — 1, поставить в области Пт + Гт следующую дифференциально-разностную граничную задачу: и„+і(х)-2ип(х) + и„-і(х) м„+і(х)-м„(х) 2 + 8"+' W = Ап+1"«+1 + /н+1 (х), ио(х) =ф(х), Ц'(Х);'<0(Х) =Ф(х), (20) и..+ і(0) =Щ{1„+1), и„+і(1) = Vi(r«+i). Здесь использованы обозначения, совершенно аналогичные применяемым в параболическом случае. Процесс решения аппроксимирующей задачи (20) очевидным образом может быть расчленен на следующие этапы. Сначала по заданной функции мо(х) = ф(х) непосредственно получают мі(х) = ф(х) +гФ(х). После этого последовательно находят функции мг(х), «з(х), -¦¦, UN(X); при этом нахождение каждой из них сводится к решению одного обыкновенного дифференциального уравнения второго порядка с граничными условиями первого рода. Построенная вычислительная схема характеризуется только первым порядком погрешности аппроксимации. Если при аппроксимации диффе-ренциального уравнения для замены производной и, использовать более точное равенство , Зк(дг, /,1+1) — 4и(х, t„) +и(х, /„-і) ut\x>'n+i) ~ ' погрешность которого есть 0(т2), то при малых значениях т > О, как правило, удается уменьшить часть погрешности аппроксимации исходного уравнения, зависящую от неточности замены первой производной по времени. Процесс построения поперечных схем метода прямых для гиперболических уравнений во многом схож со случаем уравнений параболического типа. Подобные аналогии еще более сильно проявляются в случае продольных схем метода. Так, если в рассматриваемой области вместо прямых, параллельных оси х, провести прямые дг = дг,,= nh, л = 0,1,..., N (Nh = 1), то для задачи (19) можно построить следующую продольную схему метода прямых: , ,ЧЛ ,п(Л > (л _ Л (Л и>.+ і(0 -2и„(г) + и„-і(г) "„ (0 = a„{t) + KO(0=VO(0. Kyv(r) = Vi(0. f >0, и.(0) = ф(х„), м,',(0) =Ф(х„), л= 1,2,...,W-1. При записи этой схемы использованы обозначения, вполне аналогичные принятым ранее. Аппроксимирующая задача (21) может быть истолкована как задача Коши для системы N — 1 линейных обыкновенных дифференциальных уравнений второго порядка относительно неизвестных функций u„(t), п = = 1,2,. ..,N—1, приближенно представляющих искомую функцию и(х, г) на прямых х = хп, л = 1,2,...,N - 1. 2.2.3. Метод прямых для эллиптических уравнений. Рассмотрим данный метод применительно к эллиптической задаче Дирихле частного вида. Пусть в прямоугольнике П = {0 < дг < X, 0 < >> < К} существует и единственно решение задачи Дирихле u>y + Au=f(x, v), и(х,0) = ф(х), и(х, У) = Ф(х), (22) «(o,j) = vM. ^М, где Аи = а(х, у)ихх + Ь(х, у)их + с(х, у)и, а(х, у) > а > 0. Предполагаем решение достаточно гладким и будем искать его приближенно методом прямых. Прямоугольник П аппроксимируют решеточной областью, проведя прямые y = y„ = nh, п = 0,1,...,N (Nh = Y). Производную иуу при у = у,„ я = 1,2, 1, выражаем приближенно через значения функции м, опи раясь на известную формулу ^ _ и(дг,;у„+і) -2и(х,у„) +u(x,y„-i) h2 д4и(х, у„ + в/г) иуу(х,у„)- /г2 12 Э/ где — 1 < 0 < 1. Тогда, используя обозначения, совершенно аналогичные применяемым ранее, можно записать '<«+1 С*) - 2и„(х) + u„-i(x) /г2 + A„un=f„(x), 0 < х < X, и0(х) = ф(лт), мдг(дг) = Ф(дг), 0 < дг < X, (23) w»(0)=\|/(y„), M„(X) = T(y„), п= l,2,...,N-l. Дифференциально-разностная задача (23) аппроксимирует исходную задачу Дирихле с погрешностью порядка Л2. Она сводится к решению системы обыкновенных дифференциальных уравнений второго порядка с граничными условиями первого рода. Вопрос построения вычислительных схем метода прямых обсуждался выше на примере ориентированных по осям прямоугольных областей. В случае областей более сложной формы применение этого метода может быть сопряжено с определенными осложнениями. В самом деле, проекции g,„ п = 0,1,...,N, на ось координат внутренних сечений области избранными прямыми совпадут при любом числе прямых только для соответствующим образом ориентированного прямоугольника. В общем случае п-с уравнение рассматриваемой вычислительной схемы метода, связывающее приближенные значения искомого решения на соседних прямых с номерами от п — р до n + q, без дополнительных предположений определено лишь на пересечении С/. = П gn+r -p<> Заметим также, что метод прямых применяется для численного решения дифференциальных уравнений более высоких порядков, нелинейных задач и систем уравнений. 13 В И Агошков и др 2.3. Метод сеток для интегральных уравнений (метод квадратур). Метод сеток в применении к интегральным уравнениям называют также методом квадратур (методом механических квадратур). Рассмотрим основные идеи этого метода для вычисления решения интегрального уравнения Аи = и(х) — J К(х, у)и(у) dy = f(x), 0 <х< I, с гладким ядром К(х,у), предполагая, что \К(х, у)| < р < 1. Зададим N, положим h= 1/N и будем искать таблицу [и]/, значений решения на сетке x„ = nh (я = 0, l,...,N). Для получения разностной схемы в равенстве и(хп) - JК{х„, у)и(у)dy = /(х„), л = 0,1,...,М, приближенно заменим интеграл суммой, пользуясь, например, квадратур-ной формулой трапеций погрешность которой есть величина 0(h2), если ф(у) является дважды непрерывно дифференцируемой функцией. После указанной замены интеграла получим м„ - Л ( К{хп, 0) uo + K{x„,h)u\ +... + K{x,„yN-i)uN-i + , 1) + ;; uN или (24) где go, 8 ь /СО = AhuW = m, m, . 8n, K(xn, 0) . rim, = un-h [ ¦uN] ¦ K(x,„l) 8n uo + K(xn,h)u\ + ...+ 2 V-H.-V-. 2 Построенная разностная схема A/,u^ — f^ аппроксимирует задачу Аи — = / на решении и со вторым порядком относительно шага Л, поскольку квадратурная формула трапеций имеет второй порядок точности. Пусть и('0 = (ко,«ь •.. ,UN) — какое-нибудь решение системы (24), и пусть us — один из тех компонентов решения, которые по модулю не меньше каждого из остальных > |и„|, п = 0,1,..., N. Из уравнения с? номером п = s системы (24) следует неравенство I f(xs\ = |и, - A (fei^) мо + K(xs, h) их + Поэтому схема (24) устойчива: Ц^Ік =шах|Ип| = к| < -i-|/(x,)| < -1-||/(")||,-А. п і — р і — р п В частности, при /(хп) = 0 отсюда следует, что система (24) не имеет нетривиальных решений, а следовательно, однозначно разрешима при любой правой части. Решение и^ задачи А/,иW = /М в силу теоремы о сходимости удовлетворяет неравенству IIИл - "(Л)1к = тах|и(тЛ) - и„| < ch2, где с — некоторая постоянная. Аналогично изложенному метод сеток применяется для решения многомерных интегральных уравнений.
Если же считать такие уравнения выполняющимися и вне G„, то даже для очень простых областей аппроксимирующая задача при некоторых N может оказаться неразрешимой. Чтобы преодолеть подобные затруднения, иногда предпочитают продолжать непрерывно внутрь области функцию, задающую граничные условия, а уравнения аппроксимирующей системы рассматривают лишь на соответствующих пересечениях G„, при этом вы-числительный алгоритм приходится строить специальным образом. Часто удается избежать подобных затруднений, разбивая рассматриваемую об-ласть на полосы не прямыми, а кривыми линиями, избранными с учетом формы области.