1.2 Реализация граничных условий
1) Граничные условия первого рода.
Для явной разностной схемы значение функции на границах находится непосредственно из граничного условия:
; 
. (1.3)
2) Граничные условия второго рода.
Поскольку порядок аппроксимации разностной схемы (1.2) по пространству равен двум, то и граничные условия должны быть аппроксимированы со вторым порядком точности. Для этого поступают следующим образом.
а)
Заменяют производную в точке
конечно-разностным аналогом со вторым
порядком точности:
или
. (1.4)
Здесь
– значение температуры, записанное в
фиктивном узле
б)
Из предположения, что решение непрерывно
записывают уравнение (1.2) для
:
.
Из
полученного соотношения исключают
с помощью (1.4):
. (1.5)
Соотношение
(1.6) служит для определения значения
температуры на левой границе (в точке
)
в случае граничного условия второго
рода.
Проводя
подобные рассуждения для правого
граничного условия (только в этом случае
уравнение (1.2) записывается для точки
),
получим:
. (1.6)
3) Граничные условия третьего рода.
Для нахождения температуры в крайней левой и крайней правой точках в случае граничного условия третьего рода поступают так же, как и с граничным условием второго рода.
Разностные аналоги граничных условий имеют вид:
; 
(1.7)
Из (1.2), записанного для и , и (1.7) получим, соответственно:
(1.8)
(1.9)
Условие
устойчивости разностной схемы в этом
случае будет иметь вид:
1.3 Особенности построения разностной схемы для уравнения, записанного в цилиндрической системе координат
Рассмотрим параболическое уравнение, записанное в цилиндрической системе координат:
. (1.10)
При построении разностной схемы будем использовать следующий шаблон (рис. 3):
Рис. 3
Производная по времени записывается так же, как и п.1.1: , а производная по пространству:
Здесь
.
Тогда разностная схема для уравнения (1.10) принимает вид:
,
(1.11)
Граничные
условия реализуются так же, как п. 1.2, за
исключением того случая, когда левая
граница есть центр симметрии, т.е.
.
На такой границе должно быть задано
условие симметрии:
.
Для нахождения решения на верхнем
временном слое в узле i
= 0, уравнение (1.11) нельзя применять, как
это делалось в п.1.2 с уравнением (1.2),
поскольку в данном случае существует
особенность при
.
Для
устранения этой особенности перепишем
уравнение (1.10) в виде
.
Устремляя
к нулю и раскрывая по Лопиталю
неопределенность
,
получим следующее уравнение в точке
:
(1.12)
Записывая
разностную схему для (1.12) и для граничного
условия
,
получим следующее уравнение для
нахождения функции на верхнем временном
слое в узле i
= 0:
(1.13)
Условие
устойчивости
.
1.4 Пример решения параболического уравнения методом конечных разностей
Решить краевую задачу:
НУ: 
,
ГУ: 
; 
.
Вывести
график функции
в момент времени
.
Для решения данной задачи используем разностную схему (1.2) и (1.3).
Ниже приведены примеры программ для численного решения поставленной задачи, реализованной для системы Mathcad и написанной на языке Паскаль.
1.4.1 Программа на языке Паскаль
Program Therm;
uses crt,gra1,graph;
Const
{Задание конечного момента времени}
time_end=2.0;
{Задание числа отрезков разбиения области решения}
N=50;
{Расчет шага интегрирования по пространству}
h=pi/2/N;
{Задание числа Куранта}
C=0.5;
{Расчет шага по времени из условия устойчивости}
tau=C*h*h;
Var
i:integer;
Uk1,Uk,X:array [0..N] of real;
time:real;
begin
{Расчет координат узлов по пространству разностной сетки }
for i:=0 to N do X[i]:=h*i;
{Задание начальных уловий}
for i:=1 to N-1 do Uk[i]:=0;
time:=0;
{Построение координатных осей}
CoordOs(0,pi/2,10,0,2.5,10,'Time','Uk');
while time<=time_end do begin
time:=time+tau;
{Вычисление функции в граничных точках на верхнем временном слое}
Uk1[0]:=0; Uk1[N]:=time;
{Вычисление функции во внутренних точках на верхнем временном слое}
for i:=1 to N-1 do
Uk1[i]:=C*Uk[i+1]+(1-2*C)*Uk[i]+C*Uk[i-1]+
(1+time)*sin(i*h)*tau;
{Присвоение переменным Uk, хранящим значения функции
на старом временном слое, значений функции на новом временном слое}
for i:=0 to N do Uk[i]:=Uk1[i];
end;
{Установление графического курсора в начальную точку}
FirstRis; Ris(X[0],Uk[0],13);
{Графический вывод распределения функции вдоль X}
for i:=1 to N do Ris(X[i],Uk[i],13);
repeat
until keypressed;
CloseGraph;
end.
Графический модуль.
Unit Gra1;
INTERFACE
uses Dos,graph;
var gd,gm,i,j,CX,CY,NewX,NewY :integer;
X_0,Y_0,Dx,Dy :real;
Otr,MasX,MasY :real;
St :string[7];
First :boolean;
procedure FirstRis;
procedure Coordos(X0,Xk:real;Nx:integer;Y0,Yk:real;Ny:integer;StrX,StrY:string);
procedure ris(x,y:real;p1:integer);
IMPLEMENTATION
procedure FirstRis;
begin
First:=true;
end;
function CoorX(Xtek:real):integer;
begin
CoorX:=70+round(MasX*(Xtek-X_0));
end;
function CoorY(Ytek:real):integer;
begin
CoorY:=GetMaxY-100-round(MasY*(Ytek-Y_0));
end;
procedure CoordOs(X0,Xk:real;Nx:integer;Y0,Yk:real;Ny:integer;StrX,StrY:string);
begin
gd:=detect;
InitGraph(gd,gm,'C:\BP\bgi');
{Определение графического окна}
SetViewPort(0,0,GetMaxX,GetMaxY,ClipOn);
ClearViewPort;
SetBkColor(0);
SetColor(11);
X_0:=X0;Y_0:=Y0;
Dx:=(Xk-X0)/Nx;
MasX:=(GetMaxX-120)/(Xk-X0);
Dy:=(Yk-Y0)/Ny;
MasY:=(GetMaxY-150)/(Yk-Y0);
{Вывод сетки и надписей по осям}
SetColor(11);
for i:=0 to Nx do
begin
Otr:=X0+i*Dx;
CX:=CoorX(Otr);
Line(CX,50,CX,GetMaxY-100);
Str(Otr:6:2,St);
OutTextXY(CX-40,GetMaxY-90,St);
end;
for j:=0 to Ny do
begin
Otr:=Y0+j*Dy;
CY:=CoorY(Otr);
Line(70,CY,GetMaxX-50,CY);
Str(Otr:6:0,St);
OutTextXY(15,CY-5,St);
end;
OutTextXY(GetMaxX-40,GetMaxY-105,Strx);
OutTextXY(60,30,Stry);
end{CoordOs};
procedure Ris(x,y:real;p1:integer);
begin
SetColor(p1);
SetLineStyle(0,1,1);
Newx:=CoorX(x);
Newy:=CoorY(y);
if First then
begin
MoveTo(NewX,NewY);
First:=false;
end
else
LineTo(NewX,NewY);
end{Ris};
end.
Результат расчета по приведенной программе приведен на рис. 4.
Рис.
4