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

ВВЕДЕНИЕ
Значительнаое число задач физики и техники приводят к дифференциальным уравнениям в частных прозводных (уравнения математической физики). Установившиеся процессы различной физической природы описываются уравнениями эллиптического типа.
Точные решения краевых задач для эллиптических уравнений удаётся получить лишь в частных случаях. Поэтому эти задачи решают в основном приближённо. Одним из наиболее универсальных и эффективных методов, получивших в настоящее время широкое распространение для приближённого решения уравнений математической физики, является метод конечных разностей или метод сеток.
Суть метода состоит в следующем. Область непрерывного изменения аргументов, заменяется дискретным множеством точек (узлов), которое называется сеткой или решёткой. Вместо функции непрерывного аргумента рассматриваются функции дискретного аргумента, определённые в узлах сетки и называемые сеточными функциями. Производные, входящие в дифференциальное уравнение и граничные условия, заменяются разностными производными, при этом краевая задача для дифференциального уравнения заменяется системой линейных или нелинейных алгебраических уравнений (сеточных или разностных уравнений). Такие системы часто называют разностными схемами. И эти схемы решаются относительно неизвестной сеточной функции.
Далее мы будем рассматривать применение итерационного метода Зейделя для вычисления неизвестной сеточной функции в краевой задаче с неоднородным бигармоническим уравнением. ПОСТАНОВКА ЗАДАЧИ Пусть у нас есть бигармоническое уравнение : 2 U = f
Заданное на области G={ (x, y) : 0 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) отличны от нуля (aii0) записывается в следующем виде : 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) Так как a110 то отсюда найдём 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 0a12 a13 ... . a1M a210 0 0a23 ... . a2M a31 a320 0 . L = . U= . ... . aM-1M aM1 aM2 ... . aMM-10 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 для ij - прямоугольные матрицы. В этом случае 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г.