1612725465-eb5831e2a489b993bf8c6c33f84310db (828847), страница 21
Текст из файла (страница 21)
, Ny − 1 из уравнения(D − ωL)nnun+1un+1j−1,m − 2ūj,m + uj+1,mj,m−1 − 2ūj,m + uj,m+1+= −fj,mh2xh2yв цикле по j = 1, . . . , Nx − 1 определяется вспомогательная величинаūj,m , которая сразу подправляется по формуле релаксации (12.6), поэтому для реализации метода требуется только один двумерный массив.Оптимальное значение итерационного параметра τ в случае квадратной сетки (hx = hy = h, Nx = Ny = N ), покрывающей квадрат(lx = ly = l), вычисляется по формуле2τопт =1 + sinπhl.(12.8)Остальные три метода: СПФ — схема приближенной факторизации, СПН — схема переменных направлений и ССП — схемастабилизирующей поправки — подробно описаны в § 6.
Для программной реализации этих методов достаточно двух двумерных массивов.12.2. Тестовые задачи. Задача (5.3) рассматривается в единичном(lx = ly = 1) квадрате Ω, при этом функция f задана в области Ω, афункция µ — на ее границе γ = ∂Ω.135Задача 1:f (x, y) = 4,µ(0, y) = µ(1, y) = y(1 − y),µ(x, 0) = µ(x, 1) = x(1 − x).Задача имеет точное решение (рис. 6, а):u(x, y) = x(1 − x) + y(1 − y),(12.9)которое является и решением разностной задачи (5.36). Таким образом,на этом решении разностная схема (5.36) аппроксимирует дифференциальную задачу (5.3) с бесконечным порядком.Задача 2:f (x, y) = −4 + 12(x + y) − 12(x2 + y 2 ),µ(0, y) = µ(1, y) = y 2 (1 − y)2 ,µ(x, 0) = µ(x, 1) = x2 (1 − x)2 .Задача имеет точное решение (рис. 6, б):u(x, y) = x2 (1 − x)2 + y 2 (1 − y)2 .Задача 3:(12.10)f (x, y) = 8π 2 sin(2πx) sin(2πy),µ(0, y) = µ(1, y) = 0,µ(x, 0) = µ(x, 1) = 0.Задача имеет точное решение (рис.
7, а):u(x, y) = sin(2πx) sin(2πy).Задача 4:(12.11)f (x, y) = 32π 2 sin(4πx) sin(4πy),µ(0, y) = µ(1, y) = 0,µ(x, 0) = µ(x, 1) = 0.Задача имеет точное решение (рис. 7, б):u(x, y) = sin(4πx) sin(4πy).(12.12)Для следующих задач точное решение не известно, поэтому они решались численно. Графики полученных решений изображены нарис. 8, 9, а.136uu0.60.10.40.050.2000.750.25x0.50.50.250.75y0000.750.25аx0.50.50.250.75yб0Рис. 6.
Графики решений: а — задача 1; б — задача 2uu1100-10-100.750.25x0.50.50.250.750.750.25y0xа0.50.50.250.75y0бРис. 7. Графики решений: а — задача 3; б — задача 4uu110.50.500000.750.25x0.50.50.250.7500.750.25yxа0.50.50.250.750Рис. 8. Графики решений: а — задача 5; б — задача 6137yбЗадача 5:f (x, y) = sin2 (πxy),µ(0, y) = µ(1, y) = sin(πy),µ(x, 0) = µ(x, 1) = x(1 − x).Задача 6:µ¶x+1f (x, y) = arctan,y+1µ(0, y) = µ(1, y) = 0,µ(x, 0) = sin2 (πx), µ(x, 1) = cosh [x(1 − x)] − 1.Задача 7:f (x, y) = |2x − y|,µ(0, y) = µ(1, y) = y(1 − y),µ(x, 0) = | sin(2πx)|, µ(x, 1) = | sin(2πx)|e2x .10-1u5610-24zn3210-31000.750.25x0.50.50.250.750510-4y100а23416200n300400500бРис. 9.
а — график решения задачи 7; б — поведение нормы разностиz n = kun+1 − un kC ; метод ПВР для задачи 7 при различных значенияхпараметра верхней релаксации: τ = 1.5 (1), τ = 1.6 (2), τ = 1.7 (3),τ = 1.8 (4), τ = 1.9 (5), τ = 1.95 (6). Размер сетки Nx = Ny = 5012.3. Задания к практическим занятиям. Для каждого из двухуказанных в задании итерационных методов определить экспериментально оптимальное значение τопт итерационного параметра τ , при котором итерации сходятся наиболее быстро.
В расчетах использоватьквадратную сетку с шагом h. Выяснить характер зависимости τопт от h.Экспериментально проверить сходимость разностной схемы на последо138вательности измельчающихся сеток. На основе численных экспериментов выявить, какой из двух итерационных методов предпочтительнеедля рассматриваемого класса тестовых задач.В процессе решения задачи выводить в файл данные, характеризующие сходимость итерационного процесса, например, номер итерации nи норму разности двух последовательных итераций z n = kun+1 − un kC .Эту же информацию можно визуализировать (см. рис.
9, б). Если известно точное решение u(x, y) тестовой задачи, то дополнительно в файлвыводить и норму kun − (u)h kC отклонений итераций от точного решения. Графическая информация должна включать графики численногорешения на каждой итерации в одном или нескольких заданных сечениях x = const или y = const.Задание 1. Метод простой итерации и СПФ.Задание 2. Метод Зейделя и СПН.Задание 3. Метод ПВР и ССП.Задание 4. Метод простой итерации и СПН.Задание 5. Метод Зейделя и ССП.Задание 6.
Метод ПВР и СПФ.Задание 7. Метод простой итерации и ССП.Задание 8. Метод Зейделя и СПФ.Задание 9. Метод ПВР и СПН.Задание 10. Метод простой итерации и метод Зейделя.Задание 11. Метод Зейделя и метод ПВР.Задание 12. Метод простой итерации и метод ПВР.12.4.
Следующая группа заданий предусматривает решение нелинейных уравнений эллиптического типа. Таковыми являются, например, уравнения метода эквираспределения (7.12). В заданиях требуетсяпостроить адаптивную сетку с помощью двух итерационных методови провести анализ скорости сходимости сравниваемых методов. Управляющую функцию вычислять по формуле (7.22), где u — одно из точныхрешений (12.9)—(12.12) задач 1—4. Выяснить влияние параметров a и βна вид адаптивной сетки.Таким образом, следующие 12 заданий, с 13-го по 24-е, заключаютсяв сравнении итерационных методов, перечисленных в первых 12 заданиях, но не на задаче Дирихле для уравнения Пуассона, а на нелинейнойзадаче (7.13) для уравнений метода эквираспределения.139Ответы, указания, решенияτ1= .2h6τττ11.2. ν 2 = ; ϕnj = f (xj , tn ) + ft (xj , tn ) + ν fxx (xj , tn ).h622√215h1.3.
σ = −;τ= √ .262ν 5Р е ш е н и е. Пусть u(x, t) — решение однородного уравнения теплопроводностиut = νuxx .(13.1)1.1. νРассмотрим на этом решении локальную погрешность аппроксимациив узле (xj , tn ) дифференциального уравнения (13.1) разностным уравнением (1.101):µ¶¯ττ2¯nψh,j= (Lh (u)h − fh ) ¯=u+u+u(xj , tn )−tttttt26(xj ,tn )µ¶h4h2uxxxxxx (xj , tn+1 )−−σν uxx + uxxxx +12360µ¶24hh−(1 − σ)ν uxx + uxxxx +uxxxxxx (xj , tn ) + O(τ 3 + h6 ) =12360³ττ2τ2= ut + utt + uttt − σν uxx + τ uxxt + uxxtt +262´224hhh+ uxxxx + τ uxxxxt + O(τ 2 h2 ) +uxxxxxx + O(τ h4 ) −1212360µ¶24hh−(1 − σ)ν uxx + uxxxx +uxxxxxx + O(τ 3 + h6 ).12360В последнем равенстве все производные вычисляются в одном и том жеузле (xj , tn ). Далее для сокращения записи формул этот узел указываться не будет.После приведения подобных получи쵶ττ2τh2nψh,j= ut − νuxx + utt + uttt − σντ uxxt + uxxtt + uxxxxt −26212µ 2¶hh4−νuxxxx +uxxxxxx + O(τ 3 + h6 + τ 2 h2 + τ h4 ).12360140Учитывая, что для решения уравнения (13.1) выполняются равенстваut − νuxx = 0,utt = νuxxt ,utxxxx = νuxxxxxx ,utxx = νuxxxx ,uttxx = ν 2 uxxxxxx ,utt = ν 2 uxxxx ,uttt = ν 3 uxxxxxx ,перепишем выражение для погрешности аппроксимации в следующемвиде:µ¶h2τ 2nν − σν 2 τ − νuxxxx +ψh,j=212¶µ 2h4h2τ2τ 3ν − τ σν 2− σν 3−νuxxxxxx + O(τ 3 + h6 + τ 2 h2 + τ h4 ).+6122360Пусть1h2−.(13.2)2 12τ νобратится в нуль и, следовательно,σ=Тогда коэффициент при uxxxxµ¶ν h4nψh,j=− τ 2 ν 2 uxxxxxx + O(τ 3 + h6 + τ 2 h2 + τ h4 ).12 20Если теперь приравнять нулю коэффициент при uxxxxxx , т.
е. взятьh2τ= √ ,ν 20то для погрешностиnψh,j= O(h6 ).аппроксимации(13.3)будемиметьвыражение1.4. O(τ 2 + h2 ).У к а з а н и е. Покажите, что разностное уравнение схемы (1.102)можно переписать в следующем виде:¢ un+1un+1+ unj− unj1 ¡ 2 n+1jjh Λuj − h2 Λunj += νΛ,6ττ2илиun+1− unjj= νΛτ·µ1h2−2 6τ ν¶µun+1+j1411h2+2 6τ ν¶¸unj .(13.4)Таким образом, схема (1.102) является схемой с весами, причемσ=1h2−.2 6τ ν1.5. O(τ 2 + h4 ).У к а з а н и е. Покажите, что схема (1.103) является схемой с весамиповышенного порядка аппроксимации.1.6.
У к а з а н и е. При использовании трехточечного шаблонадля аппроксимации производной ux (0, t) в левом граничном узле явнаясхема с порядком аппроксимации O(τ + h2 ) запишется так:un+1− unjunj−1 − 2unj + unj+1j=ν+ fjn ,τh2−3un0 + 4un1 − un2+ γun0 = µ0 (tn ),2hunN = µl (tn ), n = 0, . . . , M,u0j = u0 (xj ),j = 1, . . . , N − 1,(13.5)j = 0, . . . , N.Алгоритм получения решения на временно́м слое tn+1 включает в себя два шага:1) вычисление с помощью разностного уравнения значений un+1воjвнутренних узлах сетки xj , j = 1, .
. . , N − 1;2) определение значений un+1и un+1в граничных узлах x0 = 00Nи xN = l на основе разностных уравнений в этих узлах с использованием при необходимости значений un+1во внутренних узлах, ужеjвычисленных на первом шаге. Например:un+1=04un+1− un+1− 2hµ0 (tn+1 )12.3 − 2hγ(13.6)Условие γ ≤ 0 исходной задачи (1.104) гарантирует, что 3 − 2hγ 6= 0 прилюбом h.Теперь для аппроксимации той же производной ux (0, t) будем использовать двухточечный шаблон.
Чтобы получить на решении задачиаппроксимацию второго порядка по h, потребуем выполнения дифференциального уравнения вплоть до границы x = 0. Тогдаuxx (0, t) =1(ut (0, t) − f (0, t))ν142(13.7)и разностная схемаun+1− unjunj−1 − 2unj + unj+1j=ν+ fjn , j = 1, . . . , N − 1,2τhµ¶un1 − un0h un+1− un0n0−− f0 + γun0 = µ0 (tn ),h2ντunN = µl (tn ), n = 0, . . . , M,u0j = u0 (xj ),(13.8)j = 0, . . . , Nбудет аппроксимировать задачу (1.104) с порядком O(τ + h2 ).Алгоритм решения по этой схеме совпадает с алгоритмом реализации схемы (13.5), только вместо (13.6) теперь используем другую формулу:¢h2 ¡ nun+1+u0 + τ f0n+1 − hµ0 (tn+1 )12ντun+1=.(13.9)0h2− hγ1+2ντ1.7. У к а з а н и е.
Соответствующие полностью невные схемыимеют видn+1un+1+ un+1un+1− unjjj−1 − 2ujj+1=ν+ fjn+1 , j = 1, . . . , N − 1,2τh−3un+1+ 4un+1− un+1012+ γun+1= µ0 (tn+1 ),02hun+1= µl (tn+1 ), n = 0, . . . , M − 1,Nu0j = u0 (xj ),(13.10)j = 0, . . .
, N,иn+1un+1− unjun+1+ un+1jj−1 − 2ujj+1=ν+ fjn+1 , j = 1, . . . , N − 1,2τhµ¶un+1− un+1h un+1− un0n+1100−− f0+ γun+1= µ0 (tn+1 ),0h2ντun+1= µl (tn+1 ), n = 0, . . . , M − 1,Nu0j = u0 (xj ),(13.11)j = 0, . . . , N.В обеих схемах разностные уравнения во внутренних узлах являютсятрехточечными:n+1aj un+1+ bj un+1j−1 − cj ujj+1 = dj ,143j = 1, . . . , N − 1,aj = bj = r,dj = −unj − τ fjn+1 ,cj = 1 + aj + bj ,r=ντ,h2поэтому для вычисления значений un+1можно использовать метод проjгонкиun+1= µl (tn+1 ),Nun+1= ξj un+1jj+1 + ηj ,j = N − 1, .
. . , 0,прогоночные коэффициенты которого задаются следующими формулами:ξj =bj,cj − aj ξj−1ηj =aj ηj−1 − dj,cj − aj ξj−1j = 1, . . . , N − 1.Чтобы воспользоваться этими формулами, необходимо предварительноопределить начальные коэффициенты ξ0 и η0 . Для схемы (13.11) изграничного условия, записанного в виде (13.9), получаемun+1= ξ0 un+1+ η0 ,01(13.12)где1ξ0 =21+h− hγ2ντ,¢h2 ¡ nu0 + τ f0n+1 − hµ0 (tn+1 )η0 = 2ντ.h21+− hγ2ντДля схемы (13.10) граничное условие при x = 0 также можно записать в виде (13.12). Для этого надо исключить из формулы (13.6) значение un+1, использовав для этого разностное уравнение во внутреннем2узле с номером j = 1:a1 un+1− c1 un+1+ b1 un+1= d1 .012В результате получимξ0 =1.8. σ ≥ 1 −4b1 − c1,b1 (3 − 2hγ) − a1η0 = −d1 + 2hb1 µn+10.b1 (3 − 2hγ) − a1h2.2τ ν1.9.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
