Численные методы. Соснин (2005) (1160462), страница 20

Просмтор этого файла доступен только зарегистрированным пользователям. Но у нас супер быстрая регистрация: достаточно только электронной почты!

Текст из файла (страница 20)

Тогданевязка посчитается так:n+1n+1ψjn+1 = (−un+1) + ϕn+1− fjn+1 + O(τ + h2 ).t,j + uxx,j + fjjЕсли функцию ϕn+1взять равной fjn+1 с точностью O(τ + h2 ), то выражение для невязки сильноjсократится и примет вид:ψjn+1 = O(τ + h2 ).Как видно, эта разностная схема имеет первый порядок аппроксимации по τ и второй по h.Для исследования на устойчивость воспользуемся методом гармоник. Сопоставим нашему разностному уравнению однородное уравнение:n+1n+1yjn+1 − yjnyj+1− 2yjn+1 + yj−1=.τh2Проделаем те же действия, что и в предыдущем случае. Подставив в качестве решения yjn = q n eijhϕи сократив множители, получим такое уравнение относительно параметров этого решения:eihϕ − 2 + e−ihϕq−1=qτh24q−1hϕ= −q 2 sin2.τh2=⇒Выразим отсюда q :4τhϕq 1 + 2 sin2=1h2=⇒q=14τ1 + 2 sin2hhϕ2.Очевидно (так как знаменатель всегда больше или равен единице), что эта неявная разностнаясхема абсолютно устойчива (устойчива при любых значениях τ и h ).Посмотрим, как можно получить решение разностного уравнения из системы (5.33).

Перепишемего так:τ n+12ττ n+1y−1+yjn+1 + 2 yj+1= −yjn + τ ϕn+1.(5.35)j−1j22hhhОбозначим для удобстваAj =τ,h2Bj =τ,h2Cj =1+2τh2,Fjn = yjn − τ ϕn+1.jТогда (5.35) будет выглядеть постройнее:n+1n+1Aj yj−1− Cj yjn+1 + Bj yj+1= −Fjn .Рассмотрим эту систему при n = 0 :11Aj yj−1− Cj yj1 + Bj yj+1= −Fj0 ,j = 1, N − 11— это система линейных алгебраических уравнений относительно y 1 = (y01 , y11 , . . . , yN). Перепишем еев компактном виде:M y1 = F 0 ,104Глава 5. ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ КРАЕВЫХ ЗАДАЧгде матрица F 0 состоит из элементов Fj0 (которые полностью определяются начальными условиями),а M — из коэффициентов Aj , Bj , Cj , расположенных на трех диагоналях (то есть матрица имеет трехдиагональный вид).

Следовательно, применим метод прогонки, и мы можем найти сеточнуюфункцию на первом временном слое.Поступая так дальше, мы сможем определить искомую сеточную функцию на всех временных слоях(«послойно» применяя метод прогонки).5.11Разностная схема с весами для уравнения теплопроводностиРассмотрим теперь не минимально возможный шаблон, а «избыточный». Будем аппроксимироватьпроизводные в шести узлах. «Избыточность» схемы скомпенсируем введением некоторого параметра —весового множителя.Соответствующая этому шаблону разностная схема такова: n+1yj − yjnτyj0y0n+1n+1yN= σn+1n+1nnyj+1− 2yjn+1 + yj−1yj+1− 2yjn + yj−1+(1−σ)+ ϕnj ,h2h20j = 1, N − 1,n = 0, K − 1;(5.36)j = 0, N ;= u (xj ),= µ1 (tn+1 ),n = 0, K − 1;= µ2 (tn+1 ),n = 0, K − 1.— это так называемая разностная схема с весами для уравнения теплопроводности.Список вопросов остается тем же.Рассмотрим вопрос об аппроксимации.

Выясним, как ведет себя невязка:n+ 1 defψj 2 =un+1− unjjnn−+ σun+1xx,j + (1 − σ)uxx,j + ϕj .τКак мы уже поступали, разложим функции unj и un+1в ряд Тейлора:jn+ 12un+1j= ujunj= ujn+ 12τ1 n+ 1 τ 2+ utt,j 2+ O(τ 3 );2 221 n+ 1 τ 2n+ 1 τ− ut,j 2 + utt,j 2+ O(τ 3 ).2 22n+ 12+ ut,jТогда получим:n+ 12ψjn+ 12= −ut,jnn+ O(τ 2 ) + σun+1xx,j + (1 − σ)uxx,j + ϕj .(5.37)5.11. РАЗНОСТНАЯ СХЕМА С ВЕСАМИ ДЛЯ УРАВНЕНИЯ ТЕПЛОПРОВОДНОСТИ105Теперь в представлении второй разностной производной разложим все вхождения функции в рядТейлора с членами до пятого порядка включительно:1(u(xj+1 , t) − 2u(xj , t) + u(xj−1 , t)) =h2h2h3h4h5= {u(xj±1 , t) = u(xj , t) ± ux (xj , t)h + uxx (xj , t) ± uxxx (xj , t) + uxxxx (xj , t) ± uxxxxx (xj , t) + O(h6 )} =264!5!h24= uxx (xj , t) + uxxxx (xj , t) + O(h ).12uxx (xj , t) =Воспользуемся этим разложением для слагаемых в выражении (5.37):h2n+ 1+ O(h4 ) = {Разложение в ряд Тейлора с центром в точке xj 2 } =12τh2h2 τn+ 21n+ 12n+ 21n+ 12= uxx,j + uxxt,j · + uxxxx,j··+ uxxxxt,j· + O(τ 2 + h4 );21212 2h2n+ 1= unxx,j + unxxxx,j ·+ O(h4 ) = {Разложение в ряд Тейлора с центром в точке xj 2 } =12τh2h2 τn+ 21n+ 21n+ 12n+ 12= uxx,j − uxxt,j · + uxxxx,j− uxxxxt,j· + O(τ 2 + h4 ).··21212 2n+1n+1= un+1uxx,jxx,j + uxxxx,j ·nuxx,jТаким образом, выражение для невязки принимает вид:n+ 21ψjn+ 21= −ut,jn+ 1n+ 12Добавляя и вычитая fjn+ 12ψjn+ 1n+ 12+ uxx,j2 + uxxxx,jn+ 1n+ 21= −ut,j 2 + uxx,j2 + fjh21 n+ 121 n+ 12 h2+ ϕnj + (σ − )uxxt,jτ + O(τ 2 + h4 ).τ + (σ − )uxxxxt,j122212, получим эквивалентное выражение:n+ 12+ uxxxx,jh21 n+ 121 n+ 12 h2n+ 1+ ϕnj − fj 2 + (σ − )uxxt,jτ + O(τ 2 + h4 ).τ + (σ − )uxxxxt,j122212Согласно уравнению теплопроводности, первые три слагаемых обращаются в нуль:n+ 12ψjn+ 12= uxxxx,jh21 n+ 211 n+ 21 h2n+ 1+ ϕnj − fj 2 + (σ − )uxxt,jτ + O(τ 2 + h4 ).τ + (σ − )uxxxxt,j122212При σ = 21 схема (5.36) называется симметричной.

Тогда в последнем равенстве последние слагаемые обнулятся, и с помощью условия на параметр ϕnjn+ 12ϕnj = fj+ O(τ 2 + h2 )мы можем достичь такого порядка аппроксимации:n+ 12ψj= O(τ 2 + h2 ).Теперь вернемся на шаг назад и воспользуемся тем, чтоut = uxx + f =⇒ uxxt = uxxxx + fxx .Тогда формула для невязки будет несколько иной:21h2 n+ 121 n+ 12h2n+ 1n+ 1n+ 1 hψj 2 = ϕnj − fj 2 − fxx,j2 ·+ (σ − )τ +uxxt,j + (σ − )uxxxxt,j· τ + O(τ 2 + h4 ).122122121 h2Взяв σ = −, мы обнулим четвертое слагаемое, а коэффициент при пятом оценится как O(h4 ).2 12τОсталось потребовать, чтобыn+ 21ϕnj = fjn+ 1+ fxx,j2 ·h2+ O(τ 2 + h4 ),12106Глава 5. ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ КРАЕВЫХ ЗАДАЧтогда порядок аппроксимации будет таков:n+ 12ψj= O(τ 2 + h4 ).При данном σ выражение (5.36) называется разностной схемой повышенного порядка аппроксимации. При всех остальных σ порядок аппроксимации будет меньше:n+ 12ψj= O(τ + h2 )при условии, чтоn+ 21ϕnj = fj+ O(τ + h2 ).Теперь исследуем схему на устойчивость методом гармоник.

Для начала запишем однородноеуравнение:n+1n+1nnyjn+1 − yjnyj+1− 2yjn+1 + yj−1yj+1− 2yjn + yj−1=σ+(1−σ).τh2h2Подставим в качестве решения yjn = q n eijhϕ . Сокращая степени q, получим:q−1q1= σ 2 (eihϕ − 2 + e−ihϕ ) + (1 − σ) 2 (eihϕ − 2 + e−ihϕ ) ⇐⇒τhhq−1q1hϕhϕ= −4σ 2 sin2− 4(1 − σ) 2 sin2.τh2h2Отсюда получаем выражение для q :ττhϕhϕsin2) = 1 − (1 − σ)4 2 sin2=⇒h22h2τhϕ1 − (1 − σ)4 2 sin2h2 .q=τ2 hϕ1 + 4σ 2 sinh2q(1 + 4σДля получения условий на устойчивость мы требуем, чтобы |q| 6 1. В данном случае это дает дванеравенства:ττhϕ2 hϕ6 1 + 4σ 2 sin2; 1 − (1 − σ)4 2 sinh2h2hϕτhϕτ6 1 − (1 − σ)4 2 sin2.−1 − 4σ 2 sin2h2h2Первое выполнено всегда, так как τ > 0. Второе перепишется так:−8στhϕτhϕsin26 −4 2 sin2+2h22h2⇐⇒σ>1−2h24τ sin2hϕ2— оно должно быть выполнено при любом ϕ.

Взяв максимум по правой части, приходим к окончательному условию для σ :1 h2.σ> −2 4τ1Значение σ =удовлетворяет этому неравенству. Это означает, что соответствующий метод яв2ляется абсолютно устойчивым.Последним нашим долгом будет обосновать возможность вычисления приближения по этой схеме.Для этого перепишем (5.36) в таком виде:στ n+1ττ n+1ny− (1 + 2σ 2 )yjn+1 + σ 2 yj−1= −yjn − (1 − σ)τ yxx,j− ϕnj .h2 j+1hh5.12. РАЗНОСТНЫЕ СХЕМЫ ДЛЯ УРАВНЕНИЯ ТЕПЛОПРОВОДНОСТИ ОСОБОГО ТИПА 107Теперь почти очевидно, что мы можем получить все yjk (j = 0, N , k = 0, K). Действительно,фиксируем n = 0.

Тогда мы получим систему линейных уравнений с трехдиагональной матрицей.Правые части уравнений мы можем найти, используя начальные условия (заметим, что ϕnj мы задалипри исследовании порядка аппроксимации). После этого применяется метод прогонки, после которогостановятся известны все yj1 . Теперь можно увеличить n и снова получить СЛАУ, подставив в правуючасть уравнений только что найденные yj1 . Так действуем, пока не найдем все yjk .Мы рассмотрели схему для нахождения численного решения простейшей краевой задачи. Для неесуществует более простое аналитическое решение, но в общем случае его может и не быть. В то жевремя вся методика построения разностных схем и их решение достаточно легко переносятся на болеесложные задачи, к которым мы и перейдем.5.12Разностные схемы для уравнения теплопроводности особого типаРазностная схема для уравнения теплопроводности с переменными коэффициентамиРассмотрим такую краевую задачу:∂u∂∂uρ(x,t)=k(x,t)+ f (x, t),∂t∂x∂xu(0, t) = µ1 (t), 0 6 t 6 T ;u(1, t) = µ2 (t), 0 6 t 6 T ;u(x, 0) = u0 (x), 0 6 x 6 1.0 < x < 1, 0 < t 6 T ;Аналитического выражения для решения нет.

Тем не менее, известно, что если всюду верны неравенства0 < c1 6 ρ(x, t);0 < c2 6 k(x, t),то оно существует и единственно.Теперь будем приближать ut соответствующей разностной производной:∂u≈ unt,j ,∂tа производные по x с использованием интегро-интерполяционного метода можно представить следующим образом:∂∂u1uj+1 − ujuj − uj−1k(x, t)≈ aux x,j =aj+1− aj,∂x∂xhhhгде aj вычисляются по такой формуле:1aj = hZxj−1dx k(x, t)≈ k(xj− 21 , t).xj−1Используя шаблон из шести точек, мы построим разностную схему. Выкладки аналогичны предыдущим схемам, поэтому опустим их и приведем только окончательный результат:n+1nyjn+1 − yjn= σ ayx x,j + (1 − σ) ayx x,j + ϕnj ; j = 1, N − 1, n = 0, K − 1ρ(x,t)jτyj0 = u0 (xj ), j = 0, N ;(5.38)ny0 = µ1 (tn ), n = 1, K;nyN= µ2 (tn ), n = 1, K.108Глава 5. ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ КРАЕВЫХ ЗАДАЧЗдесь остаются нефиксированными момент времени t в первом уравнении (от него зависят aj и1ρ) и параметр метода σ.

Характеристики

Список файлов лекций