2.1 Аппроксимация

Для того чтобы написать разностную схему, приближенно описывающую данное дифференциальное уравнение, нужно:

2.1.1 Разностная аппроксимация простейших дифференциальных операторов

Заменим область Ω изменения аргумента x сеткой ωh, т.е. конечным множеством точек xi. Вместо функции u(x) непрерывного аргумента x Ω будем рассматривать сеточные функции y(xi), которую можно представить в виде вектора. Если пронумеровать все узлы в некотором порядке x1,x2,...,xN, то значения сеточной функции в этих узлах можно рассматривать как компоненты вектора

(y1,y2,...yN).

Дифференциальный оператор L, заданный в классе функций непрерывного аргумента, может быть приближенно заменен (аппроксимирован) разностным оператором Lh, заданным на сеточных функциях. Для этого каждая из производных заменяется разностным отношением, содержащим значения сеточной функции в нескольких узлах сетки. Рассмотрим пример аппроксимации первых и вторых производных функции одного переменного.

Пусть wh = {xi = ih} – сетка с шагом h на отрезке 0 x 1. Рассмотрим производную Lv = vфункции v(x). Заменить ее разностным выражением можно бесчисленным множеством способов. Простейшими являются замены

     vi --vi--1   -
Lv ~    h    = L hvi
(2.1)

– левая разностная производная или левое разностное отношение,

Lv ~ vi+1 --vi= L+ v
        h        h i
(2.2)

– правая разностная производная и

     vi+1 - vi-1   0
Lv ~ ----2h----= L hvi
(2.3)

– центральная разностная производная.

При замене Lv = vразностным выражением Lhvi допускается погрешность Lhvi -(Lv)i = ψih, называемая погрешностью аппроксимации оператора L разностным оператором Lh. Естественно требовать, чтобы при стремлении h к нулю эта погрешность стремилась к нулю. Разложим v(x) по формуле Тейлора в окрестности точки x = xi

                 2
vi1 = vi  hv′i + h-v′′i + O(h3)
                2

(предполагая при этом, что функция v(x) – достаточно гладкая в некоторой окрестности (x - h0,x + h0) точки x и h < h0, h0 – фиксированное число) Подставляя это разложение в (2.1)-(2.3), получим

  vx = vi+1---vi =   v′i + hv′i′ + O (h2),         (2.4)
          h             2
  vx = vi --vi-1 =   v′i - hv′i′ + O (h2),         (2.5)
          h             2
v˚x = vi+1---vi-1- =   v′i + O(h2).               (2.6)
        2h
Отсюда видно, что
ψh  =   v - v′= O(h),                 (2.7)
 ih      x    i′
ψi  =   vx- vi = O(h),                (2.8)
ψhi  =   v˚x - v′i = O(h2).               (2.9)

Будем говорить, что разностный оператор Lh аппроксимирует дифференциальный оператор L с порядком m > 0 в точке x, если max˜w h|ψih| = O(hm). Таким образом, левая и правая разностные производные аппроксимируют Lv = vс первым порядком, а центральная разностная производная – со вторым порядком.

Рассмотрим теперь вторую производную Lv = v′′ = d2v2
dx. Чтобы написать разностную аппроксимацию второй производной, надо использовать трехточечный шаблон, состоящий из узлов xi-1,xi,xi+1, в этом случае

             1-            vi+1---2vi +-vi-1
Lhvi = vxx,i = h (vx,i - vx,i) =    h2
(2.10)

Пользуясь разложением функции v(x) по формуле Тейлора, можно показать, что порядок аппроксимации в этом случае равен двум, т.е.

v  - v′′ = O(h2),
 xx

так как

       ′′  h2 (4)     4
vxx = v + 12v   + O(h ).
(2.11)

2.1.2 Аппроксимация уравнения Пуассона

Рассмотрим разностную схему для решения задачи Дирихле для уравнения Пуассона

     ∑3 ∂2u
Δu =    ∂x2-= - f(x), x ∈ Ω,
     α=1   α
u(x) = μ (x), x ∈ δΩ,
(2.12)

Пусть


Ω = (0 ≤ x1 ≤ l1,0 ≤ x2 ≤ l2,0 ≤ x3 ≤ l3,)
(2.13)

- параллелепипед со сторонами l1,l2,l3, δΩ - его граница. Рассмотрим в Ω = Ω + δΩ задачу Дирихле для уравнения Пуассона:

     ∑p   2
Δu =     ∂-u2-= - f (x),x = (x1,x2,x3) ∈ Ω,u(x) = μ(x),x ∈ δΩ.
     α=1 ∂xα
(2.14)

Построим в Ω сетку ωh с шагами h1 = l1∕N1, h2 = l2∕N2 и h3 = l3∕N3, где N1,N2,N3 > 0 - целые числа. Для этого построим два семейства прямых

xi1 = i1h1,i1 = 0,1,...,N1,
 1
(2.15)

xi22 = i2h2,i2 = 0,1,...,N2,
(2.16)

xi33 = i3h3,i3 = 0,1,...,N3.
(2.17)

Точки пересечения этих прямых x = (i1h1,i2h2,i3h3) назовем узлами. Если x = (i1h1,i2h2,i3h3) лежит внутри параллелепипеда (т.е. 0 < i1 < N1,0 < i2 < N2,0 < i3 < N3), то такой узел назовем внутренним. Пусть ωh - множество всех внутренних узлов.

Узлы, лежащие на границе прямоугольника (при i1 = 0,N1 или i3 = 0,N3 или i3 = 0,N3), кроме восьми узлов (0,0,0), (l1,0,0), (0,l2,0), (l1,l2,0), (0,0,l3), (l1,0,l3), (0,l2,l3), (l1,l2,l3), назовем граничными. Они образуют множество γh = {(i1h1,i2h2,i3h3)}. Совокупность всех внутренних и граничных узлов назовем сеткой ωh = ωh + γh в прямоугольнике Ω. В каждом внутреннем узле x ωh может быть построен семиточечный регулярный шаблон, все узлы которого x1α= 1,2,3 принадлежат ωh (т.е. либо ωh, либо γh). Поэтому во всех внутренних узлах можно заменить оператор Лапласа Δu разностным оператором

Λu = ux1x1 + ux2x2 +u x3x3.
(2.18)

Правую часть - f(x) уравнения (2.12) можно аппроксимировать сеточной функцией - ϕ(x) так, чтобы ϕ(x) - f(x) = O(|h|2),f(x) C2. Считая f(x) непрерывной функцией, полагаем ϕ(x) = f(x).

В результате задаче (2.12) ставим в соответствие разностную задачу Дирихле: найти сеточную функцию y(x), определенную на ωh, удовлетворяющую во внутренних узлах (на ωh) уравнению

Λy = - f(x),  Λy = yx1x1 + yx2x2 + yx3x3,x ∈ ω
(2.19)

и принимающую на границе γh заданные значения

y(x) = μ(x ),  x ∈ γh.
(2.20)

После замены дифференциальных уравнений разностными уравнениями на некоторой сетке ωh мы получаем систему линейных алгебраических уравнений, которую можно записать в матричной форме

ΘY = Φ,
(2.21)

где Θ - квадратная матрица, Y = (y1,y2,...,yN) - искомый вектор, Φ = (ϕ12,...,ϕN) - известная правая часть, включающая и правые части краевых условий.

Каждой матрице Θ можно поставить в соответствие некоторый линейный оператор A, отображающий пространство RN в RN. Тогда уравнение (2.21) примет вид

Ay = ϕ,
(2.22)

где y - искомый, а ϕ - известный векторы пространства RN. Оператор A отображает в себя пространство сеточных функций, заданных на ωh и удовлетворяющих однородным граничным условиям.

Матричная запись для случая нулевых граничных условий. Рассмотрим уравнение Пуассона с нулевыми граничными условиями:

Δu = - f, 0 < x < l
u(x) = 0, x = 0,
u(x) = 0,  x = l.
(2.23)

Проведем конечно-разностную аппроксимацию, и получим следующую систему линейных уравнений, для узлов i = 0,...,N:

-1
h2(yi- 1 - 2yi + yi+1) = - ϕi,i = 1,...,N - 1,
y0 = 0,
y  = 0,
 N
(2.24)

для внутренних узлов i = 1,...,N - 1 система будет выглядеть следующим образом:

  1-
- h2(- 2y1 + y2) = ϕ1,i = 1
  1-
- h2(yi-1 - 2yi + yi+1) = ϕi,i = 2,...,N - 2,
  1-
- h2(- 2yN -1 + yN-2) = ϕN-1,i = N - 1
(2.25)

Запишем полученную систему в матричном виде:

Ay = ϕ,
(2.26)

где

       (                             )
          2   - 1  0   0   ...  0   0
       || - 1  2   - 1  0   ...  0   0 ||
A = -12 ||  0   - 1  2  - 1  ...  0   0 || ,
    h  |(  ...  ...  ...   ...  ... ... ... |)
          0   0    0   0   ... - 1  2

    (       )     (       )
        y1            ϕ1
    ||   y2  ||     ||   ϕ2  ||
y = ||   y3  || ,ϕ = ||  ϕ3  || .
    |(   ...  |)     |(   ...  |)
       yN -1          ϕN- 1

Полученная система имеет симметричную и положительно определенную матрицу (symmetric positive define, SPD).

Матричная запись для случая граничных условий первого рода. Рассмотрим случай уравнение Пуассона с граничными условиями Дирихле

Δu = - f,0 < x < l
u(x) = a,x = 0,
u(x) = b,x = l.
(2.27)

Запишем аппроксимацию уравнения для внутренних узлов

  1-
- h2(a- 2y1 + y2) = ϕ1,i = 1
  1-
- h2(yi-1 - 2yi + yi+1) = ϕi,i = 2,...,N - 2,
  -1
- h2(yN-2 - 2yN-1 + b) = ϕN -1,i = N - 1
(2.28)

Перебросим константы a и b в правую сторону

 -1                 -a
-h2 (- 2y1 + y2) = ϕ1 + h2 ,i = 1
  1-
- h2(yi-1 - 2yi + yi+1) = ϕi,i = 2,...,N - 2,
 -1                        b-
-h2(yN- 2 - 2yN-1) = ϕN -1 + h2,i = N - 1
(2.29)

и запишем полученное СЛАУ с симметричной положительно определенной матрицей в следующем виде

Ay = ϕ,
(2.30)

где

       (                     )
          2  - 1   0  ...  0
       || - 1  2   - 1 ...  0  ||
    -1 ||  0  - 1   2  ...  0  ||
A = h2 ||  ...  ...  ...  ...  ... || ,
       |(  0   0    0  ...  - 1 |)
          0   0    0  ...  2
    (       )     (            )
        y1            ϕ1 + ah2-
    ||   y2  ||     ||     ϕ2     ||
    ||   y3  ||     ||     ϕ3     ||
y = ||   ...  || ,ϕ = ||     ...     || .
    |( yN -2 |)     |(    ϕN- 2   |)
      y              ϕ    + b-
       N -1           N-1   h2

Матричная запись для случая граничных условий второго рода. Пусть на одной границе задано граничное условие Неймана:

Δu = - f,0 < x < l
du-
dx = a,x = 0,
u(x) = b,x = l.
(2.31)

Система линейных уравнений, полученная после конечно-разностной аппроксимации, будет выглядеть следующим образом

--1 (y0 - 2y1 + y2) = ϕ1,i = 1
 h2
- 1-(yi-1 - 2yi + yi+1) = ϕi,i = 2,...,N - 2,
  h2
- -1(yN-2 - 2yN-1 + yN ) = ϕN- 1,i = N - 1,
  h2
(2.32)

где граничные условия аппроксимируются следующим образом

y1 --y0= a,
  h
yN = b,
(2.33)

откуда

y0 = y1 - ah,
yN = b.
(2.34)

Подставляя в уравнения для i = 1, получаем

- -1(- y1 + y2 - ah) = ϕ1, i = 1
  h2
(2.35)

и перебрасываем константу ah в правую сторону

 -1                -a
-h2 (- y1 + y2) = ϕ1 - h , i = 1
- 12(yi-1 - 2yi + yi+1) = ϕi, i = 2,...,N - 2,
  h
  1                        b
--2 (yN -2 - 2yN-1) = ϕN -1 +-2, i = N - 1
 h                         h
(2.36)

Получили СЛАУ с симметричной положительно определенной матрицей

Ay = ϕ,
(2.37)

где

       (  1  - 1   0  ...  0  )
       | - 1  2   - 1 ...  0  |
     1 ||                     ||
A = -2 ||  0  - 1   2  ...  0  || ,
    h  ||  ...  ...  ...  ...  ... ||
       (  0   0    0  ...  - 1 )
          0   0    0  ...  2
    (   y1  )     (   ϕ1 - ah   )
    ||   y2  ||     ||     ϕ2     ||
    ||   y   ||     ||     ϕ      ||
y = ||    3  || ,ϕ = ||      3     || .
    |   ...  |     |     ...     |
    ( yN -2 )     (    ϕN- 2b- )
      yN -1          ϕN-1 + h2

Матричная запись для случая граничных условий третьего рода. Рассмотрим третью краевую задачу

Δu = - f,0 < x < l
du
dx + σu = a,x = 0,
u(x) = b,x = l.
(2.38)

Граничные условия аппроксимируются следующим образом

y1 --y0
  h    + σy0 = a,
yN = b
(2.39)

откуда

y0 = y1∕(1- σh )- ah∕(1 - σh),
yN = b
(2.40)

Запишем конечно-разностную аппроксимацию уравнения

-1
h2((2- 1∕(1- σh))y1 - y2) = ϕ1 - ah∕(1- σh ),i = 1
-1
h2(- yi-1 + 2yi - yi+1) = ϕi,i = 2,...,N - 2,
 1
h2(- yN -2 + 2yN- 1) = ϕN-1 + bh,i = N - 1
(2.41)

Получили СЛАУ с симметричной положительно определенной матрицей

Ay = ϕ,
(2.42)

где

       (                               )
         2 - 1∕(1 - σh ) - 1 0   ...  0
       ||      - 1       2   - 1 ...  0  ||
    -1 ||       0       - 1  2   ...  0  ||
A = h2 ||       ...       ...  ...  ...  ... || ,
       |(       0        0   0   ...  - 1 |)
               0        0   0   ...  2
    (       )     (                 )
    |   y1  |     |  ϕ1 - ah∕(1 - σh )|
    ||   y2  ||     ||        ϕ2       ||
y = ||   y3  || ,ϕ = ||       ϕ3       || .
    ||   ...  ||     ||        ...       ||
    (  yN-2 )     (       ϕN- 2     )
       yN-1             ϕN -1 + bh

2.1.3 Разностные схемы для эллиптических уравнений с переменными коэффициентами

Рассмотрим задачу Дирихле для эллиптического уравнения в области G + δG = G:

Lu = - f(x), x ∈ G,
u = μ(x),  x ∈ δG.
(2.43)

Пусть G = {0 xα lα= 1,2,3} - параллелепипед, а

      3
     ∑
Lu =    Lαu,
     α=1
(2.44)

         (          )
L α =--∂-  kα(x) ∂u-
     ∂x α       ∂xα,
(2.45)

kα(x) - кусочно-непрерывная функция, 0 < c1 kα(x) c2= 1,2,3, где c1 > 0 и c2 > 0 постоянные.

Обозначим ωh = ωh + γh сетку с шагами h1,h2,h3. Оператор Lα аппроксимируем разностным оператором

                           (         )
Λαu = (aαux )  ≈ Lαu = -∂--  kα(x ) ∂u- ,
           α xα        ∂x α       ∂xα
(2.46)

где aα(x) выбирается так, чтобы

Λαu - Lαu = O(h2α),
(2.47)

0 < c1 ≤ aα(x) ≤ c2,
(2.48)

т.е. Λα имеет второй порядок аппроксимации, например,

     - 0.5α
aα = kα
(2.49)

или

a1(x1,x2,x3) = k1(x1 - 0.5h1,x2,x3),
a2(x1,x2,x3) = k2(x1,x2 - 0.5h2,x3),
a (x,x ,x ) = k (x ,x ,x - 0.5h ).
 3  1 2  3    3  1  2  3     3
(2.50)

Оператору Lu ставим в соответствие оператор Λu = α=13Λαu, а задаче (2.43) – разностную задачу Дирихле

Λy = - ϕ,  x ∈ ωh,
y = μ(x), x ∈ γh,
Λy = (Λ1 + Λ2 + Λ3)y, Λαy = (aαyx )x ,
                               α  α
(2.51)

имеющую погрешность аппроксимации (невязку) ψ = Λu + ϕ = O(|h|2).

Лемма. Оператор

              3
Ay = - Λy = - ∑ Λ y,  y ∈ H.
             α=1 α
(2.52)

самосопряжен и положительно определен и для него верны оценки

c1(A˙v,v) ≤ (Av, v) ≤ c2(A˙v, v),
      2                 2
c1δ||v||≤ (Av,v) ≤ c2Δ ||v||
(2.53)

для любого v H, где

 ˙     ˙     ∑3
Ay = -Λy = -    yxαxα,
             α=1
(2.54)

а δ > 0,Δ > 0 - постоянные, определяемые согласно

          ∑3
λmin = δ =    4--sin2 πh-α≥ 8-+ 8-+ 8-= δ0
          α=1 h2α     2lα    l21   l22   l23
           ∑3 -4-  2 πhα-  -4   4-   4-
λmax = Δ =    h2α cos 2lα < h21 + h22 + h23 = Δ0
           α=1
(2.55)

Неравенство (2.53) можно записать в операторной форме

  ˙         ˙
c1A ≤ A ≤ c2A,
c1δE ≤ A ≤ c2ΔE,
(2.56)

где E- единичный вектор.

Пример матричной записи. Рассмотрим эллиптическое уравнение с переменными коэффициентами

∇ (k∇u ) = - f,0 < x < l
 du-
kdx = a,x = 0,
u(x) = b,x = l.
(2.57)

Запишем конечно-разностную аппроксимацию для i = 1,..,N - 1

-1(- ki-1∕2yi-1 +(ki-1∕2 + ki+1∕2)yi - ki+1∕2yi+1) = ϕi,
h2
(2.58)

и для граничных узлов

k1∕2y1 --y0 = a, yN = 0.
      h
(2.59)

Поскольку

k1∕2(y1 - y0) = ah,
k  y = k   y - ah,
 1∕2 0   1∕2 1
(2.60)

то для уравнения (2.58) при i = 1 получим

-1 (- k  y + (k   + k    )y - k    y ) = ϕ - a∕h,
h2    1∕2 1    1∕2   1+1∕2  1   1+1∕2 2    1
 1
h2 (k3∕2y1 - k3∕2y2) = ϕ1 - a∕h.
(2.61)

Для i = N - 1 запишем

1-(- kN -3∕2yN- 2 + (kN-3∕2 + kN-1∕2)yN -1 - kN -1∕2b) = ϕN -1
h2
-1(- kN-3∕2yN -2 + (kN-3∕2 + kN- 1∕2)yN-1) = ϕN -1 + kN -1∕2b∕h2.
h2
(2.62)

В матричном виде полученная система будет выглядеть следующим образом

Ay = ϕ,
(2.63)

       (                                                 )
       |  k3∕2    - k3∕2        0      ...        0        |
     1 || - k3∕2  k3∕2 + k5∕2   - k5∕2    ...        0        ||
A = h2 ||   0      - k5∕2    k5∕2 + k7∕2 ...        0        || ,
       (   ...       ...         ...     ...        ...       )
           0        0          0      ...  kN-3∕2 + kN-1∕2
    (       )      (               )
        y1               ϕ1 - ah
    ||   y2  ||      ||       ϕ2      ||
y = ⋅||  y3  || ,ϕ = ||       ϕ3      || .
    |(   ...  |)      |(       ...      |)
                            kN-1∕2b
       yN-1          ϕN -1 +  h2
Матрица у представленной СЛАУ является симметричной и положительно определенной.

2.1.4 Разностные схемы для стационарного уравнения конвекции-диффузии

Рассмотрим стационарное уравнение конвекции-диффузии в недивергентном виде:

                        (       )
 ∑  v (x)-∂u--  ∑   -∂-- k(x) ∂u- = f(x),c ∈ Ω
α=1,2 α   ∂xα   α=1,2 ∂xα      ∂xα
(2.64)

в прямоугольнике

Ω = {x|x = (x1,x2),0 < x α < lα,α = 1,2}.
Это уравнение дополняется простейшим граничным условием
u (x) = 0,x ∈ ∂Ω.
Также рассмотрим уравнение конвекции-диффузии в дивергентной форме:
 ∑    ∂            ∑    ∂  (     ∂u )
     ∂x--(vα (x)u)-      ∂x-- k(x)∂x-- = f (x),c ∈ Ω.
α=1,2  α          α=1,2  α        α
(2.65)

Стационарную задачу конвекции-диффузии запишем в виде операторного уравнения

Au = f,  A = C + D.
(2.66)

Здесь D - оператор диффузионного переноса, который определяется выражением

        ∑       (        )
Du  = -     -∂-- k(x)-∂u- .
       α=1,2∂xα      ∂xα
(2.67)

Оператор конвективного переноса записывается в соответствии различным формам. Для оператора конвективного переноса в недивергентной форме положим C = C1, где

       ∑
C1u =      vα(x)-∂u-
      α=1,2    ∂x α
(2.68)

Аналогично имеем для дивергентной формы C = C2

       ∑
C2u =      -∂--(vα(x)u).
      α=1,2∂xα
(2.69)

Оператор диффузионного переноса D самосопряжен и положительно определен (SPD):

      *
D  = D ,  D ≤ k1δE,
(2.70)

где δ,Δ - минимальное и максимальное собственное значение оператора Лапласа, k(x) : k1 < k < k2 и справедлива оценка:

γ1E ≤ D ≤ γ2E,
(2.71)

γ1 = k1δ, γ2 = k2Δ,
δ =-82+  82,
   l1   l2
Δ = 4( 12+-12 ).
      h1  h2

Операторы конвективного переноса в дивергентной и недивергентной (C1,C2) формах являются ограниченными:

|(Cαu,u)| ≤ M1 ||u||2,
M  =  1||divv||   ,  ||w ||    = max    |w (x)|.
  1   2      C(Ω)      C(Ω )     x∈Ω
(2.72)

Приведем также оценки подчиненности оператора конвективного переноса оператору диффузионного переноса:

||Cu ||2 ≤ M2 (Du, u),
              2       2
(C = C1)M2 = k1max {||vα||C(Ω)},
             -2         2                2          1
(C = C2)M2 = k1(2max {||vα||C(Ω )} + M0||divv||C(Ω)),M0  = δ,
(2.73)

Рассмотрим в области Ω стационарное уравнение конвекции-диффузии, которое соответствует дивергентной записи конвективных слагаемых

∑3               ∑3   2
    -∂-(vα(x)u)-    k∂-u-= f(x),x ∈ Ω,
α=1 ∂xα          α=1 ∂x2α
(2.74)

где vα - компоненты скорости v = (v1,v2,v3), которые мы считаем заданными.

Уравнение (2.74) дополняется обычными условиями первого рода на границе

u(x) = g(x),x ∈ δΩ.
(2.75)

Для приближенного решения задачи введем равномерную разностную сетку ω с шагом h. Приведем разностную схему с центральными разностями

∑3            ∑3
   (bα (x)y)˙xα -   kyxαxα = ϕ(x ),x ∈ ω.
α=1           α=1
(2.76)

Семиточечное разностное уравнение во внутренних узлах сетки будем записывать в следующей форме

A (x)y(x) = B+1 (x)yi+1,j,k + B-1 (x )yi-1,j,k
   +             -
+B 2 (x )yi,j+1,k + B2 (x)yi,j-1,k
+B+3 (x)yi,j,k+1 + B -3 (x)yi,j,k-1 + F(x),x ∈ ω,
(2.77)

где

A (x) = 2- + 2-+ -2 ,
       h21   h22  h23
B+ (x) = 1-- bx1+h1,x2,x3,B- (x) = 1- + bx1-h1,x2,x3,
  1      h21      2h      1      h21      2h
         1   bx ,x +h ,x          1    bx,x-h ,x
B+2 (x) = h2---1-22h-2-3,B-2 (x) = h2 + -1-22h-2-3,
          2                      2
B+ (x) = 1-- bx1,x2,x3+h3,B- (x) = 1- + bx1,x2,x3-h3,
  3      h23      2h      3      h23      2h
F (x ) = f(x).
(2.78)

Обычные достаточные условия выполнения принципа максимума для разностной схемы имеют вид

          +-
A(x) > 0,Bα-(x) > 0,
D(x) ≥ 0,
(2.79)

где

              3          3
D (x ) = A(x)- ∑ B+ (x)- ∑  B -(x).
             α=1  α     α=1  α
(2.80)

Условие D(x) 0 соответствует диагональному преобладанию по строкам. Эти условия для схемы с центральными разностями будут выполнены лишь при достаточно малых шагах сетки (малых скоростях в задачах конвективного переноса)

    -1
hαk   maxx∈ω|bα(x)| < 2,α = 1,2,3.
(2.81)

Безусловно монотонные разностные схемы можно построить при использовании аппроксимации первого порядка направленными разностями. Введем обозначения

bα(x) = b+α(x)+ b-α(x),
b+α(x) = 1 (bα(x)+ |bα(x)|) ≥ 0,b-α(x) = 1(bα(x)- |bα(x)|) ≤ 0,
       2                          2
(2.82)

так что b+(x) = b(x), если b(x) > 0 и b-(x) = 0 при b(x) < 0.

Будем рассматривать следующее разностное уравнение

∑3  +         ∑3  -
  (bα(x)y)xα +    (bα (x)y)xα
α=1           α=1
  ∑3
-    yxαxα = ϕ(x),x ∈ ω.
 α=1
(2.83)

Разностное уравнение во внутренних узлах сетки запишем в следующей форме

A (x)y(x) = B+1 (x)yi+1,j,k + B-1 (x )yi-1,j,k
+B+  (x )y      + B-(x)y
   2    i,j+1,k    2    i,j-1,k
+B+3 (x)yi,j,k+1 + B -3 (x)yi,j,k-1 + F(x),x ∈ ω,
(2.84)

где

      ∑3      ∑3   +     ∑3  -
A(x) =    22+     bα(x)-    bα(x),
      α=1 hα  α=1  hα    α=1  hα
              -                       +
B+ (x) = 1-- bx1+h1,x2,x3,B- (x) = 1- + bx1-h1,x2,x3,
  1      h21      2h      1      h21      2h
         1   b-x ,x +h ,x          1    b+x,x-h ,x
B+2 (x) = h2---1-22h-2-3,B-2 (x) = h2 + -1-22h-2-3,
          2                      2
              -                       +
B+ (x) = 1-- bx1,x2,x3+h3,B- (x) = 1- + bx1,x2,x3-h3,
  3      h23      2h      3      h23      2h
F (x ) = f(x).
(2.85)

При таких коэффициентах достаточные условия монотонности выполнены при любых шагах сетки.

2.1.5 Разностные схемы для нестационарных уравнений

Для выяснения методов построения разностных схем в случае нестационарных задач, а также методов их исследования рассмотрим первую краевую задачу для уравнения

∂u-  ∂2u-
∂t = ∂x2 +f,
(2.86)

в прямоугольнике

D= (0 ≤ x ≤ 1,0 ≤ t ≤ T).

Требуется найти непрерывное в D решение u = u(x,t) задачи

 ∂u   ∂2u
 ∂t-= ∂x2 + f(x,t),0 < x < 1,0 < t ≤ T,
      u(x,0) = u(x),0 ≤ x ≤ 1,
               0
u (0,t) = u1(t),u(1,t) = u2(t),0 ≤ t ≤ T.
(2.87)

Введем сетки

wh = xi = ih,i = 0,1,...,N,wτ = tj = jτ,j = 0,1,...,j0

и сетку в D:

whτ = wh × wτ = (ih,jτ),i = 0,1,...,N,j = 0,1,...,j0

с шагами h = 1∕N и τ = T∕j0. Обозначим через yij значение в узле (xi,tj) сеточной функции y, определенной на w. Заменяя производную ∂u∕∂t первой разностной производной, а 2u∕∂x2 – второй разностной производной uxx = Λu и вводя произвольный вещественный параметр σ, рассмотрим однопараметрическое семейство разностных схем

 j+1   j     (               )
yi---- yi = Λ σyj+1+ (1- σ)yj  + φj,0 < i < N, 0 ≤ j < j0
    τ           i           i     i
(2.88)

Схему (2.88) будем называть иногда схемой с весами.

Краевые и начальные условия аппроксимируем точно:

 j    j  j    j
y0 = u1,yN = u 2,
(2.89)

y0 = y(x,0) = u (x ).
 i      i     0  i
(2.90)

Здесь φij – сеточная функция, аппроксимирующая правую часть f уравнения (2.87), например,

φji = f(xi,tj+0,5),tj+0,5 = tj + 0,5τ,

                             2
Λyi = yxx,i = (yi- 1 - 2yi + yi+1)∕h

От выбора параметра σ, зависят точность и устойчивость схемы (2.88). Рассмотрим схемы, соответствующие частным значениям σ. При σ = 0 получаем четырехточечную схему

 j+1    j
yi-----yi = Λyji + φji,
    τ

или

 j+1          j     j     j       j
yi  = (1- 2γ)yi + γ(yi- 1 + yi+1)+ τφ i,γ = τ∕h2,
(2.91)

определенную на шаблоне (xi,tj+1),(xi,tj),(xi1,tj). Значение yij+1 в каждой точке слоя t = tj+1 (нового слоя) выражается по явной формуле (2.91) через значения yij на слое t = tj (на старом слое).Так как при t = 0 задано начальное значение yi0 = u0(xi), то формула (2.91) позволяет последовательно определить значения y на любом слое. Схема (2.91) называется явной.

Если σ0, то схема (2.88) неявной двухслойной схемой. При σ0 для определения yij+1 на новом слое получаем систему алгебраических уравнений

    j+1  1-j+1      j
σ Λyi  - τyi   = - Fi ,
F j=  1yj+ (1- σ)Λyj+ φj,
  i   τ i          i    i
(2.92)

i = 1,...,N - 1

с краевыми условиями y0j+1 = u1j+1,yNj+2 = u2j+1.

При σ = 1 имеем схему с опережением или чисто неявную схему

 j+1    j
yi-----yi = Λyj+1 + φj.
    τ        i      i
(2.93)

В явных схемах налагаются строгие ограничения на величины временных шагов, а неявные схемы безусловно устойчивы.

2.1.6 Нелинейные уравнения

Для решения нелинейных уравнений использование явных схем нецелесообразно. Условие устойчивости явной схемы требует мелкого шага по времени. Поэтому применяются безусловно устойчивые неявные схемы с итерационными процессами. Итерационными методы предполагают задание достаточно близких к искомому решению начальных данных.

Рассмотрим различные неявные схемы для одномерного квазилинейного параболического уравнения

∂u    ∂ (    ∂u )
---= --- k (u)---  +f (u),0 < x < 1,0 < t ≤ T,
∂t   ∂x      ∂x
   u(x,0) = u0(x),u(0,t) = u1(t),u(1,t) = u2(t),
(2.94)

где k(u) > 0.

Схема a):

          [                             ]
ˆyi --yi = 1-a  (y)ˆyi+1---ˆyi- a(y)ˆyi --ˆyi--1 + f(y)
   τ     h  i+1      h       i     h          i
(2.95)

Схема б):

          [                             ]
ˆyi --yi  1-       ˆyi+1---ˆyi      ˆyi --ˆyi--1
   τ   = h ai+1(ˆy)   h    - ai(ˆy)   h     + f(ˆyi)
(2.96)

где ŷi = yij+1,yi = yij,ai(v) = a(vi-1,vi), например,

ai(v) = 0,5[k(vi-1)+ k(vi)],
(2.97)

        ( vi-1 + vi)
ai(v) = k  ---2---- ,
(2.98)

       2k(vi--1)k(vi)-
ai(v) = k(vi-1)+ k(vi).
(2.99)

Сравним схемы (2.95) и (2.96). Погрешность аппроксимации этих схем O(τ + h2). Обе они абсолютно устойчивы. Схема а) линейна относительно значения функции yj+1 на слое tj+1, и значения функции yj+1 находятся по значению функции yj на слое tj. Поскольку схема абсолютна устойчива шаг τ выбирается из соображений точности. Схема б) нелинейна относительно функции yj+1 и для нахождения ее решения используется метод итераций. Итерационный процесс строится следующим образом:

(s+1)       ⌊       (s+1)  (s+1)       (s+1) (s+1)⌋
-yi---yi= 1-⌈ai+1 ((ys)) yi+1----yi-- ai (s(y)) -yi---yi-1⌉ + f((sy)i).
   τ     h             h                h
(2.100)

Относительно (s+1)
 y разностная схема оказывается линейной. В качестве начальной итерации берется функция y предыдущего шага по времени: (0y) = yj. Итерационный процесс для большинства встречающихся на практике коэффициентов k и f сходится. Практически оказывается достаточным сделать две-три итерации. Даже в том случае, если итерации не сходятся, для повышения точности схемы оказывается полезным сделать две итерации. При счете по итерационной схеме (2.96),(2.100) задают число итераций или точность сходимости итераций ε и требуют выполнения условия

    (s+1)  (s)
maix|  yi - yi | ≤ ε.

Недостаток схемы (2.96),(2.100) в том, что счет итераций требует удвоения числа занимаемых в память ячеек по сравнению со схемой а), т.к. для вычисления (s+1)
 y нужно помнить (s)
y и y. Для нахождения значения yj+1 по yj по схеме(2.96),(2.100) нужно сделать несколько итераций, а ппри счете по схеме а) значение yj+1 находится сразу.

Поскольку обе схемы абсолютно устойчивы и имеют одинаковый порядок аппроксимации, то казалось бы, что схема а) имеет преимущества перед итерационной схемой б). Но, как показывает практика, для получения одинаковой точности счета схема б) позволяет использовать настолько более крупный шаг по времени, что, несмотря на необходимость итераций, это приводит к уменьшению объема вычислительной работы.

Можно использовать схемы, имеющие второй порядок аппроксимации по пространству и времени. Однако такие схемы немонотонны, что приводит часто к появлению “ряби”, поэтому приходится выбирать достаточно мелкий шаг по времени.

Для решения задачи (2.96), помимо метода (2.100), можно применять метод Ньютона.

Рассмотрим метод Ньютона на следующей задаче. Пусть задана функция f(x) = 0 действительного переменного. Требуется найти корни уравнения

f(x) = 0,
(2.101)

где начальное приближение x0 известно. Заменим f(x) отрезком ряда Тейлора

f(x) ≈ H1(x) = f(x0)+ (x- x0)f′(x0)

и за следующее приближение x1 возьмем корень уравнения H1(x) = 0, т.е.

         -f(x0)-
x1 = x0 - f′(x0).

Вообще, если итерация xk известна, то следующее приближение xk+1 в методе Ньютона определяется по правилу

           -f(xk)
xk+1 = xk - f ′(xk),k = 0,1,...
(2.102)

Метод Ньютона называют также методом касательных, так как новое приближение xk+1 является абциссой точки пересечения касательной, проведенной в точке (xk,f(xk)) к графику функции f(x), с осью Ox.

Метод имеет квадратичную сходимость, т.е. в отличие от линейных задач погрешность на следующей итерации пропорциональна квадрату погрешности на предыдущей итерации: xk+1  - x*  = O((xk - x*)2). Быстрая сходимость метода Ньютона гарантируется лишь при очень хороших, т.е. близких к точному решению, начальных приближениях. Если начальное приближение выбрано неудачно, то метод может сходиться медленно, либо не сойдется вообще.