Автоматихация производства ЭВА, страница 3

Метод решения задачи, записанной в виде ДУ, с помощью разностного уравнения называют методом конечных разностей. При таком подходе решение ДУ заменяется решением системы линейных алгебраических уравнений с количеством неизвестных, равных количеству дискрет разбиения области определения функции F.

Например, рассмотрим уравнение Пуассона: f (x,y)

д2F	+	д2F	= f(x,y)
дХ2		дY2

Переходя от вторых производных к конечным разностям в точке (mh, nh) области Р получим разностное уравнение вида:

Fm+1n-2Fmn + Fm-1,n	+	Fm,n+1-2Fmn + Fm,n-1	= f(x,y)	(12)
h2		h2

Формируя уравнение (12) для всех точек области Р, получим систему алгебраических уравнений с числом неизвестных равным числу дискрет области:

f₂₁ + f₁₂ – 4f₁₁ + f₀₁ + f₁₀ = f(1,1)

f₃₁ + f₂₂ – 4f₂₁ + f₂₁ + f₂₁ = f(2,1) (13)

  

f_m+1,n + f_m,n+1 – 4f_mn + f_m-1,n + f_m,n-1 = f(m,n)

При заданных на границе области Р значениях функции f данное уравнение может иметь единственное решение, которое и определит дискретную модель непрерывной величины f в области Р.

Система (13) позволяет определить приближенное значение функции F в области Р. Необходимо отметить, что количество уравнений может быть весьма велико – несколько сотен и более. Решать подобные уравнения без помощи ЭВМ не возможно. Для решения подобных систем линейных уравнений на ЭВМ используются известные методы Эйлера и Гаусса.

Отметим также, что свести решение ДУ в частных производных к решению систем алгебраических уравнений удается не всегда, а только в случае стационарных процессов (установившихся во времени). При моделировании нестационарных процессов в ДУ появляются члены, зависящие от времени. Методы решения таких задач будут рассмотрены далее (раздел ?).

Рассмотрим построение разностной схемы на примере следующего уравнения:

д2F	+	д2F	+ e - x	дF	- F = f(x,y)	(14)
дХ2		дY2		дX

в квадрате Р = {0х1; 0 y1} с краевым условием:

(^x, ^y) = exp(^x-^y) (15)

То есть функция F в точках периметра квадрата Р, или, что то же, - в точках краевой линии S, - изменяется по закону (15).

Последовательно находим:

1. Условие (15) позволяет определить функцию f (x,y) во всех точках линии S. Действительно, как бы не изменялась функция F внутри области Р она в точках периметра должна принимать значения, определенные функцией (15), поэтому, подставляя (15) в (14) и проводя дифференцирование, получим вид f (x,y) в точках на линии S:

f(^x, ^y) = e^{-^y} (e^{^x}+1) (16)

2. Выбираем шаг h=1/3. Перенумеровываем узлы области Р так, как показано на рисунке 2, и в точках периметра последовательно вычисляем:

Рис.2	Рис. 3

F _1,1 = 1; F _1,2 = exp (1/3) 1,4; F _1,3 = exp (2/3)  1,95;

F _1,4 = exp (1)  2,7; F _2,1 = exp (-1/3) 0,72; F_3,1 = exp (-2/3) 0,51;

F _4,1 = exp (-1)  0,36; F _4,2 = exp (-2/3) 0,51; F _4,3 = exp (-1/3) 0,72;

F _4,4 = 1; F _3,4 = exp (1/3) 1,4; F _2,4 = exp (2/3)  1,95; `

3. Для значений во внутренних узлах согласно (15) составляем систему из 4-х уравнений с четырьмя неизвестными:

F23 – 2F22 + F21	+	F32 – 2F22 + F12	+ e – 1/3	F23– F21	- F22 = e – 1/3 (e 1/3 + 1)
(1/3)2		(1/3)2		(2/3)
F24 – 2F23 + F23	+	F33 – 2F23 + F13	+ e – 2/3	F24– F22	- F23 = e – 1/3 (e 2/3 + 1)
(1/3)2		(1/3)2		(2/3)
F34 – 2F33 + F32	+	F43 – 2F33 + F23	+ e – 2/3	F34– F32	- F33 = e – 2/3 (e 2/3 + 1)
(1/3)2		(1/3)2		(2/3)
F33 – 2F32 + F31	+	F42 – 2F32 + F22	+ e – 1/3	F33– F31	- F32 = e – 2/3 (e 2/3 + 1)
(1/3)2		(1/3)2		(2/3)

Введем для краткости обозначения: F_2,2  x, F_2,3  y, F_3,2  z, F_3,3  u. Подставляя их в полученную систему и проводя предварительные вычисления, получим систему из 4-х следующих линейных уравнений: –37x + 10,1y +9z = –16,5

+8,2x –37y +9u = –35,45 +9y +8,2z – 37u = –18,7 +9x – 37z+10,1u = –7,42

Программа вычисления корней данной системы приведена ниже:

uses crt; const n=4;

type qw=array[1..n,1..n] of real; Linia=array[1..n]of real;

const MotL:qw= ((-37,10.1,9,0),(8.2,-37,0,9),(0,9,8.2,-37),(9,0,-37,10.1));

BotL:Linia= (-16.5,-35.45,-18.7,-7.42);

var m1,m2:qw; x,b1,b2:Linia; aa,ss,zz:real;i,j,k,q,tt:integer;

Procedure CoeFA(i,j,k:byte);

begin m2[i,j]:=m1[i,j] - m1[i,k]*m1[k,j]/m1[k,k] End;

Procedure FreeB(i,k:byte);begin b2[i]:=b1[i]-b1[k]*m1[i,k]/m1[k,k] End;

BEGIN clrscr; For i:=1 To n Do x[i]:=0; For i:=1 to n Do b1[i]:=BotL[i];

For i:=1 to n Do For j:=1 To n Do m1[i,j]:=MotL[i,j];

For i:=1 to n Do For j:=1 to n Do m2[i,j]:=0;

For j:=1 to n Do m2[1,j]:=m1[1,j]; b2[1]:=b1[1];

For tt:=2 To n Do

Begin For i:=tt to n Do m2[i,tt]:=0;

For i:=tt To n Do For j:=2 to n Do CoeFA(i,j,tt-1);

For i:=tt to n Do FreeB(i,tt-1);

For i:=1 to n Do For j:=1 To n Do m1[i,j]:=m2[i,j];

For i:=1 to n Do b1[i]:=b2[i] End;

zz:=0; X[n]:=b2[n]/m2[n,n];

For i:=n-1 DownTo 1 Do

Begin zz:=b2[i]; q:=n;

For j:=i To n-1 Do Begin zz:=zz-m2[i,q]*x[q]; dec(q); End;

x[i]:=zz/m2[i,i] End;

For i:=1 to n Do WriteLn('X',i,'=',X[i]:6:4); Repeat Until KeyPressed END.

Система имеет следующее решение:

F_2,2  x=1.0119 F_2,3  y=1,4288 F_3,2  z=0,7233 F_3,3  u=1,0133.

2.3. Оценка погрешности дискретной модели

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

Погрешности дискретизации зависят от следующих факторов:

• способа замены дифференциальных уравнений разностными;

• от конфигурации элементов конструкции (формы рассматриваемой области);

• внешних воздействий (граничных условий);

• длительности рассчитываемого процесса.

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

Покажем, как определить порядок аппроксимации на примере замены производных конечными разностями. Допустим, что мы хотим заменить первую производную в точке 0 (рисунок 3) и для этого наметим два узла сетки в точках x=-a и x=h-a. Будем считать функцию F и ее производную в точке 0 известными. Воспользовавшись разложением в ряд Тейлора, находим значения функции на концах отрезка при x=-a и x=h-a:

F(–a) = F –

дF



a

+

д2F



a2

-

д3F



a2

+ …

дX

1!

дX2

2!

дX3

3!

F(h–a) = F +	дF		(h-a)	+	д2F		(h-a)2	+ …
	дX		1!		дX2		2!

Далее определим значение конечной разности:

F(h–a) – F(–a)	=	дF	+	д2F		(h-2a)	+	д3F		h2 – 3ah +3a2	+ …
h		дX		дX2		2!		дX3		3!

Погрешность от замены первой производной конечной разностью будет равна:

F(h–a) – F(–a)	-	дF	=	д2F		(h-2a)	+	д3F		h2 – 3ah +3a2	+ …
h		дX		дX2		2!		дX3		3!

При а=0 разность будет правой, при a=h - левой (9-б и 9-а соответственно). При этом погрешности соответственно составят:

для правой разности:	{	д2F		h	} + {	д3F		h2	}
		дX2		2!		дX3		3!
для левой разности:	{	д2F		h	} + {	д3F		h2	}
		дX2		2!		дX3		3!

В том и в другом случае погрешность пропорциональна шагу сетки, то есть имеет место первый порядок аппроксимации производной конечной разностью. Условно это можно записать в виде:

Fm+1n- Fmn	=	дF	+ O(h) и	Fm,n- Fm-1,n	=	дF	+ O(h)
h		дY		h		дY

Для вторых разностей ошибка замены второй производной может быть определена аналогично. Используя разложение функции F в ряд Тейлора вблизи точки X = mh, можно показать, что здесь имеет место второй порядок аппроксимации.