Решить краевую задачу

. (3.13)

Граничные условия:

Решение уравнение (3.1) будем находить с помощью разностной схемы (3.3) - (3.5). Область решения покрывается разностной сеткой из N ´ M =400 точек, координаты которых будут вычисляться по формулам: , , , , где , .

В соответствие со схемой (3.3-3.4) функция во внутренних точках области на k +1 итерации будет определяться по формуле:

,

а в граничных точках по формулам:

.

Ниже приведены примеры программ для численного решения поставленной задачи, реализованной для системы Mathcad, и написанной на языке Паскаль.

3.5.1. Программа на языке Паскаль.

Program Two_Therm_Dirichle;

Uses Crt;

Const

N=20; M=20;

u01=100; u02=400;

u03=100; u04=400;

u00=100;

Omega=1.75;

{Вычисление шага разностной сетки}

Hx=1.0/N; Hy=1.0/M;

Type

Arr2=Array [0..N,0..M] of real;

ArrX=Array [0..N] of real;

Arry=Array [0..M] of real;

Var u0,u1,u:Arr2;

X:ArrX;

Y:ArrY;

i,j,N_iter:integer;

Alfa,Betta,Gamma,Delta:real;

f:text;

Begin

{ClrScr; }

{Задание краевых условий}

for i:=1 to N-1 do begin u[i,0]:=u01;u[i,M]:=u03 end;

for j:=1 to M-1 do begin u[0,j]:=u04;u[N,j]:=u02 end;

{Задание начального приближения}

for i:=1 to N-1 do

for j:=1 to M-1 do u[i,j]:=u00;

{Вычисление координат сетки}

for i:=0 to N do X[i]:=i*Hx;

for j:=0 to M do Y[j]:=j*Hy;

Alfa:=0.5/(1.0+Sqr(Hx/Hy));

Betta:=0.5/(1+Sqr(Hy/Hx));

Gamma:=0.5*Sqr(Hx)/(1.0+Sqr(Hx/Hy));

N_iter:=0;

repeat

{Присваивание значения температуры на k-ой итерации}

for j:=0 to M do

for i:=0 to N do u0[i,j]:=u[i,j];

{Вычисление значения температуры на k-ой итерации}

for j:=1 to M-1 do

for i:=1 to N-1 do

u[i,j]:=((u[i+1,j]+u[i-1,j])*Alfa+

(u[i,j+1]+u[i,j-1])*Betta)*Omega+u0[i,j]*(1-Omega);

{Вычисление критерия останова итерационного процесса}

Delta:=0;

for j:=1 to M-1 do

for i:=1 to N-1 do Delta:=Delta+abs(u[i,j]-u0[i,j]);

Delta:=Delta/(M-1)/(N-1);

N_iter:=N_iter+1;

until delta<1e-5;

{Запись результатов на диск}

writeln(Omega:6:2,' ',N_iter);

Assign(f,'result.dat');

Rewrite(f);

for j:=0 to M do

for i:=0 to N do writeln(f,X[i],y[j],u[i,j]);

close(f);

End.

Результаты расчетов по приведенной программе – изолинии температуры внутри пластины приведены на рис.10.

Рис.10. Изолинии функции

Поиск по сайту: