Методы вычислений
 Скачать 2.39 Mb.
|
|
τ 2 Λ xx )(E − τ 2 Λ yy ) u n+1 j,m − u n j,m τ = Λu n j,m + f j,m , (x j , y m ) ∈ ω h , n = 0, 1 . . . , u n j,m = µ(x j , y m ), (x j , y m , t n ) ∈ Γ h , u 0 j,m = u 0 (x j , y m ), (x j , y m ) ∈ ¯ ω h . (6.54) СПФ имеет несколько реализаций в дробных шагах, но все эти реализации используют прогонки с переменой направления с продольного на поперечное, поэтому итерационный метод (6.54) будем называть методом переменных направлений. Покажем, что этот итерационный метод сходится, причем намного быстрее, чем рассмотренный выше явный метод простой итерации Теорема 6.1. Метод переменных направлений (6.54) сходится, при этом для квадратной области ¯ Ω (l x = l y = l) и квадратной сетки h y = h = l/N ) τ опт = h 2 sin πh l . (6.55) Д ока за тел ь ст во. Пусть z n j,m = u n j,m − u j,m — погрешность на й итерации в узле (x j , y m ). Для нее получаем следующую задачу − τ 2 Λ xx ´ ³ E − τ 2 Λ yy ´ z n+1 j,m − z n j,m τ = Λz n j,m , (x j , y m ) ∈ ω h , n = 0, 1 . . . , z n j,m = 0, (x j , y m , t n ) ∈ Γ h , z 0 j,m = u 0 (x j , y m ) − u j,m , (x j , y m ) ∈ Поскольку эта задача однородная, то ее решение можно представить в виде конечного ряда Фурье (Коэффициенты этого ряда найдем, подставив решение (6.57) враз- ностное уравнение задачи (6.56): T n+1 (k,l) − T n (k,l) τ ³ 1 + τ 2 λ (x) k ´ ³ 1 + τ 2 λ (y) l ´ = − ³ λ (x) k + где и λ (y) l — собственные значения операторов A x = и A y = соответственно. Отсюда получаем, что где − τ 2 λ (x) k ´ ³ 1 − τ 2 λ (y) l ´ ³ 1 + τ 2 λ (x) k ´ ³ 1 + τ 2 λ (y) l ´ Поэтому решение (6.57) можно записать в виде (где T 0 (k,l) — коэффициенты Фурье функции начальной ошибки z 0 102 Тогда для относительной ошибки на ой итерации будем иметь такую же неулучшаемую оценку (5.45), как в теореме поэтому о сходимости и скорости сходимости итераций можно судить по поведению правой части этого неравенства. Из свойств собственных значений операторов и следует, что max k,l ¯ ¯q (k,l) ¯ ¯ < 1, поэтому итерационный метод переменных направлений сходится. Найдем значение итерационного параметра τ , при котором в случае квадратной области (l x = l y = l) и квадратной сетки (h x = h y = итерации будут сходиться наиболее быстро. Для квадратной сетки собственные значения и совпадают, поэтому для них можно использовать общее обозначение λ k , с учетом которого выражение (примет вид − τ 2 λ k 1 + τ 2 λ k · 1 − τ 2 λ l 1 +Пусть = min k λ k = 4 h 2 sin 2 πh 2l , β = На квадрате [α, β] × [α, β] рассмотрим функцию, λ, µ) = 1 − τ 2 λ 1 + τ 2 λ · 1 − τ 2 µ 1 +Тогда q (k,l) = q(τ, λ k , λ l ). Легко убедиться, что функция |q(τ, λ, см. ее график на риса) будет принимать максимальное значение водном из углов квадрата [α, β] × [α, β], те При любом τ > 0 среднее из чисел меньше одного из крайних, поэтому max α≤λ,µ≤β |q(τ, λ, µ)| = max ³ 1 − τ 2 α ´ 2 ³ 1 + τ 2 α ´ 2 , ³ 1 − τ 2 β ´ 2 ³ 1 + τ 2 β ´ 2 = = max ³ ¯ ¯ ¯1 − τ 2 α ¯ ¯ ¯ 1 + τ 2 α , ¯ ¯ ¯1 − τ 2 β ¯ ¯ ¯ 1 + τ 2 β ´ 2 ≡ Q 2 (τ График функции Q(τ ) показан на рис. 1, б 0.5 1.0 0 100 200 300 0 100 200 300 |q| λ µ 0 0.1 0.2 0.50 0.75 1.00 Q τ а б Рис. 1. Графики функций а — |q(τ, λ, µ)| при τ = h; б — Q(τ ); l x = 1; l y = 1; N x = 10; N y = Нетрудно проверить, что ) = 1 − τ 2 α 1 +при < τ ≤ 2 √ αβ , τ 2 β − 1 1 +при поэтому функция Q(τ ) принимает при τ = 2/ √ αβ минимальное значение Итак = 1 − tg πh 2l 1 + tg πh 2l 2 ; (6.61) τ опт = h 2 sin πh l , (6.62) что и доказывает теорему. Чтобы при оптимальном значении итерационного параметра τ относительная погрешность встала меньше заданного положительного числа ε, необходимо и достаточно выполнение следующего неравенства те При малых h/l отсюда получаем ε 2 ln ¡ 1 − 2 πh 2l ¢ ≈ ln ε −2 πh l = 1 2π N Видим, что количество итераций в методе переменных направлений пропорционально числу N узлов, расположенных на стороне квадрата, в то время как в методе простой итерации для достижения точности требовалось проделать число итераций, пропорциональное N 2 . Таким образом, метод переменных направлений сходится гораздо быстрее, чем метод простой итерации. Например, при N = 100 получим, что при использовании метода переменных направлений будет достигнута относительная погрешность, не превосходящая ε = 0.45 · 10 −4 ≈ для произвольного начального приближения, если проделать n min ≈ 160 итераций. Ранее было показано, что метод простой итерации требует для этого около 000 итераций ЗАДАЧИ. Докажите, что операторы A x = и A y = попарно перестановочны. 6.2. Докажите корректность и устойчивость метода прогонки, используемого для решения уравнений СПН (6.30), (6.31). 6.3. Докажите, что СПН (6.34)—(6.36) для трехмерной задачи (экономична, имеет порядок аппроксимации O(τ + h 2 x + h 2 y + h 2 z ), не обладает свойством полной аппроксимации и необходимым условием ее устойчивости является выполнение неравенств (6.37). 6.4. Далее мы будем строить экономичные схемы на основе полностью неявной схемы u n τ = Λu n+1 + Рассмотрим СПФ (E − τ Λ xx ) (E − τ Λ yy ) u n+1 − u n τ = Λu n + которая аппроксимирует уравнение u xx + u yy + f (x, y, и реализуется методом дробных шагов − τ Λ xx ) ξ n+1/2 j,m = Λu n j,m + f n+1 j,m ; (6.69) (E − τ Λ yy ) ξ n+1 j,m = ξ n+1/2 j,m ; (6.70) u n+1 j,m = u n j,m + τ Исследовать свойства схемы (6.67) и выписать ее обобщение для случая уравнения теплопроводности стремя пространственными переменными u xx + u yy + u zz + f (x, y, z, t). (6.72) 6.5. Используя замену u n τ , (6.73) 106 записать метод дробных шагов (6.69)—(6.71) в виде схемы переменных направлений u n τ = Λ xx u n+1/2 + Λ yy u n + f n+1 ; (6.74) u n+1 − u n τ = Λ xx u n+1/2 + Λ yy u n+1 + аналогичной СПН (6.30), (6.31). Исследовать свойства схемы (6.74), (6.75) и выписать ее обобщение для уравнения (6.72). 6.6. Используя замену (6.73), записать метод дробных шагов) в виде схемы стабилизирующей поправки u n τ = Λ xx u n+1/2 + Λ yy u n + f n+1 ; (6.76) u n+1 − u n+1/2 τ = Λ yy ¡ u n+1 − аналогичной ССП (6.47), (6.48). Исследовать свойства схемы (6.76), (6.77) и выписать ее обобщение для уравнения (6.72). 6.7. Для решения уравнения (6.68) предлагается схема расщепления u n τ = Λ xx u n+1/2 + f n+1 ; (6.78) u n+1 − u n+1/2 τ = Показать, что она абсолютно устойчива, экономична, аппроксимирует уравнение (6.68) с порядком O(τ + h 2 x + h 2 y ) и не обладает свойством полной аппроксимации. Обладает ли она этим свойством при f ≡ 0? 6.8. Будем использовать схему расщепления (6.78), (6.79) на двух первых дробных шагах схемы предиктор – корректор u n τ /2 = Λ xx u n+1/4 + f n+1/2 ; (6.80) u n+1/2 − u n+1/4 τ /2 = Λ yy u n+1/2 ; (6.81) u n+1 − u n τ = Λu n+1/2 + Показать, что она аппроксимирует уравнение (6.68) с порядком O(τ 2 + h 2 x + h 2 y ), экономична, абсолютно устойчива, ноне обладает свойством полной аппроксимации. Обладает ли она этим свойством при f ≡ 0? 107 6.9. Пусть для решения разностной задачи Дирихле (5.36) используется итерационный процесс, основанный на схеме приближенной факторизации, y m ), (x j , y m , t n ) ∈ Γ h , u 0 j,m = u 0 (x j , y m ), (x j , y m ) ∈ где u 0 — начальное итерационное приближение, удовлетворяющее граничному условию (5.38). Докажите аналог теоремы 6.1 для метода переменных направлений (6.83). Какой из итерационных методов сходится быстрее — (6.54) или (6.83)? § 7. Метод адаптивных сеток. До сих пор конечно-разностные схемы для многомерных задач мы рассматривали лишь на равномерных сетках. Схемы на неравномерных сетках рассматривались в § 2.7 для стационарного уравнения теплопроводности с одной пространственной переменной. При этом было показано, что применение неравномерных сеток может привести к повышению точности приближенного решения по сравнению с решением, полученным на равномерной сетке. В настоящем параграфе мы на примере задачи Дирихле (5.3) для уравнения Пуассона+ u yy + f (x, y) = 0, (x, y) ∈ Ω, u(x, y) = µ(x, y), (x, y) ∈ γ = покажем, как строится разностная схема на неравномерной сетке и укажем один из методов построения адаптивных сеток. Адаптивными будем называть сетки, при построении которых используется информация о поведении решения дифференциальной задачи. Например, в области быстрого изменения решения сетку можно сгустить, а там, где решение меняется медленно, взять разреженную сетку. Это может привести к более точному представлению решения и обеспечить для заданного числа узлов большую точность, чем на равномерной сетке. Кроме того, адаптивные сетки подстраиваются к криволинейным границам области решения Ω. 108 Как пояснялось в § 2.7, для того чтобы воспользоваться методом адаптивных сеток, необходимо предварительно задать управляющую функцию w(x, y) ≥ 1, в которой каким-то образом должна быть учтена информация о поведении погрешности численного решения. В общем случае управляющую функцию задают зависящей от производных до некоторого порядка, причем функция w должна принимать большие значения в подобластях больших значений производных отрешения. Например, для одномерных задач управляющая функция может быть задана формулой (2.7.54). Далее будем считать, что управляющая функция w(x, y) задана и может быть вычислена в любой точке, y) ∈ Метод адаптивных сеток позволяет получать численное решение задачи для области Ω с криволинейной границей γ сложной формы, однако ради простоты изложения, будем предполагать, что область Ω является, как ив предыдущих параграфах, прямоугольником в плоскости, те. Пусть область Ω требуется покрыть неравномерной сеткой с количеством узлов (причем на левой и правой сторонах прямоугольника Ω должно размещаться одинаковое количество (N 2 + 1) узлов, а на нижней и верхней (N 1 + 1) узлов. Наряду с физической областью Ω будет использоваться вычислительная область Ξ = (0, 1) × (0, 1) — единичный квадрат в плоскости. Пусть = x(ξ), (7.2) — произвольное достаточно гладкое взаимно-однозначное невырожденное (для определенности с положительным якобианом) отображение замкнутой области ¯ Ξ на ¯ Ω, такое, что левая, нижняя, правая и верхняя стороны прямоугольника Ω являются образами левой, нижней, правой и верхней сторон квадрата Ξ соответственно. Здесь x = (x, y), ξ = (ξ, η), x(ξ) = (x(ξ, η), y(ξ, η)). Покроем ¯ Ξ равномерной прямоугольной сеткой с количеством узлов (N 1 + 1) и шагом в направлении оси Oξ и с количеством узлов (N 2 + 1) и шагом в направлении оси Oη. Тогда 1/N 1 , h 2 = 1/N 2 , а координаты узлов в вычислительной области задаются формулами j 1 h 1 , η j 2 = j 2 h 2 , j 1 = 0, . . . , N 1 , j 2 = 0, . . . , Для узлов в вычислительной области будет использоваться обозначение, где j = (j 1 , j 2 ) — мультииндекс. Совокупность всех узлов обозначим ¯ Ξ h , множество граничных узлов — ∂ ¯ Ξ h . Множество ¯ Ξ h \∂ представляет собой совокупность внутренних узлов сетки. Криволинейной сеткой назовем совокупность узлов x j = (x j , являющихся образами узлов при отображении (7.2). Таким образом, координаты узлов неравномерной сетки вычисляются по формулам x(ξ j 1 , η j 2 ), y j 1 ,j 2 = y(ξ j 1 , η j 2 ), j 1 = 0, . . . , N 1 , j 2 = 0, . . . , с помощью отображения (7.2). Отметим, что координаты узлов равномерной на Ω сетки однозначно определяются заданием числа узлов и N 2 . Для неравномерной сетки это не так. Для заданных и существует бесконечное множество неравномерных сеток ¯ Ω h , а отдельный представитель этого семейства будет определяться заданием конкретного отображения (7.2). Отображение (7.2) заранее неизвестно и наша задача — определить его. Для этого будем использовать метод эквирас- пределения [16]. 7.2. Построение сетки в многомерной области начинается с построения сетки на ее границе. Для расстановки узлов на отрезках прямых, из которых составлена граница γ, будем использовать метод эквирас- пределения [16], при этом сетка будет сгущаться или разрежаться в соответствии с заданной управляющей функцией |