Методы решения дифференциальных уравнений с частными производными

Для решения уравнений с частными производными самым распространенным являются метод конечных разностей и метод конечных элементов.

Метод конечных разностей. В методе конечных разностей производные первого и второго порядка заменяются на конечные разности функцией:

; (4.6.1)

;

где h-шаг изменения по величине x.

Применительно к уравнению теплопроводности

(4.6.2)

где T=f(x, y, t)

Формулы конечных разностей для явной схемы будут иметь вид:

;

;

; (4.6.3)

;

;

где -температура в точке с координатами i, j по осям x, y соответственно, n-координата по времени.

Уравнение (4.6.2) после подстановки формул (4.6.3) будет иметь вид:

(4.6.4)

где Dx и Dy шаги по величинам x и y. Для реализации схемы (4.6.4) с равномерным шагом плоскость по которой ищут решение T(x, y, t) разбивается квадратной сеткой с шагом Dx=Dy.

Для решения задачи (4.6.2) для прямоугольной области задают начальные и граничные условия:

;

;

где и , а

È и È ;

область 1- незаштрехованная часть плоскости, область 2- заштрехованная часть плоскости.

Пример 1. Раписать схему решения уравнения

для плоскости и начальных данных для области 1- =500°С и для области 2- =100°С.

Исходя из (4.6.4) для каждой узла итерационные уравнения будут иметь вид:

; ;

; ;

; ;

; ;

; ;

;

; ;

; ;

; ;

; ;

; ;

; ;

; ;

; ;

; ;

; ;

; ;

; ;

; ;

; ;

; ;

; ;

; ;

; ;

; ;

Пример 2. Написать алгоритм и программу решения двухмерного уравнения теплопроводности конечно-разностным методом для пространства:

1 1 1 1 1 1 1 1 1 1 1

1 1 1 1 1 1 1 1 1 1 1

1 1 1 1 1 1 1 1 1 1 1

1 1 1 1 2 2 2 1 1 1 1

1 1 1 1 2 2 2 1 1 1 1

1 1 1 1 2 2 2 1 1 1 1

1 1 1 2 2 2 2 2 1 1 1

1 1 1 2 2 2 2 2 1 1 1

1 1 1 2 2 2 2 2 1 1 1

1 1 1 1 2 2 2 1 1 1 1

1 1 1 1 2 2 2 1 1 1 1

1 1 1 1 2 2 2 1 1 1 1

1 1 1 1 1 1 1 1 1 1 1

1 1 1 1 1 1 1 1 1 1 1

1 1 1 1 1 1 1 1 1 1 1

используя следующие параметры для типов клеток:

для типа 1: теплоемкость c=519, теплопроводность l=44.4, плотность r=1600, начальная температура t=100;

для типа 2: теплоемкость c=343, теплопроводность l=110, плотность r=8900, начальная температура t=1100;

с шагом по пространству 0.002 и по времени 0.00005; рассчитать температурное поле до времени 0.1.

uses crt;

var

f: text;

pole, tip, temp: array[1..15, 1..11] of real;

h, t, delta, A, B, C, D, E: real;

l1, l2, c1, c2, r1, r2, t1, t2, point: real;

i, j: integer;

Begin

clrscr;

assign(f, 'tip.txt');

reset(f);

for i: =1 to 15 do begin

for j: =1 to 11 do

read(f, tip[i, j]);

readln(f);

end;

close(f);

assign(f, 'dat.txt');

reset(f);

readln(f, c1);

readln(f, c2);

readln(f, r1);

readln(f, r2);

readln(f, l1);

readln(f, l2);

readln(f, h);

readln(f, t);

readln(f, delta);

readln(f, t1);

readln(f, t2);

close(f);

for i: =1 to 15 do

for j: =1 to 11 do

if tip[i, j]=1 then pole[i, j]: =t1 else pole[i, j]: =t2;

point: =t;

while point< =delta do begin

clrscr;

for i: =2 to 14 do

for j: =2 to 10 do begin

if tip[i, j]=1 then begin

A: =(l1*t)/(c1*r1*h*h);

B: =A;

C: =-4*A+1;

D: =A;

E: =A;

end else begin

A: =(l2*t)/(c2*r2*h*h);

B: =A;

C: =-4*A+1;

D: =A;

E: =A;

end;

temp[i, j]: =A*pole[i+1, j]+B*pole[i, j+1] +C*pole[i, j]+D*pole[i-1, j]+E*pole[i, j-1];

end;

for i: =2 to 14 do

for j: =2 to 10 do

pole[i, j]: =temp[i, j];

point: =point+t;

for i: =1 to 15 do begin

for j: =1 to 11 do

write(' ', pole[i, j]: 5: 1);

writeln(' ');

end;

readln;

end;

End.

Результаты задачи 6:

Файл tip.txt:

1 1 1 1 1 1 1 1 1 1 1

1 1 1 1 1 1 1 1 1 1 1

1 1 1 1 1 1 1 1 1 1 1

1 1 1 1 2 2 2 1 1 1 1

1 1 1 1 2 2 2 1 1 1 1

1 1 1 1 2 2 2 1 1 1 1

1 1 1 2 2 2 2 2 1 1 1

1 1 1 2 2 2 2 2 1 1 1

1 1 1 2 2 2 2 2 1 1 1

1 1 1 1 2 2 2 1 1 1 1

1 1 1 1 2 2 2 1 1 1 1

1 1 1 1 2 2 2 1 1 1 1

1 1 1 1 1 1 1 1 1 1 1

1 1 1 1 1 1 1 1 1 1 1

1 1 1 1 1 1 1 1 1 1 1

Файл dat.txt:

44.4

0.002

0.00005

0.1

Поле температур после окончания расчета:

100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0

100.0 109.4 128.2 161.6 200.4 217.1 200.4 161.6 128.2 109.4 100.0

100.0 123.9 171.0 255.2 354.0 396.6 354.0 255.2 171.0 123.9 100.0

100.0 144.2 229.1 379.0 563.0 641.8 563.0 379.0 229.1 144.2 100.0

100.0 167.0 288.0 488.3 719.0 815.5 719.0 488.3 288.0 167.0 100.0

100.0 193.4 349.6 584.2 821.1 914.2 821.1 584.2 349.6 193.4 100.0

100.0 219.8 409.8 676.7 900.1 979.2 900.1 676.7 409.8 219.8 100.0

100.0 231.0 435.1 715.3 930.5 1002.7 930.5 715.3 435.1 231.0 100.0

100.0 219.8 409.8 676.7 900.1 979.2 900.1 676.7 409.8 219.8 100.0

100.0 193.4 349.6 584.2 821.1 914.2 821.1 584.2 349.6 193.4 100.0

100.0 167.0 288.0 488.3 719.0 815.5 719.0 488.3 288.0 167.0 100.0

100.0 144.2 229.1 379.0 563.0 641.8 563.0 379.0 229.1 144.2 100.0

100.0 123.9 171.0 255.2 354.0 396.6 354.0 255.2 171.0 123.9 100.0

100.0 109.4 128.2 161.6 200.4 217.1 200.4 161.6 128.2 109.4 100.0

100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0 100.0