Численные методы. Соснин (2005) (1160462), страница 21
Текст из файла (страница 21)
Взяв σ =и t = tn+ 21 , мы можем получить такую оценку на порядок2аппроксимации:n+ 12ψj= O(τ 2 + h2 ).В противном случае оценка будет похуже:ψjn = O(τ + h2 ).Наконец, общая формула для получения y будет такова:n+1n+1Aj yj+1− Cj yjn+1 + Bj yj−1= −Fj .При этом правая часть уравнения зависит только от yn . Это означает, что, решая данное уравнениепослойно (с n = 0 ) и используя начальные условия, мы можем найти все yjk . Все эти умозаключенияабсолютно идентичны тем, что были в конце предыдущего параграфа.Разностная схема для нелинейного уравнения теплопроводностиИсследуем случай, когда уравнение теплопроводности в краевой задаче имеет такой вид:∂∂u∂u=k(u, x, t)+ f (u, x, t).∂t∂x∂x— это означает, что коэффициенты при производных зависят еще и от искомой функции u, причем,вообще говоря, нелинейно.В этом случае рекомендуется использовать неявные разностные схемы, так как они чаще всегоабсолютно устойчивы.
Приведем пример:n+1yjn+1 − yjn= ayx x,j + f (yjn+1 ),τпри этом коэффициенты a зависят еще и от y. Расписав разностные производные, получим такуюформулу:"#n+1n+1yjn+1 − yjnyj+1− yjn+1yjn+1 − yj−11n+1n+1=aj+1 (yj ) ·− aj (yj ) ·+ f (yjn+1 ).τhhh— для каждого слоя это система нелинейных уравнений относительно yjn+1 . Решается она итераци(k)онным методом следующего вида.
Если обозначить за yj приближение для yjn+1 , то формула дляполучения следующего приближения будет такова:(k+1)(k+1)(k+1)(k+1)− yjn(k)(k)(k)1yj+1 − yjyj − yj−1 =aj+1 ( yj ) ·− aj ( yj ) ·+ f ( yj )τhhh(k+1)yj— это уже система линейных уравнений с трехдиагональной матрицей. Она, в свою очередь, решаетсяметодом прогонки, и мы получаем (k + 1)-е приближение к yjn+1 . Обычно ограничиваются пятьюприближениями:(5)yjn+1 = yj .5.13. РАЗНОСТНАЯ СХЕМА ДЛЯ УРАВНЕНИЯ КОЛЕБАНИЙ5.13109Разностная схема для уравнения колебанийРассмотрим стандартную краевую задачу на уравнение колебаний:utt = uxx + f (x, t), 0 < x < 1, 0 < t 6 T ;0 6 t 6 T; u(0, t) = µ1 (t),u(1, t) = µ2 (t),0 6 t 6 T;u(x, 0) = u0 (x),0 6 x 6 1;ut (x, 0) = ψ(x),0 6 x 6 1.(5.39)Решение такой задачи существует и единственно.Введем дискретную сетку на рассматриваемой области:ωk =xj = jh, h = N1 , j = 0, N ;T, n = 0, K .ωτ =tn = nτ, τ = KЗаметим, что начальные условия дают нам значения искомой функции на границе прямоугольника.Будем использовать уже привычные нам обозначения: u(xj , tn ) = unj — точное решение в узлахсетки, yjn — искомое приближение.Выпишем аналоги краевых и начальных условий: ny= µ1 (tn ); 0nyN = µ2 (tn ); 0yj = u0 (xj ).Сопоставим уравнению дискретный аналог в каждом узле.
Будем использовать шаблон из пятиточек (это минимально необходимый шаблон для аппроксимации вторых производных по t и x).Приближения построим так:1nutt ≈ utt,j= 2 un+1− 2unj + un−1;jjτ1uxx ≈ unxx,j = 2 unj+1 − 2unj + unj−1 .hТогда получим следующий дискретный аналог исходной задачи:nnytt,j= yxx,j+ ϕnj , j = 1, N − 1, n = 1, K − 1;y0n = µ1 (tn ),n = 1, K;nyN = µ2 (tn ),n = 1, K;(5.40)00yj = u (xj ),j = 0, N ;10 yj − yj = ψ(xj ),j = 0, N .τШаблон будет выглядеть следующим образом:Стоит отметить, что для того, чтобы разностная схема была сбалансированной (то есть, чтобыне делать лишних вычислений в одном месте, а потом загрублять их в другом), необходимо, чтобыпорядки аппроксимации в уравнении и в краевом условии второго типа были согласованы.
Иначе,если использовать низкий (первый) порядок аппроксимации в краевом условии, то вся схема будетаппроксимировать исходную задачу с первым порядком аппроксимации.Постараемся все данные в задаче приблизить со вторым порядком точности. Для этого разложимu1j в ряд Тейлора в точке (xj , 0) с остаточным членом второго порядка малости:u1j − u0jτ= u0t,j + u0tt,j + O(τ 2 ).τ2110Глава 5. ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ КРАЕВЫХ ЗАДАЧПостараемся избавиться от u0tt,j . Пусть на границе также выполняется уравнение колебаний, тогдаutt = uxx + f (x, 0), и краевое условие примет вид:τu1j − u0j= u0t,j + u0xx,j + f (xj , 0)+ O(τ 2 ),τ2где выражение u0xx,j + f (xj , 0) точно вычислимо из начальных условий.Тогда разностную схему можно переписать в виде:nnytt,j+ ϕnj ,j = 1, N − 1, n = 1, K − 1;= yxx,jy0n = µ1 (tn ),n = 1, K;nyN = µ2 (tn ),n = 1, K;yj0 = u0 (xj ),j = 0, N ; 10 yj − yj = ψ(xj ) + τ u0 (xj ) + f 0 , j = 0, N .jτ2 xxПорядок аппроксимации краевого условия будет O(τ 2 ).
Потребуем в задаче второй порядок аппроксимации и в первом уравнении. Для этого достаточно взять ϕnj = f (xj , tn ) + O(h2 + τ 2 ).Таким образом, разностный аналог исходного ДУ будет его аппроксимировать со вторым порядкомпо τ и по h.Рассмотрим вопрос о поиске решения. Найдем из построенной разностной схемы выражение дляyjn+1 :nyjn+1 = 2yjn − yjn−1 + τ 2 yxx,j+ τ 2 ϕnj .(5.41)Покажем,что значение сеточной функции вычислимо на всей сетке. Рассмотрим формулу (5.41)при n = 1 :τ2 11yj2 = 2yj1 − yj0 + 2 yj+1− 2yj1 + yj−1+ τ 2 ϕ1j .hТак как yj0 известно из начальных условий, а yj1 можем найти из второго краевого условия, то2yj находится явно на всем втором временном слое.
Аналогично можно найти yj3 и так далее, то естьопределить значение искомой сеточной функции на всей сетке (как и раньше, делая это послойно).Как видно из процесса поиска решения, эта схема — явная, и построенное решение будет единственным.Рассмотрим вопрос устойчивости. Будем использовать привычный нам метод гармоник. Запишемоднородное уравнение:τ2 nnyjn+1 − 2yjn + yjn−1 = 2 yj+1− 2yjn + yj−1.hПодставим в него решение вида yjn = q n eijhϕ , тогда получим:τ2τ22 hϕihϕ−ihϕ22q − 2q + 1 = 2 q(e−2+e) ⇐⇒ q − 2 − 4 2 sinq + 1 = 0.hh25.14. РАЗНОСТНАЯ АППРОКСИМАЦИЯ ЗАДАЧИ ДИРИХЛЕ111Решения этого квадратного уравнения будут такими: s2τ2hϕτ22 hϕq1,2 = 1 − 2 2 sin− 1.±1 − 2 2 sin2h2h2Выделим два случая на дискриминант (подкоренное выражение):1. D > 0 :Так как свободный член равен единице, то |q1 ||q2 | = 1. В этом случае (можно показать, что корнине могут быть противоположными по знаку) |q1 | или |q2 | больше единицы, поэтому устойчивостине будет при данных значениях параметров.2.
D 6 0 :Аналогично, |q1 ||q2 | = 1. Найдем значение абсолютной величины корня:2|q1 | =τ2hϕ1 − 2 2 sin2h22τ2hϕ+ 1 − 1 − 2 2 sin2h22= 1.Тогда |q1 | = |q2 | = 1, и, следовательно, решение будет устойчивым.Распишем условие неположительности дискриминанта:|1 − 2hϕτ2sin2| 6 1,h22или, расписав:τ2hϕsin26 1.h22Правое неравенство выполнено всегда, перепишем левое:−1 6 1 − 2τ2hϕsin26 1.2h2τ2τ6 1, что эквивалентно6 1 (напомним, что в уравнении теплоh2h1τпроводности мы получали ограничение вида 2 6 ). Получается, что построенная схема являетсяh2условно устойчивой.Итак, можно выделить несколько особенностей для уравнения колебаний:1) дополнительное условие на аппроксимацию (из начального условия);2) решение строится послойно (впрочем, как и для многих других задач);3) условная устойчивость с простыми ограничениями.В худшем случае получим5.14Разностная аппроксимация задачи ДирихлеКратко опишем численное решение задачи Дирихле на уравнение Пуассона.
Пусть область, для которой задано уравнение — прямоугольник в E 2 :ux1 x1 + ux2 x2 = f (x1 , x2 ), 0 < x1 < l1 , 0 < x2 < l2 ;u(0, x2 ) = µ1 (x2 ),0 6 x2 6 l2 ;u(l1 , x2 ) = µ2 (x2 ),0 6 x2 6 l2 ;u(x,0)=µ(x),06 x1 6 l1 ;1 11u(x1 , l2 ) = µ2 (x1 ),0 6 x1 6 l1 .112Глава 5. ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ КРАЕВЫХ ЗАДАЧДискретизируем область:ωhx1ωhx2l1=x1i = ih1 , i = 0, N1 , h1 =;N1l2.=x2j = jh2 , j = 0, N2 , h2 =N2Выглядеть это будет следующим образом (обозначения аналогичны предыдущим рисункам):`2 ×××××××××××××lmnopqrstuvwxyz{|}~NOPQRSTUVWXYZ[\]^_`abcdefg456789:;<=>?@ABCDEFGHIJKLM!"#$%&'()*+,-./0123××××××××××××××××jk××× hi××××××××××`1Запишем дискретные аналоги граничных условий:y0,jyN1 ,jyi,0yi,N2====µ1 (x2j ),µ2 (x2j ),µ1 (x1i ),µ2 (x1i ),j = 0, N2 ;j = 0, N2 ;i = 0, N1 ;i = 0, N1 .(5.42)Производные будем приближать разностными производными со вторым порядком аппроксимации:yx1 x1 ,ij + yx2 x2 ,ij = ϕij ,i = 1, N1 − 1, j = 1, N2 − 1,(5.43)где для достижения второго порядка аппроксимации во всем уравнении ϕij считается по формуле:ϕij = f (x1i , x2j ) + O(h21 + h22 ).Найдем формулу для поиска приближенного решения.
Явно получить выражение для приближенийна следующем слое из (5.43) довольно сложно, поэтому сделаем хитрее.Введем вектор всех неизвестных (заметим, всех, то есть получим огромный вектор), используясамый простой, регулярный способ нумерации:Y = y11 , y21 , . . . , y(N1 −1)1 , y12 , y22 , . . . , y(N1 −1)2 , .
. . , y1(N2 −1) , . . . , y(N1 −1)(N2 −1)T.Аналогично введем вектор правой части:F = ϕ11 , ϕ21 , . . . , ϕ(N1 −1)1 , ϕ12 , ϕ22 , . . . , ϕ(N1 −1)2 , . . . , ϕ1(N2 −1) , . . . , ϕ(N1 −1)(N2 −1)T.Объединив, используя эти обозначения, уравнения (5.42) и (5.43), получим уравнение:AY = F.Из структуры шаблона следует такой вид матрицы A (она будет иметь пять ненулевых «диагона-5.14. РАЗНОСТНАЯ АППРОКСИМАЦИЯ ЗАДАЧИ ДИРИХЛЕ113лей»; i, j, s — некоторые параметры, зависящие от N1 , N2 ):a11 a1200...0a1j00...0 a21 a22a230...00a2j0...0 0 a32aa...000a...0 33343j.
. . . . . . . . . . . . . . . . . . . . . . 0 ...0aaa0...00ajjj1 j(j−1)j(j+1).......................A= ai1 0...0aaa0...00iii(i−1)i(i+1)....................... 0 . . . a(s−2)(i−2)00...0aaa0(s−2)(s−3) (s−2)(s−2) (s−2)(s−1) 0 ...0a(s−1)(i−1) 00...0a(s−1)(s−2) a(s−1)(s−1) a(s−1)s 0 0...0asi00...0as(s−1)ass— так называемая матрица ленточной структуры5 .Мы, кстати, опять получили сеточное уравнение (напомним, подробные указания о том, как ихрешать, можно найти в [2]).5 названав честь ленточного червя (тоже большая гадость).Оглавление1 Решение систем линейных алгебраических уравнений1.1 Прямые методы решения СЛАУ. Метод квадратного корня .1.2 Линейные одношаговые итерационные методы .
. . . . . . .1.3 Сходимость одношаговых стационарных методов . . . . . . .1.4 Оценка погрешности одношаговых стационарных методов .1.5 Попеременно-треугольный итерационный метод . . . . . . .1.6 Чебышевский набор итерационных параметров . . . . . . .
.1.7 Одношаговые итерационные методы вариационного типа . .1.8 Примеры итерационных методов вариационного типа . . . .1.9 Двухшаговые итерационные методы вариационного типа . ..........................................................................................................................................................449131721252931342 Задачи на собственные значения362.1 Поиск собственных значений методом вращений . . . . . . . . . . .