Метод конечных разностей или метод сеток
ВВЕДЕНИЕ
Значительнаое
число задач физики и техники приводят к дифференциальным уравнениям в частных
прозводных (уравнения математической физики). Установившиеся процессы различной
физической природы описываются уравнениями эллиптического типа.
Точные
решения краевых задач для эллиптических уравнений удаётся получить лишь в
частных случаях. Поэтому эти задачи решают в основном приближённо. Одним из
наиболее универсальных и эффективных методов, получивших в настоящее время
широкое распространение для приближённого решения уравнений математической
физики, является метод конечных разностей или метод сеток.
Суть
метода состоит в следующем. Область непрерывного изменения аргументов,
заменяется дискретным множеством точек (узлов), которое называется сеткой или решёткой.
Вместо функции непрерывного аргумента рассматриваются функции дискретного
аргумента, определённые в узлах сетки и называемые сеточными функциями.
Производные, входящие в дифференциальное уравнение и граничные условия,
заменяются разностными производными, при этом краевая задача для
дифференциального уравнения заменяется системой линейных или нелинейных
алгебраических уравнений (сеточных или разностных уравнений). Такие системы
часто называют разностными схемами. И эти схемы решаются относительно неизвестной
сеточной функции.
Далее
мы будем рассматривать применение итерационного метода Зейделя для вычисления
неизвестной сеточной функции в краевой задаче с неоднородным бигармоническим
уравнением.
ПОСТАНОВКА
ЗАДАЧИ
Пусть у нас есть
бигармоническое уравнение :
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
(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
Аппроксимируем
нашу задачу т.е. заменим уравнение и краевые условия на соответствующие им
сеточные уравнения.
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
- 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 точек (узлов сетки)
Рассмотрим вид матрицы А - для данной задачи.
Матрица метода
получается следующим образом : все узлы сетки перенумеровываются и
размещаются в
матрице
Так что
все
узлы попадают на одну
строку и поэтому
матрица метода для нашей
задачи будет тринадцатидиагональной .
ОПИСАНИЕ
ПРОГРАММЫ.
Константы используемые в программе :
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г.