ВВЕДЕНИЕ

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

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

 
Пусть у нас есть бигармоническое уравнение :
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)=ih x , i=0,1...N, h x N=a }

---

W y ={ y(j)=jh y , j=0,1...M, h y M=b }
Множество узлов U ij =(x(i),y(j)) имеющих координаты на плоскости х(i) , y(j) называется сеткой в прямоугольнике G и обозначается :

---

---

W={ U ij =(ih x ,jh 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=A U называется разностным или сеточным оператором. Множество узлов сетки используемое при написании разностного оператора в узле сетки называется шаблоном этого оператора.
Простейшим разностным оператором является оператор дифференцирования сеточной функции, который порождает разностные производные. Пусть 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х точечный шаблон.
Анологично не представляет труда определить разностные производные от сеточных функций нескольких переменных.
Аппроксомируем нашу задачу с помощью разностных производных. И применим к получившейся сеточной задаче метод Зейделя.

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

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

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

M
a ijU 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 элементов, а именно эта ситуация имеет место для сеточных эллиптических уравнений, то на реализацию итерационного шага потребуется 2Mm-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 +h f (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)

---

---

Построим на области 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 oj =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г.