Автоматихация производства ЭВА (1075779), страница 3
Текст из файла (страница 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) для всех точек области Р, получим систему алгебраических уравнений с числом неизвестных равным числу дискрет области:
f21 + f12 – 4f11 + f01 + f10 = f(1,1)
f31 + f22 – 4f21 + f21 + f21 = f(2,1) (13)
fm+1,n + fm,n+1 – 4fmn + fm-1,n + fm,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 y1} с краевым условием:
(^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; F3,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) |
Введем для краткости обозначения: F2,2 x, F2,3 y, F3,2 z, F3,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.
Система имеет следующее решение:
F2,2 x=1.0119 F2,3 y=1,4288 F3,2 z=0,7233 F3,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, можно показать, что здесь имеет место второй порядок аппроксимации.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
