ВВЕДЕНИЕ
Значительнаое число задач физики и техники приводят к дифференциальным уравнениям в частных прозводных (уравнения математической физики). Установившиеся процессы различной физической природы описываются уравнениями эллиптического типа.
Точные решения краевых задач для эллиптических уравнений удаётся получить лишь в частных случаях. Поэтому эти задачи решают в основном приближённо. Одним из наиболее универсальных и эффективных методов, получивших в настоящее время широкое распространение для приближённого решения уравнений математической физики, является метод конечных разностей или метод сеток.
Суть метода состоит в следующем. Область непрерывного изменения аргументов, заменяется дискретным множеством точек (узлов), которое называется сеткой или решёткой. Вместо функции непрерывного аргумента рассматриваются функции дискретного аргумента, определённые в узлах сетки и называемые сеточными функциями. Производные, входящие в дифференциальное уравнение и граничные условия, заменяются разностными производными, при этом краевая задача для дифференциального уравнения заменяется системой линейных или нелинейных алгебраических уравнений (сеточных или разностных уравнений). Такие системы часто называют разностными схемами. И эти схемы решаются относительно неизвестной сеточной функции.
Далее мы будем рассматривать применение итерационного метода Зейделя для вычисления неизвестной сеточной функции в краевой задаче с неоднородным бигармоническим уравнением.
ПОСТАНОВКА ЗАДАЧИ
Пусть у нас есть бигармоническое уравнение :
 2
U
=
f
Заданное на области G={
(x,y)
: 0<=x<=a,
0<=y<=b
}.
Пусть также заданы краевые условия на границе области G
.
  U
= 0 Y
x=0 b
  U
xxx = 0
x=0
G
 U
x = 0
x=a
 U
xxx = 0 0 a X
x=a
  U
= 0 U
= 0
y=0 y=b
  U
y = 0 U
xx + U
yy = 0
y=0 y=b y=b
Надо решить эту задачу численно.
  Для решения будем использовать итерационный метод Зейделя для решения сеточных задач.
 По нашей области G
построим равномерные сетки W
x
и W
y
с шагами h
x
и h
y
соответственно .
 W
x
={
x(i)=i
h
x
,
i=0,1...N,
h
x
N=a
}
 W
y
={
y(j)=j
h
y
,
j=0,1...M,
h
y
M=b
}
Множество узлов U
ij
=(x(i),y(j))
имеющих координаты на плоскости х(i)
,y(j)
называется сеткой в прямоугольнике G
и обозначается :
W={
U
ij
=(ih
x
,j
h
y
),
i=0,1...N,
j=0,1...M,
h
x
N=a,
h
y
M=b
}
Сетка W
очевидно состоит из точек пересечения прямых x=x(i)
и y=y(j)
.
  Пусть задана сетка W
.Множество всех сеточных функций заданных на W
образует векторное пространство с определённом на нём сложениемфункций и умножением функции на число. На пространстве сеточных функций можно определитьразностные или сеточные операторы. 0ператор A
преобразующий сеточную функцию U
в сеточную функцию f=AU
называется разностным или сеточным оператором. Множество узлов сетки используемое при написании разностного оператора в узле сетки называется шаблоном этого оператора.
Простейшим разностным оператором является оператор дифференцирования сеточной функции, который порождает разностные производные. Пусть W
- сетка с шагом h
введённая на R
т.е.
W={X
i
=a+ih, i=0, +
1, +
2...}
Тогда разностные производные первого порядка для сеточной функции Y
i
=Y(X
i
)
, X
i
из W
, определяется по формулам :
L
1
Y
i
= Y
i
- Y
i-1
,
L
2
Y
i
=
L
1
Y
i+1
h
и называются соответственно левой и правой производной. Используется так же центральная производная :
L
3
Y
i
=Y
i+1
- Y
i-1
= (
L
1
+
L
2
)Y
i
2h 2
Разностные операторы A
1
, A
2
, A
3
имеют шаблоны состоящие 2х точек и используются при апроксимации первой производной Lu=u’
. Разностные производные n
-ого порядка определяются как сеточные функции получаемые путём вычисления первой разностной производной от функции, являющейся разностной производной n-1
порядка, например :
Y
xxi
=Y
xi+1
- Y
xi
= Y
i-1
-2Y
i
+Y
i+1
2
h h
Y
xxi
= Y
xi+1
-Y
xi-1
= Y
i-2
- 2Y
i
+Y
i+ 2
2
2h 4h
которые используются при апроксимации второй производной. Соответствующие разностные операторы имеют 3х точечный шаблон.
Анологично не представляет труда определить разностные производные от сеточных функций нескольких переменных.
Аппроксомируем нашу задачу с помощью разностных производных. И применим к получившейся сеточной задаче метод Зейделя.
МЕТОД ЗЕЙДЕЛЯ
Одним из способов решения сеточных уравнений является итерационный метод Зейделя.
Пусть нам дана система линейных уравнений :
AU
= f
или в развёрнутом виде :
M
 a
ij
U
j
= f
i
, i=1,2...M
i=1
Итерационный метод Зейделя в предположении что диагональные элементы матрицы А=(
a
ij
)
отличны от нуля (
a
ii
<>0
)
записывается в следующем виде :
 i (k+1)
M
(k)
a
ij
Y
j
+ a
ij
Y
j
= f
i
, i=1,2...M
j=1 j=i+1
(k)
где Y
j
- j
ая компонента итерационного приближения номера k
. В качестве начального приближения выбирается произвольный вектор.
Определение (
k
+1)
-ой итерации начинается с i=1
(k+1) M (k)
 a
11
Y
1
= - a
1j
Y
j
+f
1
j=2
(k+1)
Так как a
11
<>0
то отсюда найдём Y
1
. И для i=2
получим :
 (
k+1
)
(k+1)
M (k)
a
22
Y
2
= - a
21
Y
1
- a
2j
Y
j
+ f
2
j=3
(k+1) (k+1) (k+1)
(
k+1
)
Пусть уже найдены Y
1
, Y
2
... Y
i-1
. Тогда Y
i
находится из уравнения :
  (k+1) i-1 (k+1) M (k)
a
ii
Y
i
= - a
ij
Y
j
- a
ij
Y
j
+ f
i
(*)
j=1 j=i+1
Из формулы (*)
видно , что алгоритм метода Зейделя черезвычайно прост. Найденное по формуле (*)
значение Y
i
размещается на месте Y
i
.
Оценим число арифметических действий, которое требуется для реализации одного итерационного шага. Если все a
ij
не равны нулю, то вычисления по формуле (*)
требуютM-1
операций умножения и одного деления. Поэтому реализация
2
одного шага осуществляется за 2M - M
арифметических действий.
Если отлично от нуля лишь m элементов, а именно эта ситуация имеет место для сеточных эллиптических уравнений, то на реализацию итерационного шага потребуется 2
Mm-M
действий т.е. число действий пропорционально числу неизвестных M
.
Запишем теперь метод Зейделя в матричной форме. Для этого представим матрицу A
в виде суммы диагональной, нижней треугольной и верхней треугольной матриц :
A = D + L + U
где
0 0 . . . 0 0
a
12
a
13
. . .
a
1M
a
21
0 0 0
a
23
. . .
a
2M
a
31
a
32
0 0 .
L = . U= .
. .
.
a
M-1M
a
M1
a
M2
. . .
a
MM-1
0 0 0
И матрица D
- диагональная.
(k) (k) (k)
Обозначим через Y
k
= ( Y
1
,Y
2
... Y
M
)
вектор k
-ого итерационного шага. Пользуясь этими обозначениями запишем метод Зейделя иначе :
(
D + L
)
Y
k+1
+ UY
k
= f
, k=0,1...
Приведём эту итерационную схему к каноническому виду двухслойных схем :
(
D + L
)
(Y
k+1
- Y
k
) +AY
k
= f , k=0,1...
Мы рассмотрели так называемый точечный или скалярный метод Зейделя, анологично строится блочный или векторный метод Зейделя для случая когда a
ii
- есть квадратные матрицы, вообще говоря, различной размерности, а a
ij
для i<>j
- прямоугольные матрицы. В этом случае Y
i
и f
i
есть векторы, размерность которых соответствует размерности матрицы a
ii
.
ПОСТРОЕНИЕ РАЗНОСТНЫХ СХЕМ
Пусть Y
i
=Y(i)
сеточная функция дискретного аргумента i
. Значения сеточной функции Y(i)
в свою очередь образуют дискретное множество. На этом множестве можно определять сеточную функцию, приравнивая которую к нулю получаем уравнение относительно сеточной функции Y(i)
- сеточное уравнение. Специальным случаем сеточного уравнения является разностное уравнение.
Сеточное уравнение получается при аппроксимации на сетке интегральных и дифференциальных уравнений.
Так дифференциальное уравнение первого порядка :
dU
= f
(x) , x > 0
dx
можно заменить разностным уравнением первого порядка :
Y
i+1
- Y
i
= f
(x
i
) , x
i
= ih, i=0,1...
h
илиY
i+1
=Y
i
+hf
(x)
, где h
- шаг сетки v
={x
i
=ih, i=0,1,2...}
. Искомой функцией является сеточная функция Yi=Y(i)
.
При разностной аппроксимации уравнения второго поряда
2
d U
= f
(x)
 2
dx
получим разностное уравнение второго порядка :
2
Y
i+1
- 2Y
i
+ Y
i+1
=
y
i
,
где
y
i
=h
f
i
f
i
= f(x
i
)
x
i
= ih
Для разностной aппроксимациипроизводных U’,
U’’,
U’’’
можно пользоваться шаблонами с большим числом узлов. Это приводит к разностным уравнениям более высокого порядка.
Анологично определяется разностное уравнение относительно сеточной функции U
ij
= U(i,j)
двух дискретных аргументов. Например пятиточечная разностная схема “крест” для уравнения Пуассона
U
xx
+ U
yy
= f
(x,y)
на сетке W
выглядит следующим образом :
U
i-1j
-
2U
ij
+U
i+1j
+
U
ij-1
-
2U
ij
+U
ij+1
=
f
ij
2 2
h
x
h
y
где h
x
- шаг сетки по X
h
y
- шаг сетки поY
Сеточное уравнение общего вида можно записать так:
N
C
ij
U
j
= f
i
i=0,1...N
j=0
Оно содержит все значения U
0
, U
1
... U
N
сеточной функции. Его можно трактовать как рзностное уравнение порядка N
равного числу узлов сетки минус единица.
В общем случае под i
- можно понимать не только индекс , но и мультииндекс т.е. вектор i = (i
1
... i
p
)
с целочисленными компонентами и тогда :
 С
ij
U
j
=f
i
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
сеточную функцию U
ij
=U(x
i
,y
j
)
,
где
x
i
=x
0
+ih
x
y
i
=y
0
+jh
y
h
x
= a/N ,
h
y
= b/M
и т.к.
x
0
=y
0
то
x
i
=ih
x
, y
i
=jh
y
, i=0...N
j=0...M
Найдём разностные производные входящие в уравнение
2
D
U = f
(т.е построим разностный аналог бигармонического уравнения).
Ux
ij
= U
i+1j
- U
ij
, Ux
i-1j
= U
ij
- U
i-1j
h
x
h
x
Uxx
ij
= U
i-1j
- 2U
ij
+ U
i+1j
h
x
Рассмотрим Uxxxx
ij
как разность третьих производных :
Uxx
i-1j
- Uxx
ij
- Uxx
ij
- Uxx
i+1j
Uxxxx
ij
= h
x
h
x
= U
i-2j
- 4U
i-1j
+ 6U
ij
- 4U
i+1j
+ U
i+2j
4
h
x
h
x
Анологично вычислим производную по y
:
Uyyyy
ij
= U
ij-2
- 4U
ij-1
+ 6U
ij
- 4U
ij+1
+U
ij+2
4
h
y
Вычислим смешанную разностную производнуюUxxyy
:
Uxx
ij-1
- Uxx
ij
- Uxx
ij
- Uxx
ij+1
(Uxx)yy
ij
= h
y
h
y
= Uxx
ij-1
- 2Uxx
ij
+Uxx
ij+1
=
2
hy hy
= U
i-1j-1
- 2U
ij-1
+ U
i+1j-1
-
2
U
i-1j
- 2U
ij
+ U
i+1j
+ U
i-1j-1
- 2U
ij+1
+ U
i+1j+1
2 2 2 2 2 2
h
x
h
y
h
x
h
y
h
x
h
y
В силу того чтоD
U =
f
имеем:
U
i-2j
- 4U
i-1j
+ 6U
ij
- 4U
i+1j
+U
i+2j
+
4
h
x
+
2
U
i-1j-1
- 2U
ij-1
+ U
i+1j-1
-
4
U
i-1j
- 2U
ij
+U
i+1j
+
2
U
i-1j+1
-2U
ij+1
+ U
i+1j+1
+
2 2 2 2 2 2
h
x
h
y
h
x
h
y
h
x
h
y
+ U
ij-2
- 4U
ij-1
+ 6U
ij
- 4U
ij+1
+ U
ij+2
=
f
ij
(*)
4
h
y
Это уравнение имеет место для
i=1,
2, ...
N-1
j=1,2, ... M-1
Рассмотрим краевые условия задачи. Очевидно следующее :
x=
0
~
i
=
0
x=a ~ x
N
=a
y=0 ~ Yo=0
y=b ~ Y
M
=b
1)
х=0
(левая граница областиG
)
Заменим условия
U = 0
x=o
Uxxx = 0
x=o
на соответствующие им разностные условия
U
o
j
=0
U
-1j
=U
2j
- 3U
1j
(1`)
2)
х=а
(правая граница областиG
)
i=N
 Ux = 0
 x=a
Uxxx = 0
x=a
из того что U
i+1j
- U
i-1j
= 0
2h
x
U
N+1j
= U
N-1j
U
Nj
= 4 U
N-1j
- U
N-2j
(2`)
3
3)
у=0
(нижняя граница области G
)
  j=0
U
i ,-1
= U
i1
U
i0
= 0 (3`)
это есть разностный аналогUy = 0
y=o
U =0
y=o
4)
у=b
i=M
U = 0
y=b
т.е.U
iM
=0
(**)
Распишем через разностные производныеUxx + Uyy =0
и учитывая чтоj=M
и (**)
получим
U
iM-1
= U
iM+1
Итак краевые условия на у=b
имеют вид
U
iM+1
= U
iM-1
U
iM
= 0 (4`)
Итого наша задача в разностных производных состоит из уравнения (*)
заданного на сетке W
и краевых условий (1
`
)-(4
`
)
заданных на границе области G
(или на границе сетки W
)
ПРИМЕНЕНИЕ МЕТОДА ЗЕЙДЕЛЯ
Рассмотрим применение метода Зейделя для нахождения приближенного решения нашей разностной задачи (*)
,(1`) - (4`).
В данном случае неизвестными являются
U
ij
= U(
x
i
,y
j
)
где x
i
= ih
x
y
j
= jh
y
при чём h
x
= a/N ,
h
y
= b/M
это есть шаг сетки по x
и по у
соответственно , а N
и М
соответственно количество точек разбиения отрезков [
0
,
а]
и [0
,
b]
Пользуясь результатами предыдущего раздела запишем уравнение
2
D
U = f
как разностное уравнение. И упорядочим неизвестные естественным образом по строкам сетки W
, начиная с нижней строки.
   1
U
i-2j
- 4
+ 4
U
i-1j
+ 6
- 8
+ 6
U
ij
- 4
+ 4
U
i+1j
+ 1
U
i+2j
+ 2
U
i-1j-1
-
  4
4 2 2
4 2 2 4 4 2 2 4 2 2
h
x
h
x
h
x
h
y
h
x
h
x
h
y
h
y
h
x
h
x
h
y
h
x
h
x
h
y
 - 4
+ 4
U
ij-1
+ 2
U
i+1j-1
+ 2
U
i-1j+1
- 4
+ 4
U
ij+1
+ 2
U
i+1j+1
+ 1
U
ij-2
+
 2 2 4 2 2 2 2 2 2 4 2 2 4
 h
x
h
y
h
y
h
x
h
y
h
x
h
y
h
x
h
y
h
y
h
x
h
y
h
y
+ 1
U
ij+2
=
f
ij
для i=1 ... N-1, j=1 ... M-1
4
h
y
и U
удовлетворяет краевым условиям (1
`
)
- (4`),
так как в каждом уравнении связаны вместе не более 13 неизвестных то в матрице А
отличны от нуля не более 13-элементов в строке. В соответствии со вторым разделом перепишем уравнение:
 (k+1) (k+1) (k+1) (k+1)
6
- 8
+ 6
U
ij
= - 1
U
ij-2
- 2
U
i-1j-1
+ 4
+ 4
U
ij-1
-
4 2 2 4 4 2 2 2 2 4
h
x
h
x
h
y
h
y
h
y
h
x
h
y
h
x
h
y
h
y
   (k+1) (k+1) (k+1) (k)
- 2
U
i+1j-1
- 1
U
i-1j
+
 4
+ 4
U
i-1j
+ 4
+ 4
U
i+1j
-
2 2
4
4 2 2
4 2 2
 h
x
h
y
h
x
h
x
h
x
h
y
h
x
h
x
h
y
(k) (k) (k) (k) (k)
- 1
U
i+2j
- 2
U
i-1j+1
+ 4
+ 4
U
ij+1
- 2
U
i+1j+1
- 1
U
ij+2
+
f
ij
4 2 2 2 2 4 2 2 4
  h
x
h
x
h
y
h
x
h
y
h
y
h
x
h
y
h
y
(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г.
|
