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