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

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 - шаг сетки ={xi=ih, i=0,1,2...}. Искомой функцией является сеточная функция Yi=Y(i).

При разностной аппроксимации уравнения второго поряда


2

d U = f(x)

2

dx


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


2

Yi+1 - 2Yi + Yi+1 = i , где i=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

jW


где сумирование происходит по всем узлам сетки 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

U = 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


В силу того что U = 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]

Пользуясь результатами предыдущего раздела запишем уравнение


U = 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г.