Численные методы. Соснин (2005) (1160462), страница 18
Текст из файла (страница 18)
Известно, что если выполняются условия k(x) > k0 > 0;q(x) > 0;β > 0,то решение задачи (5.18) существует и единственно.Задача (5.18) содержит уравнение параболического типа. Обычно такие уравнения возникают приисследовании распределения температуры в тонком стержне или в диффузионных процессах.Решение системы типа (5.18) проходит в несколько этапов.
Сначала ей сопоставляется дискретнаямодель, а на ее основе строится разностная схема. Существует несколько методов построения такихсхем, и первым мы рассмотрим интегро-интерполяционный метод. Название его происходит оттого, что в процессе построения соответствующей разностной схемы мы переходим от интегральныхсоотношений к интерполяционным уравнениям.Построение разностной схемы92Глава 5. ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ КРАЕВЫХ ЗАДАЧПерейдем к построению схемы.
Для начала введем на отрезке [0; l] равномерную сетку:l}Nωh = {xi = ih, i = 0, N , h =— очевидно, x0 = 0, xN = l. Теперь введем такое обозначение для средних точек между узлами сетки:xi± 21 = xi ±h.2Также обозначим ui = u(xi ) — значение искомой функции в узлах сетки, и W (x) = k(x)u0 (x).Применяя эти обозначения, фиксируем произвольное i ∈ [1; N −1] и проинтегрируем первое уравнениесистемы (5.18) по отрезку [xi− 12 ; xi+ 12 ] :xi+ 12Z((k(x)u0 (x))0 − q(x)u(x) + f (x)) dx = 0 ⇐⇒xi− 12xi+ 1xi+ 1ZZ2Wi+ 21 − Wi− 12 −2q(x)u(x) dx +xi− 1f (x) dx = 0,(5.19)xi− 122где Wi± 21 = W (xi± 12 ).xi+ 12ZПервый интеграл можно приблизить значением ui ·q(x) dx. Тогда (5.19) можно переписать какxi− 12приближенное равенство:xi+ 1xi+ 1ZZ2Wi+ 21 − Wi− 12 − ui ·2f (x) dx ≈ 0.q(x) dx +xi− 1(5.20)xi− 122Перейдем от интегральных выражений к линейным.
Для этого введем новые обозначения:xi+ 1ϕi =1hxi+ 12Zf (x) dx;di =xi− 1q(x) dx.(5.21)xi− 12Теперь заметим, что u0 (x) =2Z1h2W (x). Проинтегрировав это равенство на отрезке [xi ; xi+1 ], получим:k(x)xZi+1ui+1 − ui =W (x)dx.k(x)xiЗаменим это равенство приближенным:xZi+1ui+1 − ui ≈ Wi+ 21xidx.k(x)5.6. ИНТЕГРО-ИНТЕРПОЛЯЦИОННЫЙ МЕТОД ПОСТРОЕНИЯ РАЗНОСТНЫХ СХЕМ93Отсюда следует, что если обозначитьai+11=hxZi+1−1dx ,k(x)xiто будут справедливы приближенные равенства:ui+1 − ui; Wi+ 21 ≈ ai+1h W 1 ≈ a ui − ui−1 .ii− 2h(5.22)Воспользовавшись обозначениями (5.21) и (5.22), приближенное равенство (5.20) можно переписатьтак:ui − ui−1ui+1 − ui− ai− hdi ui + ϕi h ≈ 0.(5.23)ai+1hhОбозначим за yi такие числа, которые при подстановке в (5.23) вместо ui дают точное равенство:ai+1yi+1 − yiyi − yi−1− ai− hdi yi + ϕi h = 0.hh(5.24)Найденные из таких уравнений значения yi и будут считаться приближениями к ui .
Полученноеравенство и будет искомой разностной схемой, однако ее можно переписать и в более компактном виде,заметив, что первые две дроби — не что иное, как разностные производные назад. Обозначив их yx,i+1и yx,i , перепишем (5.24) в следующем виде:ai+1 yx,i+1 − ai yx,i − hdi yi + ϕi h = 0.Разделив это равенство на h, мы сможем объединить две разностные производные во вторую разностную производную вперед:(ayx )i+1 − (ayx )i≡ (ayx )x,i .hВ итоге мы получим такой вид разностной схемы:(ayx )x,i − di yi + ϕi = 0.(5.25)Мы имеем право написать такие равенства для i = 1, N − 1.
Их можно объединить в системулинейных (по построению) уравнений относительно yi . Она будет содержать N − 1 уравнение и N + 1неизвестное. Необходимые для однозначной разрешимости системы 2 уравнения добавим из краевыхусловий. Мы можем заменить в последнем равенстве в (5.18) u(l) на yN , тогда получим, чтоyN = µ2 .(5.26)Для получения последнего уравнения мы выполним те же самые действия, что и при выводе равенства (5.19), но интегрирование будем проводить на отрезке [0; h2 ].
Тогда можно получить такоеравенство:hhZ2Z2W 12 − W0 − u0 q(x) dx + f (x) dx ≈ 0.(5.27)0W120и W0 мы найдем, заменив в приближенных равенствахu1 − u0;hW0 ≈ βu0 − µ1 .W 12 ≈ a194Глава 5. ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ КРАЕВЫХ ЗАДАЧuk на yk и получив уравнения:y1 − y0;hW0 = βy0 − µ1 .W 21 = a1После этого, воспользовавшись обозначениямиhϕ0 =1h2hZ2f (x) dx,0d0 =1h2Z2q(x) dx,0мы приведем (5.27) к такому виду:hha1 yx,0 − βy0 + µ1 − d0 y0 + ϕ0 = 0 ⇐⇒22a1 yx,0 − βy0 = µ,d0 h2 ,(5.28)ϕ0 h2 .µ = µ1 +где β = β +Уравнения (5.25), (5.26) и (5.28) представляют собой окончательный вариант разностной схемы,полученной с использованием интегро-интерполяционного метода.Решение разностной схемыВторым шагом на пути решения краевой задачи численно становится выбор метода решения построенной схемы.
В нашем случае мы получили систему линейных уравнений, для которой методрешения и выбирается. Заметим, что (5.25) можно переписать так:1yi+1 − yiyi − yi−1ai+1− ai− di yi + ϕi = 0.hhhСобрав коэффициенты при yi , получим:Ai yi+1 − Ci yi + Bi yi−1 = −Fi ,i = 1, N − 1,где Ai = ai+1 , Bi = ai , Ci = ai + ai+1 + di h2 .Добавив к этим уравнениям уравнения (5.26) и (5.28), получим систему из N + 1 уравнения. Матрица, задающая эту систему уравнений, будет являться трехдиагональной, а такие системы обычнорешаются методом прогонки. Он применим, так как выполнены условия:Ai , Bi > 0, Ci > Bi + Ai— они дают существование и единственность yi , отвечающих уравнениям (5.25), (5.26) и (5.28).5.7Метод аппроксимации квадратичного функционалаЭто другой метод построения разностных схем.
Будем рассматривать его для задачи, схожей с (5.18),но с более простыми краевыми условиями:(k(x)u0 (x))0 − q(x)u(x) + f (x) = 0, 0 < x < 1;u(0) = u(1) = 0.Известно, что решение такой задачи эквивалентно поиску u, минимизирующих функционалZ1J[u] =0[k(x)(u0 (x))2 + q(x)u2 (x) − 2f (x)u(x)] dx.5.8. КОРРЕКТНОСТЬ РАЗНОСТНОЙ СХЕМЫ95Задав на отрезке [0; 1] равномерную сетку, разобьем интеграл по всему отрезку на сумму по подотрезкам:N ZxiXJ[u] =k(x)(u0 (x))2 + q(x)u2 (x) − 2f (x)u(x) dx.i=1xi−1Для упрощения поиска экстремума заменим обыкновенную производную на ее разностный аналог:N ZxiX2J[u] ≈k(x)ux,i+ q(x)u2 (x) − 2f (x)u(x) dx.i=1xi−11Воспользовавшись обозначением ai =hZxik(x) dx, получим:xi+1ZxiN X2ai ux,i h +J[u] ≈q(x)u2 (x) − 2f (x)u(x) dx .i=1xi−1Интегралы посчитаем по формуле трапеций, заменив всюду ui на yi — приближенные значения:J[u] ≈ J[y] = Jh (y0 , y1 , .
. . , yN ) == {y0 = yN = 0} =NX2h+ai yx,iN Xh22ai yx,ih + (qi yi2 − 2fi yi + qi−1 yi−1− 2fi−1 yi−1 ) =2i=1N−1X(qi yi2 − 2fi yi )h.i=1i=1Мы свели задачу о поиске элемента, минимизирующего функционал, к поиску чисел yi , доставляющих минимум функции многих переменных — при этом, правда, мы потеряли в точности.Необходимым условием экстремума будет равенство нулю всех частных производных:∂Jh= 0 ⇐⇒∂yi⇐⇒⇐⇒112ai+1 yx,i+1 (− · h) + 2ai yx,i ( · h) + (2qi yi − 2fi )h = 0 ⇐⇒hhai+1 yx,i+1 − ai yx,i− qi yi + fi = 0h(ayx )x,i − qi yi + fi = 0.Последнее уравнение — это уже часть итоговой разностной схемы (осталось добавить только краевые условия). Можно заметить, что она похожа на схему, возникающую в интегро-интерполяционномметоде, однако коэффициенты различны, да и свойства схем тоже довольно сильно различаются.5.8Корректность разностной схемыНапомним несколько определений.Определение.
Разностная схема аппроксимирует исходную дифференциальную задачу в точкеxi [на всей сетке], если погрешность аппроксимации в этой точке [соответственно, норма погрешностиаппроксимации] стремится к нулю [соответственно, тоже к нулю] с уменьшением h :hih→0h→0ψ(xi ) −−−→ 0 ||ψ||h = ||ψ||C(ωn ) = max |ψi | −→ 0 .i96Глава 5. ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ КРАЕВЫХ ЗАДАЧОпределение. Разностная схема аппроксимирует исходную дифференциальную задачу с p-мпорядком аппроксимации, если ψi = O(hp ) в точках xi , или в целом на сетке, если ||ψi ||h = O(hp ).Будем поступать так же, как и в случае выяснения порядка аппроксимации для задачи Коши вразностной задаче аппроксимации.Можно показать, что в целом на сетке схема, построенная интегро-интерполяционным методом,будет иметь второй порядок аппроксимации.
Вычисления громоздки и мы их опускаем, заметив, чтоневязка имеет порядок не хуже O(h2 ) даже в крайних узлах:ψ0 = O(h2 ),ψN = 0.Если мы требуем такой порядок аппроксимации, то можно сэкономить на вычислении параметров,вычисляя ai , di , ϕi по формуле прямоугольников, при этом получим следующие значения параметров: ai = k xi− 12 ;di = q(xi );ϕi = f (xi ).Если считать параметры по квадратурной формуле трапеций, то получим следующие выражения:111+;=ak(x)k(xii−1i)1di =q xi− 12 + q xi+ 12 ;2 ϕi = 1 f x 1 + f x 1 .i− 2i+ 22Рассмотрим вопрос о сходимости приближенного решения к точному.
Как обычно, обозначимzi = yi − ui — погрешность и напомним несколько определений.h→0Определение. Приближенное решение yi сходится к точному в точке xi , если zi −→ 0.h→0Определение. Приближенное решение yi сходится к точному на всей сетке, если ||zi ||h −→ 0.Определение. Если величина погрешности zi в каждой точке (или на всей сетке) есть O(hp ), тометод имеет p-й порядок точности.Можно установить, что сеточная норма ||z||C(ωi ) = O(h2 ).
Доказательство этого утверждения можно посмотреть в книге [1].Подставим в расчетную схему yi = ui + zi :(azx )x,i − di zi = −(aux )x,i + di ui − ϕi .Аналогичную операцию проведем для граничных условий:(azx )x,i − di zi = −ψi ;−a1 zx,0 + βz0 = −ψ;zN = 0.Как нетрудно заметить, задача для погрешности имеет ту же структуру, что и исходная разностнаясхема, с заменой правой части на невязку.После преобразований системы, которые мы снова опускаем, можно получить, что||z||h 6 M1 ||ψ||h .5.9. ЯВНАЯ РАЗНОСТНАЯ СХЕМА ДЛЯ УРАВНЕНИЯ ТЕПЛОПРОВОДНОСТИ97Разностная схема имеет 2-й порядок аппроксимации (||ψ||h = O(h2 )), а, следовательно, и второйпорядок точности ||z||h = O(h2 ). Подробные указания на то, как это получить, можно найти в [1].Проделав те же самые действия, можем получить оценку на приближенное решение в сеточнойнорме (по аналогии — уравнения очень похожи).||y||h 6 M1 (||ϕ||h + |µ1 | + |µ2 |),(5.29)— такой оценки и следовало ожидать.Теперь мы можем дать несколько определений.Определение.