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

ВВЕДЕНИЕ





Значительнаое число задач физики и техники приводят к дифференциальным уравнениям в частных прозводных (уравнения математической физики). Установившиеся процессы различной физической природы описываются уравнениями эллиптического типа.

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

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

Далее мы будем рассматривать применение итерационного метода Зейделя для вычисления неизвестной сеточной функции в краевой задаче с неоднородным бигармоническим уравнением.









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





Пусть у нас есть бигармоническое уравнение :

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, определяется по формулам :

	

	1Yi  = Yi - Yi-1    , 2Yi=1Yi+1

	                                                               h

	

	и называются соответственно левой и правой производной. Используется так же центральная производная :

	

	 3Yi=Yi+1 - Yi-1  = (1+2)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 - шаг сетки ={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

		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

		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г.