Методы вычислений
 Скачать 2.39 Mb.
|
|
j,m−1 − 2u n+1 j,m + u n j,m+1 h 2 y + f j,m = 0, причем расчет можно проводить по строкам, те. для каждого, . . . , N y − 1 эта величина вычисляется при изменении = 1, . . . , N x − 1. Расчет начинается с узла (x 1 , y 1 ) и проводится по строкам сетки с перебором строк снизу вверх. Найденное значение u n+1 j,m размещается в массиве на месте u n j,m , поэтому для программной реализации метода Зейделя требуется только один двумерный массив. ПВР — метод последовательной верхней релаксации отличается от метода Зейделя тем, что вначале по схеме, аналогичной (вычисляется вспомогательная величина ¯ u j,m : −D¯ u j,m + Lu n+1 j,m + Ru n j,m + f j,m = а затем уже определяется u n+1 j,m : u n+1 j,m = ω¯ u j,m + (1 − где 1 ≤ ω < 2 — параметр верхней релаксации. Приведем формулы (12.5), (12.6) к общей схеме (12.1) итерационных методов. Из равенства (12.6) следует, что u n j,m ω + u n j,m . 134 Поэтому уравнение (12.5) можно переписать как u n j,m ω − Du n j,m + Lu n j,m + Ru n j,m + Lu n+1 j,m − Lu n j,m + f j,m = или u n j,m ω + Λu n j,m + L ¡ u n+1 j,m − u n j,m ¢ + f j,m = те+ Таким образом, метод ПВР — неявный итерационный метод с оператором и итерационным параметром τ , равным параметру релаксации ω. При ω = 1 метод ПВР совпадает с методом Зейделя (12.3). Расчетные формулы следуют из формул (12.5) и (12.6). При каждом фиксированном m = 1, . . . , N y − 1 из уравнения 2¯ u j,m + u n j+1,m h 2 x + u n+1 j,m−1 − 2¯ u j,m + u n j,m+1 h 2 y = в цикле по j = 1, . . . , N x − 1 определяется вспомогательная величина, которая сразу подправляется по формуле релаксации (12.6), поэтому для реализации метода требуется только один двумерный массив. Оптимальное значение итерационного параметра τ в случае квадратной сетки (h x = h y = h, N x = N y = N ), покрывающей квадрат l y = l), вычисляется по формуле τ опт = 2 1 + Остальные три метода СПФ — схема приближенной факторизации, СПН — схема переменных направлений и ССП — схема стабилизирующей поправки — подробно описаны в § 6. Для программной реализации этих методов достаточно двух двумерных массивов. Тестовые задачи. Задача (5.3) рассматривается в единичном l y = 1) квадрате Ω, при этом функция f задана в области Ω, а функция µ — на ее границе γ = ∂Ω. 135 Задача 1: f (x, y) = 4, µ(0, y) = µ(1, y) = y(1 − y), µ(x, 0) = µ(x, 1) = x(1 − Задача имеет точное решение риса которое является и решением разностной задачи (5.36). Таким образом, на этом решении разностная схема (5.36) аппроксимирует дифференциальную задачу (5.3) с бесконечным порядком. Задача 2: f (x, y) = −4 + 12(x + y) − 12(x 2 + y 2 ), µ(0, y) = µ(1, y) = y 2 (1 − y) 2 , µ(x, 0) = µ(x, 1) = x 2 (1 − Задача имеет точное решение (рис. 6, б, y) = x 2 (1 − x) 2 + y 2 (1 − Задача 3: f (x, y) = 8π 2 sin(2πx) sin(2πy), µ(0, y) = µ(1, y) = 0, µ(x, 0) = µ(x, 1) = Задача имеет точное решение риса) Задача 4: f (x, y) = 32π 2 sin(4πx) sin(4πy), µ(0, y) = µ(1, y) = 0, µ(x, 0) = µ(x, 1) = Задача имеет точное решение (рис. 7, б, y) = sin(4πx) Для следующих задач точное решение неизвестно, поэтому они решались численно. Графики полученных решений изображены на риса 0.75 y u а б Рис. 6. Графики решений а — задача 1; б — задача 2 -1 0 1 0 0.25 0.5 0.75 x 0 0.25 0.5 0.75 y u -1 0 1 0 0.25 0.5 0.75 x 0 0.25 0.5 0.75 y u а б Рис. 7. Графики решений а — задача 3; б — задача 4 0 0.5 1 0 0.25 0.5 0.75 x 0 0.25 0.5 0.75 y u 0 0.5 1 0 0.25 0.5 0.75 x 0 0.25 0.5 0.75 y u а б Рис. 8. Графики решений а — задача 5; б — задача 6 137 Задача 5: f (x, y) = sin 2 (πxy), µ(0, y) = µ(1, y) = sin(πy), µ(x, 0) = µ(x, 1) = x(1 − Задача 6: f (x, y) = arctan µ x + 1 y + 1 ¶ , µ(0, y) = µ(1, y) = 0, µ(x, 0) = sin 2 (πx), µ(x, 1) = cosh [x(1 − x)] − Задача 7: f (x, y) = |2x − y|, µ(0, y) = µ(1, y) = y(1 − y), µ(x, 0) = | sin(2πx)|, µ(x, 1) = | sin(2πx)|e 2x . 0 1 2 3 4 0 0.25 0.5 0.75 x 0 0.25 0.5 0.75 y u n z n 100 200 300 400 500 10 -4 10 -3 10 -2 10 -1 1 2 4 3 5 5 6 6 а б Рис. 9. а — график решения задачи 7; б — поведение нормы разности ku n+1 − u n k C ; метод ПВР для задачи 7 при различных значениях параметра верхней релаксации τ = 1.5 (1), τ = 1.6 (2), τ = 1.7 (3), τ = 1.8 (4), τ = 1.9 (5), τ = 1.95 (6). Размер сетки N x = N y = 50 12.3. Задания к практическим занятиям. Для каждого из двух указанных в задании итерационных методов определить экспериментально оптимальное значение τ опт итерационного параметра τ , при котором итерации сходятся наиболее быстро. В расчетах использовать квадратную сетку с шагом h. Выяснить характер зависимости τ опт от Экспериментально проверить сходимость разностной схемы на последовательности измельчающихся сеток. На основе численных экспериментов выявить, какой из двух итерационных методов предпочтительнее для рассматриваемого класса тестовых задач. В процессе решения задачи выводить в файл данные, характеризующие сходимость итерационного процесса, например, номер итерации и норму разности двух последовательных итераций z n = ku n+1 − Эту же информацию можно визуализировать (см. рис. 9, б. Если известно точное решение u(x, y) тестовой задачи, то дополнительно в файл выводить и норму ku n − отклонений итераций от точного решения. Графическая информация должна включать графики численного решения на каждой итерации водном или нескольких заданных сечениях или y = Задание 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 заданий, с го по е, заключаются в сравнении итерационных методов, перечисленных в первых 12 заданиях, ноне на задаче Дирихле для уравнения Пуассона, а на нелинейной задаче (7.13) для уравнений метода эквираспределения. 139 Ответы, указания, решения. ν τ h 2 = 1 6 1.2. ν τ h 2 = 1 6 ; ϕ n j = f (x j , t n ) + τ 2 f t (x j , t n ) + ν τ 2 f xx (x j , t n ). 1.3. σ = 1 2 − √ 5 6 ; τ = h 2 Решение. Пусть u(x, t) — решение однородного уравнения теплопроводности Рассмотрим на этом решении локальную погрешность аппроксимации в узле (x j , t n ) дифференциального уравнения (13.1) разностным уравнением 2 u xxtt + + h 2 12 u xxxx + τ h 2 12 u xxxxt + O(τ 2 h 2 ) + h 4 360 u xxxxxx + O(τ h 4 ) ´ − −(1 − σ)ν µ u xx + h 2 12 u xxxx + h 4 360 u xxxxxx ¶ + O(τ 3 + В последнем равенстве все производные вычисляются водном и том же узле (x j , t n ). Далее для сокращения записи формул этот узел указываться не будет. После приведения подобных получим u t − νu xx + τ 2 u tt + τ 2 6 u ttt − σντ µ u xxt + τ 2 u xxtt + h 2 12 u xxxxt ¶ − −ν µ h 2 12 u xxxx + h 4 360 u xxxxxx ¶ + O(τ 3 + h 6 + τ 2 h 2 + τ h 4 ). 140 Учитывая, что для решения уравнения (13.1) выполняются равенства перепишем выражение для погрешности аппроксимации в следующем виде σν 2 τ − ν h 2 12 ¶ u xxxx + + µ τ 2 6 ν 3 − τ σν 2 h 2 12 − σν 3 τ 2 2 − ν h 4 360 ¶ u xxxxxx + O(τ 3 + h 6 + τ 2 h 2 + τ Пусть = 1 2 − h 2 12τ Тогда коэффициент при обратится в нуль и, следовательно 20 − τ 2 ν 2 ¶ u xxxxxx + O(τ 3 + h 6 + τ 2 h 2 + τ Если теперь приравнять нулю коэффициент прите. взять то для погрешности аппроксимации будем иметь выражение O(h 6 ). 1.4. O(τ 2 + Указание. Покажите, что разностное уравнение схемы (можно переписать в следующем виде 6τ ¡ h 2 Λu n+1 j − h 2 Λu n j ¢ + u n+1 j − u n j τ = νΛ u n+1 j + или u n j τ = νΛ ·µ 1 2 − h 2 6τ ν ¶ u n+1 j + µ 1 2 + h 2 6τ ν ¶ u n j ¸ . (13.4) 141 Таким образом, схема (1.102) является схемой с весами, причем = 1 2 − h 2 6τ ν . 1.5. O(τ 2 + Указание. Покажите, что схема (1.103) является схемой с весами повышенного порядка аппроксимации. Указание. При использовании трехточечного шаблона для аппроксимации производной u x (0, t) в левом граничном узле явная схема с порядком аппроксимации O(τ + h 2 ) запишется так u n j τ = ν u n j−1 − 2u n j + u n j+1 h 2 + f n j , j = 1, . . . , N − 1, −3u n 0 + 4u n 1 − u n 2 2h + γu n 0 = µ 0 (t n ), u n N = µ l (t n ), n = 0, . . . , M, u 0 j = u 0 (x j ), j = 0, . . . , Алгоритм получения решения на временном слое включает все- бя два шага) вычисление с помощью разностного уравнения значений во внутренних узлах сетки x j , j = 1, . . . , N − 1; 2) определение значений u n+1 ив граничных узлах x 0 = и x N = l на основе разностных уравнений в этих узлах с использованием при необходимости значений во внутренних узлах, уже вычисленных на первом шаге. Например 0 = 4u n+1 1 − u n+1 2 − 2hµ 0 (t n+1 ) 3 − Условие γ ≤ 0 исходной задачи (1.104) гарантирует, что 3 − 2hγ 6= 0 при любом Теперь для аппроксимации той же производной u x (0, t) будем использовать двухточечный шаблон. Чтобы получить на решении задачи аппроксимацию второго порядка по h, потребуем выполнения дифференциального уравнения вплоть до границы x = 0. Тогда, t) = 1 ν (u t (0, t) − f (0, t)) (13.7) 142 и разностная схема u n j τ = ν u n j−1 − 2u n j + u n j+1 h 2 + f n j , j = 1, . . . , N − 1, u n 1 − u n 0 h − h 2ν µ u n+1 0 − u n 0 τ − f n 0 ¶ + γu n 0 = µ 0 (t n ), u n N = µ l (t n ), n = 0, . . . , M, u 0 j = u 0 (x j ), j = 0, . . . , будет аппроксимировать задачу (1.104) с порядком O(τ + Алгоритм решения по этой схеме совпадает с алгоритмом реализации схемы (13.5), только вместо (13.6) теперь используем другую формулу. Указание. Соответствующие полностью невные схемы имеют вид u n j τ = ν u n+1 j−1 − 2u n+1 j + u n+1 j+1 h 2 + f n+1 j , j = 1, . . . , N − 1, −3u n+1 0 + 4u n+1 1 − u n+1 2 2h + γu n+1 0 = µ 0 (t n+1 ), u n+1 N = µ l (t n+1 ), n = 0, . . . , M − 1, u 0 j = u 0 (x j ), j = 0, . . . , и u n j τ = ν u n+1 j−1 − 2u n+1 j + u n+1 j+1 h 2 + f n+1 j , j = 1, . . . , N − 1, u n+1 1 − u n+1 0 h − h 2ν µ u n+1 0 − u n 0 τ |