Математика

Метод конечных разностей или метод сеток


                                  ВВЕДЕНИЕ


    Значительнаое число задач физики и техники приводят к  дифференциальным
уравнениям  в  частных   прозводных   (уравнения   математической   физики).
Установившиеся   процессы   различной   физической    природы    описываются
уравнениями эллиптического типа.
    Точные  решения  краевых  задач  для  эллиптических  уравнений  удаётся
получить лишь в частных  случаях.  Поэтому  эти  задачи  решают  в  основном
приближённо.  Одним  из  наиболее  универсальных  и   эффективных   методов,
получивших в  настоящее  время  широкое  распространение  для  приближённого
решения уравнений математической физики, является метод  конечных  разностей
или метод сеток.
    Суть  метода  состоит  в  следующем.  Область  непрерывного   изменения
аргументов,  заменяется  дискретным  множеством   точек   (узлов),   которое
называется  сеткой  или  решёткой.  Вместо  функции  непрерывного  аргумента
рассматриваются функции дискретного аргумента, определённые в узлах сетки  и
называемые сеточными функциями.  Производные,  входящие  в  дифференциальное
уравнение и граничные  условия,  заменяются  разностными  производными,  при
этом краевая задача  для  дифференциального  уравнения  заменяется  системой
линейных или нелинейных алгебраических уравнений  (сеточных  или  разностных
уравнений). Такие системы часто называют разностными схемами.  И  эти  схемы
решаются относительно неизвестной сеточной функции.
    Далее мы будем рассматривать применение  итерационного  метода  Зейделя
для вычисления неизвестной сеточной функции в краевой задаче с  неоднородным
бигармоническим уравнением.



                              ПОСТАНОВКА ЗАДАЧИ


    Пусть у нас есть бигармоническое уравнение :
                                      2
                                    U = f

    Заданное на области G={ (x,y) : 0<=x<=a,  0<=y<=b }. Пусть также заданы
краевые условия на границе области G .


    U     = 0                            Y
           x=0
                                b
    Uxxx       = 0
                   x=0
                                                         G

    Ux         = 0
        x=a
    Uxxx  = 0                                                         0
                                              a      X
          x=a

    U    = 0                      U   = 0
       y=0
                  y=b
    Uy        = 0                      Uxx  + Uyy    = 0
         y=0
                  y=b                  y=b



Надо решить эту задачу численно.
    Для решения будем использовать итерационный метод Зейделя  для  решения
сеточных задач.
    По нашей области G построим равномерные сетки Wx и Wy с шагами hx  и hy
соответственно .
                     Wx={ x(i)=ihx,  i=0,1...N,  hxN=a }
                     Wy={ y(j)=jhy,  j=0,1...M,  hyM=b }
    Множество  узлов  Uij=(x(i),y(j))  имеющих  координаты   на   плоскости
х(i),y(j) называется сеткой в прямоугольнике G  и обозначается :

       W={ Uij=(ihx,jhy),   i=0,1...N,   j=0,1...M,   hxN=a,  hyM=b }

Сетка W очевидно состоит из точек пересечения прямых x=x(i)  и y=y(j).
    Пусть задана сетка W.Множество всех  сеточных  функций  заданных  на  W
образует векторное пространство с определённом  на  нём  сложениемфункций  и
умножением  функции  на  число.  На  пространстве  сеточных  функций   можно
определитьразностные  или  сеточные  операторы.  0ператор   A  преобразующий
сеточную функцию  U  в  сеточную  функцию  f=AU  называется  разностным  или
сеточным  оператором.  Множество  узлов  сетки  используемое  при  написании
разностного оператора в узле сетки называется шаблоном этого оператора.
    Простейшим разностным оператором  является  оператор  дифференцирования
сеточной функции, который порождает разностные производные. Пусть W -  сетка
с шагом  h введённая на R т.е.

                        W={Xi=a+ih, i=0, + 1, + 2...}

Тогда разностные производные первого порядка для сеточной функции   Yi=Y(Xi)
, Xi из W, определяется по формулам :

                     L1Yi  = Yi - Yi-1    , L2Yi=L1Yi+1
                                                               h

и называются соответственно левой и правой производной. Используется так  же
центральная производная :

                        L3Yi=Yi+1 - Yi-1  = (L1+L2)Yi
                                                                          2h
        2

Разностные операторы   A1,  A2,  A3  имеют  шаблоны  состоящие  2х  точек  и
используются  при  апроксимации  первой  производной  Lu=u’   .   Разностные
производные  n-ого порядка  определяются  как  сеточные  функции  получаемые
путём  вычисления  первой  разностной  производной  от  функции,  являющейся
разностной производной n-1 порядка, например :

                      Yxxi=Yxi+1 - Yxi = Yi-1-2Yi+Yi+1
                                          2
                                h                  h

                    Yxxi= Yxi+1-Yxi-1 = Yi-2 - 2Yi+Yi+ 2
                                          2
                             2h                     4h

которые используются при апроксимации  второй  производной.  Соответствующие
разностные операторы имеют 3х точечный шаблон.
    Анологично не представляет труда определить разностные  производные  от
сеточных функций нескольких переменных.
    Аппроксомируем нашу задачу с помощью разностных производных. И применим
к получившейся сеточной задаче метод Зейделя.



                               МЕТОД   ЗЕЙДЕЛЯ


    Одним из способов  решения  сеточных  уравнений  является  итерационный
метод Зейделя.
    Пусть нам дана система линейных уравнений :

                                   AU  = f
или в развёрнутом виде :

                         M
                          aijUj   = fi  ,              i=1,2...M
                               i=1

Итерационный  метод  Зейделя  в  предположении  что  диагональные   элементы
матрицы А=(aij) отличны от нуля (aii<>0) записывается в следующем виде :

                  i               (k+1)        M          (k)
              aijYj       +       aijYj    = fi  ,   i=1,2...M
                 j=1                                 j=i+1
           (k)
где Yj - jая компонента  итерационного  приближения  номера  k.  В  качестве
начального приближения выбирается произвольный вектор.
    Определение (k+1)-ой итерации начинается с i=1

                             (k+1)         M          (k)
                          a11Y1 = -      a1jYj +f1
                                  j=2

                             (k+1)
Так как a11<>0 то отсюда найдём Y1. И для i=2 получим :

                        (k+1)                 (k+1)      M              (k)
                    a22Y2 = - a21Y1 -         a2jYj  + f2
                                           j=3

                  (k+1)                   (k+1)                        (k+1)
            (k+1)
Пусть уже найдены  Y1        , Y2       ...   Yi-1 .    Тогда  Yi  находится
из уравнения :

                 (k+1)                          i-1                    (k+1)
 M                  (k)
             aiiYi   =   -            aijYj     -              aijYj      +
   fi                            (*)
                        j=1                                     j=i+1

Из формулы (*) видно ,  что  алгоритм  метода  Зейделя  черезвычайно  прост.
Найденное по формуле (*) значение Yi размещается на месте Yi.
    Оценим число арифметических действий, которое требуется для  реализации
одного итерационного шага. Если все aij не  равны  нулю,  то  вычисления  по
формуле (*)    требуют   M-1   операций    умножения   и   одного   деления.
Поэтому реализация


                             2
одного шага осуществляется за 2M -  M арифметических действий.
    Если отлично от нуля лишь m элементов,  а  именно  эта  ситуация  имеет
место для сеточных эллиптических уравнений, то на  реализацию  итерационного
шага потребуется 2Mm-M действий т.е. число  действий  пропорционально  числу
неизвестных M.
    Запишем теперь метод Зейделя в матричной форме.  Для  этого  представим
матрицу  A  в  виде  суммы  диагональной,  нижней  треугольной   и   верхней
треугольной матриц :

                                A = D + L + U

где

          0           0         .         .         .                     0
                0    a12  a13  .    .    . a1M
    a21                                                                   0
                0     0    a23  .    .    . a2M
    a31                               a32                                 0
                            0                    .
L                                   =                                      .
       U=                                          .
    .
                                                       .
    .
                                                 aM-1M
    aM1       aM2          .           .          .         aMM-1         0
        0                                  0


И матрица D - диагональная.
                           (k)     (k)                   (k)
    Обозначим через Yk = ( Y1 ,Y2  ...  YM  )  вектор  k-ого  итерационного
шага. Пользуясь этими обозначениями запишем метод Зейделя иначе :

                        (   D   +    L    )Yk+1    +     UYk    =    f    ,
k=0,1...

Приведём эту итерационную схему к каноническому виду двухслойных схем :

                       (   D   +   L   )(Yk+1   -   Yk)   +AYk   =   f    ,
k=0,1...

    Мы рассмотрели так называемый точечный  или  скалярный  метод  Зейделя,
анологично строится блочный или векторный метод  Зейделя  для  случая  когда
aii - есть квадратные матрицы, вообще говоря, различной размерности,  а  aij
для i<>j - прямоугольные матрицы. В  этом  случае  Yi  и  fi  есть  векторы,
размерность которых соответствует размерности матрицы aii.



                        ПОСТРОЕНИЕ  РАЗНОСТНЫХ  СХЕМ


    Пусть  Yi=Y(i)  сеточная  функция  дискретного  аргумента  i.  Значения
сеточной функции Y(i)  в свою  очередь  образуют  дискретное  множество.  На
этом множестве можно определять  сеточную  функцию,  приравнивая  которую  к
нулю получаем  уравнение  относительно  сеточной  функции  Y(i)  -  сеточное
уравнение.  Специальным  случаем  сеточного  уравнения  является  разностное
уравнение.
    Сеточное уравнение получается при аппроксимации на сетке интегральных и
дифференциальных уравнений.
    Так дифференциальное уравнение первого порядка :

                          dU  = f(x)   ,      x > 0
                          dx

можно заменить разностным уравнением первого порядка :

                 Yi+1 - Yi = f(xi) ,     xi = ih,  i=0,1...
                        h

или Yi+1=Yi+hf(x),   где  h  -  шаг  сетки  v={xi=ih,  i=0,1,2...}.  Искомой
функцией является сеточная функция Yi=Y(i).
    При разностной аппроксимации уравнения второго поряда

                               2
                                d U    = f(x)
                                 2
                              dx

получим разностное уравнение второго порядка :

                       2
                  Yi+1 - 2Yi + Yi+1 = yi ,    где  yi=h f i
                                           fi  = f(xi)
                                           xi  = ih

Для разностной aппроксимации  производных U’, U’’, U’’’  можно  пользоваться
шаблонами с большим числом  узлов.  Это  приводит  к  разностным  уравнениям
более высокого порядка.
    Анологично  определяется  разностное  уравнение  относительно  сеточной
функции Uij = U(i,j)  двух   дискретных  аргументов . Например  пятиточечная
разностная схема “крест” для уравнения Пуассона

                             Uxx + Uyy = f(x,y)

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

               Ui-1j - 2Uij+Ui+1j  +  Uij-1 - 2Uij+Uij+1 = fij
                        2
  2
                       hx                             hy

где  hx - шаг сетки по X
    hy -  шаг сетки по Y
Сеточное уравнение общего вида можно записать так:

                      N
                     CijUj = fi               i=0,1...N
                      j=0

Оно содержит все  значения  U0,  U1  ...  UN  сеточной  функции.  Его  можно
трактовать как рзностное уравнение  порядка  N  равного  числу  узлов  сетки
минус единица.
    В общем случае под  i  -  можно  понимать  не  только  индекс  ,  но  и
мультииндекс т.е. вектор i = (i1 ... ip)  с  целочисленными  компонентами  и
тогда :


                            СijUj =fi      i О W
                           jОW

где сумирование происходит по всем узлам сетки W. Если коэффициенты  Сij  не
зависят от i, тоуравнение называют уравнением с постоянными коэффициентами.
    Аппроксимируем нашу задачу т.е. заменим уравнение и краевые условия  на
соответствующие им сеточные уравнения.

    U=U(x,y)

          y

        M  b

       M-1


                      Uij                                         j
    j



    1
      0 1 2                                  i
    N-1   N=a     x

                    i
Построим  на  области  G  сетку  W  .  И  зададим  на  W  сеточную  функцию
Uij=U(xi,yj) ,
 где
xi=x0+ihx
yi=y0+jhy
hx = a/N ,
hy = b/M  и т.к.
x0=y0
то
xi=ihx,  yi=jhy,  i=0...N
                         j=0...M

Найдём разностные производные входящие в уравнение
         2
                                   DU = f

(т.е построим разностный аналог бигармонического уравнения).


         Uxij  =  Ui+1j - Uij    ,             Uxi-1j = Uij - Ui-1j
                           hx
                  hx

                       Uxxij  =   Ui-1j - 2Uij + Ui+1j
                                   hx

Рассмотрим Uxxxxij как разность третьих производных  :

                    Uxxi-1j - Uxxij - Uxxij - Uxxi+1j
 Uxxxxij =              hx                   hx          =    Ui-2j - 4Ui-1j
                           + 6Uij - 4Ui+1j + Ui+2j
                                                    4
                                                                          hx
                              hx
Анологично вычислим производную по y :

               Uyyyyij = Uij-2 - 4Uij-1 + 6Uij - 4Uij+1 +Uij+2
                                        4
                                     hy

Вычислим смешанную разностную производную Uxxyy :

                       Uxxij-1 - Uxxij   -   Uxxij - Uxxij+1
   (Uxx)yyij =           hy                              hy              =
                        Uxxij-1 - 2Uxxij +Uxxij+1   =
                                                   2
                                                                          hy
                          hy

 =   Ui-1j-1 - 2Uij-1 + Ui+1j-1   -   2 Ui-1j - 2Uij + Ui+1j    +    Ui-1j-1
                             - 2Uij+1 + Ui+1j+1
              2                                                           2
      2   2                                                        2   2
             hxhy                                                      hxhy
             hxhy

В силу того что   DU = f
имеем:

                  Ui-2j - 4Ui-1j + 6Uij - 4Ui+1j +Ui+2j  +
                                        4
                                      hx

   + 2 Ui-1j-1 - 2Uij-1 + Ui+1j-1  -  4 Ui-1j - 2Uij +Ui+1j   +  2 Ui-1j+1
                            -2Uij+1 + Ui+1j+1   +
              2  2                                                        2
 2                                                              2   2
             hxhy                                                      hxhy
                                                        hxhy

            +  Uij-2 - 4Uij-1 + 6Uij - 4Uij+1 + Uij+2     =   fij
                                             (*)
                     4
                    hy

Это уравнение имеет место для
                               i=1,2, ... N-1
                               j=1,2, ... M-1
Рассмотрим краевые условия задачи. Очевидно следующее :
                              x=0   ~    i = 0
                               x=a   ~   xN=a
                               y=0   ~   Yo=0
                               y=b   ~   YM=b

1)   х=0    (левая граница области G)

      Заменим условия
                                 U    = 0
                                   x=o
                                 Uxxx  = 0
                                          x=o

на соответствующие им разностные условия

                  Uoj=0
        U-1j=U2j - 3U1j                                         (1`)


2)  х=а      (правая граница области G)
     i=N

             Ux   = 0
                   x=a
             Uxxx  = 0
                       x=a                    из того что    Ui+1j  -  Ui-1j
= 0

2hx

                                      UN+1j = UN-1j
          UNj = 4 UN-1j - UN-2j                               (2`)
                                                           3

3)  у=0 (нижняя граница области G)
     j=0

                  Ui ,-1 = Ui1
        Ui0 = 0                                                  (3`)

это есть разностный аналог Uy   = 0
                                                      y=o
                                                U   =0
                                                    y=o

4)   у=b
      i=M

            U  = 0
                          y=b                          т.е.            UiM=0
               (**)

Распишем через разностные производные  Uxx + Uyy =0 и  учитывая  что  j=M  и
(**) получим

                                UiM-1 = UiM+1

Итак краевые условия на у=b имеют вид

                   UiM+1 = UiM-1
        UiM = 0                                                 (4`)

Итого  наша  задача  в  разностных  производных  состоит  из  уравнения  (*)
заданного на сетке  W  и  краевых  условий  (1`)-(4`)  заданных  на  границе
области G (или на границе сетки W)



                          ПРИМЕНЕНИЕ МЕТОДА ЗЕЙДЕЛЯ


Рассмотрим применение метода Зейделя для  нахождения  приближенного  решения
нашей  разностной задачи (*),(1`) - (4`).
В данном случае неизвестными являются

                               Uij = U(xi,yj)
где  xi = ihx
       yj = jhy
при чём  hx = a/N  ,
                hy = b/M
это есть шаг сетки по x и по у соответственно ,  а  N   и  М  соответственно
количество точек  разбиения отрезков [0 , а] и [0 , b]
Пользуясь результатами предыдущего раздела запишем уравнение

                                2
                                   DU = f

как разностное уравнение. И упорядочим неизвестные естественным  образом  по
строкам сетки W , начиная с нижней строки.

    1 Ui-2j -  4  +  4   Ui-1j  +  6  -  8 + 6   Uij  -  4  +   4   Ui+1j +
    1 Ui+2j +  2Ui-1j-1 -
      4            4          2   2                4         2    2       4
                      4               2      2                            4
       2  2
    hx         hx   hxhy             hx   hxhy  hy             hx      hxhy
             hx         hxhy


    -    4  +  4  Uij-1 + 2 Ui+1j-1 + 2 Ui-1j+1  -   4   +   4   Uij+1 +  2
    Ui+1j+1 +  1 Uij-2 +
            2  2          4                    2    2                     2
    2                          2  2            4                          2
     2                        4
        hxhy     hy           hxhy           hxhy              hxhy      hy
          hxhy             hy


    + 1 Uij+2  =  f ij      для   i=1 ... N-1, j=1 ... M-1
        4
           hy
и U удовлетворяет краевым условиям (1`) - (4`), так как  в каждом  уравнении
связаны вместе не более 13 неизвестных то в матрице А  отличны  от  нуля  не
более 13-элементов в строке. В соответствии со  вторым  разделом   перепишем
уравнение:

                           (k+1)                                       (k+1)
                                                                       (k+1)
(k+1)
     6   -    8   +   6   Uij      =       -  1 Uij-2     -    2    Ui-1j-1
       +      4   +    4     Uij-1   -
               4                     2             2                      4
                     4                                     2              2
    2   2              4
    hx      hxhy     hy                         hy                     hxhy
                       hxhy       hy


                                    (k+1)                              (k+1)
                                                                       (k+1)
                                        (k)
    -    2   Ui+1j-1    -     1   Ui-1j     +      4     +     4      Ui-1j
    +          4    +   4    Ui+1j  -
        2   2                   4                    4                 2  2
                                     4              2  2
         hxhy                    hx                      hx            hxhy
                hx       hxhy

                               (k)                                      (k)
                                           (k)                          (k)
                            (k)
    -   1  Ui+2j   -    2  Ui-1j+1    +     4     +   4     Uij+1    -    2
    Ui+1j+1  -   1  Uij+2     +  fij
                     4                                      2             2
         2     2                  4                                2      2
               4
           hx                   hxhy                     hxhy            hy
    hxhy               hy



                               (k)
При  чем   U   удовлетворяет  краевым  условиям  (1`)  -  (4`).   Вычисления
начинаются с i=1,  j=1 и продолжаются либо  по  строкам   либо  по  столбцам
сетки  W. Число неизвестных в задаче n = (N-1)(M-1).
    Как  видно  из  вышеизложенных  рассуждений   шаблон  в   этой   задаче
тринадцатиточечный  т.е. на каждом шаге в разностном уравнении участвуют  13
точек (узлов сетки) Рассмотрим вид матрицы А - для данной задачи.

j+2
j+1
j
j-1
Матрица  метода   получается   следующим   образом   :   все    узлы   сетки
перенумеровываются и размещаются в  матрице  Так что  все   узлы    попадают
на  одну строку    и  поэтому   матрица  метода   для   нашей  задачи  будет
тринадцатидиагональной .

j-2



i-1
i
i+1
i+2

i-2

Шаблон задачи



                             ОПИСАНИЕ ПРОГРАММЫ.


    Константы используемые в программе :

         aq = 1 - правая граница области G
         b   = 1 - левая граница области G
         N  = 8 - колличество точек разбиения отрезка [0,a]
         M = 8 - колличество точек разбиения отрезка  [0,b]
         h1 = aq/N - шаг сетки по X
         h2 = b/M - шаг сетки по Y

    Переменные :

         u0 - значения сеточной функции U на k-ом шаге
         u1 - значения сеточной функции U на (k+1)-ом шаге
         a - массив коэффициентов шаблона

    Описание процедур :

         procedure Prt(u:masa) - печать результата

         function ff(x1,x2: real):real - возвращает значение  функции  f   в
                                 узле (x1,x2)
         procedure Koef - задаёт значения коэффициентов

    Действие :

         Берётся начальое приближение u0 и с учётом краевых условий  ведётся
         вычисление с i=2 ... N , j=2 ... M.  На  каждом  итерационном  шаге
         получаем u1 по u0. По достижении заданной точности eps>0 вычисления
         прекращаются. И все элементы матрицы A, которые лежат ниже  главной
         диагонали получают итерационный шаг (k+1) , а те  элементы  которые
         лежат выше главной диагонали (исключая главную диагональ)  получают
         итерационный шаг k.

    Примечание : программа реализована на языке Borland Pascal 7.0



           Министерство общего и профессионального образования РФ
                   Воронежский государственный университет



                                факультет ПММ
                     кафедра Дифференциальных уравнении



                               Курсовой проект
            “Решение бигармонического уравнения методом Зейделя”



                                      Исполнитель : студент 4 курса 5 группы
                                                                Никулин Л.А.

                                       Руководитель : старший  преподаватель
                                                                 Рыжков А.В.



                               Воронеж 1997г.





смотреть на рефераты похожие на "Метод конечных разностей или метод сеток "