gdzie:

Jeśli zastosujemy równanie (7.55) do siatki przedstawionej na rysunku 7.3, na każdym kroku czasowym musimy rozwiązać następujący układ szesnastu równań liniowych z szesnastoma niewiadomymi wartościami funkcji w wewnętrznych węzłach siatki

=

, (7.56)

gdzie:

Jest to układ równań z trójdiagonalną macierzą blokową. Macierz ta jest diagonalnie dominująca.

Rys. 7.3

Trudności występujące przy rozwiązywaniu zagadnień początkowo-brzegowych dla wielowymiarowych równań parabolicznych, związane z bardzo rygorystycznymi ograniczeniami dla schematów jawnych i koniecznością rozwiązywania bardzo dużych układów równań w przypadku metod niejawnych, spowodowały poszukiwanie takich bezwarunkowo stabilnych schematów, których koszt realizacji jest proporcjonalny do liczby niewiadomych. Takie schematy nazywają się schematami ekonomicznymi.

Najbardziej popularnymi schematami ekonomicznymi są schematy metody naprzemiennych kierunków (Alternating-Direction Implicit Method). Za pomocą tej metody rozwiązywanie równania (7.46) odbywa się w dwu etapach:

Etap I

(7.57)

Etap II

(7.58)

gdzie i
są aproksymacjami różnicowymi pochodnych oraz
Rezultatem rozszczepienia jest konieczność rozwiązywania tylko układów równań z trójdiagonalnymi macierzami współczynników: najpierw dla każdego wskaźnika j, a następnie dla każdego wskaźnika i. Metoda naprzemiennych kierunków jest metodą drugiego rzędu dokładności z błędem aproksymacji Istnieją również schematy ekonomiczne dla trójwymiarowych równań parabolicznych [40].

*

Zagadnienie, określone równaniem

(7.59)

i warunkami granicznymi:

(7.60)

którego rozwiązaniem dokładnym jest funkcja

(7.61)

jest rozwiązywane w programie 7.1 przy wykorzystaniu czterech rodzajów schematów różnicowych: jawnej metody pierwszego rzędu, niejawnej metody pierwszego rzędu, metody Cranka-Nicolsona i metody Duforta-Frankela. Wymagany dodatkowy warunek początkowy dla schematu Duforta-Frankela jest wyznaczany z rozwiązania dokładnego (7.61).

{Program 7.1}

unit Obliczenia;

interface

uses

Windows, Messages, SysUtils, Classes, Graphics, Controls,

Forms, Dialogs, StdCtrls, Buttons;

type

Tabl = array[0..100] of Real;

. . . . . . . . . . . . . . . .

var

Form3: TForm3;

a,al,al1,al2,b,bl,dt,h,h2,

p,t,tmax,ua,ub,x: Real;

aa,bb,dd,u,un,un1: Tabl;

j,m,m1,nr: Integer;

plik: Text;

implementation

uses Ustawienia, Informacje, Grafika, Podglad;

{$R *.DFM}

function u0(x: Real): Real;

begin

u0:=3*Sin(x);

end;

function ud(x,t: Real): Real;

begin

ud:=3*Exp(-t)*Sin(x);

end;

{procedure Tridiag1(n: Integer; a,b,c,d: Tabl;

var x: Tabl);}

. . . . . . . . . . . . . . . . . . . . . . .

procedure TForm3.BitBtn1Click(Sender: TObject);

begin

Form2.Show;

AssignFile(plik,Edit8.Text);

Rewrite(plik); Writeln(plik,'PROGRAM 7.1.');

Writeln(plik,'Zagadnienie początkowo-brzegowe');

Writeln(plik,'dla równania parabolicznego.');

a:=StrToFloat(Edit1.Text);

b:=StrToFloat(Edit2.Text);

m:=StrToInt(Edit3.Text);

ua:=StrToFloat(Edit4.Text);

ub:=StrToFloat(Edit5.Text);

dt:=StrToFloat(Edit6.Text);

tmax:=StrToFloat(Edit7.Text);

if RadioButton1.Checked then nr:=1;

if RadioButton2.Checked then nr:=2;

if RadioButton3.Checked then nr:=3;

if RadioButton4.Checked then nr:=4;

case nr of

1: Writeln(plik,'Metoda jawna pierwszego rzędu.');

2: Writeln(plik,'Metoda niejawna pierwszego rzędu.');

3: Writeln(plik,'Metoda Cranka-Nicolsona.');

4: Writeln(plik,'Metoda Duforta-Frankela.');

end;

Writeln(plik);

Writeln(plik,'Liczba podprzedziałów: m = ',m:3);

Writeln(plik,'Początek przedziału: a = ',a:13);

Writeln(plik,'Koniec przedziału: b = ',b:13);

Writeln(plik,'Krok czasowy: dt = ',dt:13);

Writeln(plik,'Warunki brzegowe: ua = ',ua:13);

Writeln(plik,' ub = ',ub:13);

Writeln(plik);

h:=(b-a)/m; h2:=h*h;

m1:=m-1; al:=dt/h2;

u[0]:=ua; u[m]:=ub;

un[0]:=ua; un[m]:=ub;

un1[0]:=ua; un1[m]:=ub;

al1:=(1-2*al)/(1+2*al);

al2:=2*al/(1+2*al);

for j:=1 to m1 do begin

x:=a+j*h;

un[j]:=u0(x);

end;

t:=0;

repeat

t:=t+dt;

case nr of

1: for j:=1 to m1 do

u[j]:=un[j]+al*(un[j-1]-2*un[j]+un[j+1]);

2: begin

for j:=1 to m1 do begin

aa[j]:=-al;

bb[j]:=1+2*al;

dd[j]:=un[j];

end;

dd[1]:=dd[1]-aa[1]*ua;

dd[m1]:=dd[m1]-aa[m1]*ub;

Tridiag1(m1,aa,bb,aa,dd,u);

end;

3: begin

for j:=1 to m1 do begin

aa[j]:=-al/2; bb[j]:=1+al;

p:=un[j-1]-2*un[j]+un[j+1];

dd[j]:=un[j]+al*p/2;

end;

dd[1]:=dd[1]-aa[1]*ua;

dd[m1]:=dd[m1]-aa[m1]*ub;

Tridiag1(m1,aa,bb,aa,dd,u);

end;

4: begin

if Abs(t-dt)<1e-10 then

for j:=1 to m1 do

u[j]:=ud(a+j*h,t)

else

for j:=1 to m1 do begin

p:=un[j-1]+un[j+1];

u[j]:=al1*un1[j]+al2*p;

end;

end;

end;

for j:=1 to m1 do begin

un1[j]:=un[j];

un[j]:=u[j];

end;

until (t>=tmax);

Writeln(plik,'Wartości funkcji u(x,t) dla t = ',t:13,':');

Writeln(plik,' j x[j] u[j] błąd');

for j:=0 to m do begin

x:=a+j*h; bl:=ud(x,t)-u[j];

Writeln(plik,j:3,' ',x:13,' ',u[j]:18,' ',bl:13);

end;

Writeln(plik); CloseFile(plik);

Form2.Wyniki.Lines.LoadFromFile(Edit8.Text);

end;

. . . . . . . . . . . . . . . . . . . . . . .

procedure TForm3.BitBtn3Click(Sender: TObject);

begin

Close;

end;

end.

W programie 7.1 występują dwie procedury funkcyjne: u0(x) i ud(x,t) - przeznaczone do obliczania wartości funkcji
i rozwiązania dokładnego
oraz procedura

Tridiag1(n,a,b,c,d,x),

której zadaniem jest rozwiązywanie układu równań (2.76) przy wykorzystaniu algo-rytmu metody faktoryzacji (2.77) ÷ (2.79). Sens wczytywanych danych jest objaśniony odpowiednimi napisami ukazującymi się na formularzu Dane (rys. 7.4).

Rys. 7.4

Wyniki obliczeń numerycznych dla uzyskane za pomocą kolejnych schematów różnicowych przy przyjęciu i kroków czasowych zawierają załączone tabulogramy wyników. Zgodnie z warunkiem (7.39) maksymalny dopuszczalny krok czasowy dla schematu jawnego pierwszego rzędu wynosi

PROGRAM 7.1.

Zagadnienie początkowo-brzegowe

dla równania parabolicznego.

Metoda jawna pierwszego rzędu.

Liczba podprzedziałów: m = 36

Początek przedziału: a = 0.0000E+0000

Koniec przedziału: b = 3.1416E+0000

Krok czasowy: dt = 1.0000E-0003

Warunki brzegowe: ua = 0.0000E+0000

ub = 0.0000E+0000

Wartości funkcji u(x,t) dla t = 1.0000E+0000:

j x[j] u[j] błąd

0 0.0000E+0000 0.000000000E+0000 0.0000E+0000

1 8.7266E-0002 9.620138102E-0002 -1.2963E-0005

2 1.7453E-0001 1.916706114E-0001 -2.5828E-0005

3 2.6180E-0001 2.856811128E-0001 -3.8496E-0005

4 3.4907E-0001 3.775174084E-0001 -5.0871E-0005

5 4.3633E-0001 4.664805685E-0001 -6.2859E-0005

6 5.2360E-0001 5.518935298E-0001 -7.4368E-0005

7 6.1087E-0001 6.331062482E-0001 -8.5312E-0005

8 6.9813E-0001 7.095006457E-0001 -9.5606E-0005

9 7.8540E-0001 7.804953149E-0001 -1.0517E-0004

10 8.7266E-0001 8.455499435E-0001 -1.1393E-0004

11 9.5993E-0001 9.041694264E-0001 -1.2183E-0004

12 1.0472E+0000 9.559076341E-0001 -1.2881E-0004

13 1.1345E+0000 1.000370807E+0000 -1.3480E-0004

14 1.2217E+0000 1.037220555E+0000 -1.3977E-0004

15 1.3090E+0000 1.066176428E+0000 -1.4367E-0004

16 1.3963E+0000 1.087018054E+0000 -1.4648E-0004

17 1.4835E+0000 1.099586817E+0000 -1.4817E-0004

18 1.5708E+0000 1.103787060E+0000 -1.4874E-0004

19 1.6581E+0000 1.099586817E+0000 -1.4817E-0004

20 1.7453E+0000 1.087018054E+0000 -1.4648E-0004

21 1.8326E+0000 1.066176428E+0000 -1.4367E-0004

22 1.9199E+0000 1.037220555E+0000 -1.3977E-0004

23 2.0071E+0000 1.000370807E+0000 -1.3480E-0004

24 2.0944E+0000 9.559076340E-0001 -1.2881E-0004

25 2.1817E+0000 9.041694264E-0001 -1.2183E-0004

26 2.2689E+0000 8.455499435E-0001 -1.1393E-0004

27 2.3562E+0000 7.804953149E-0001 -1.0517E-0004

28 2.4435E+0000 7.095006457E-0001 -9.5606E-0005

29 2.5307E+0000 6.331062482E-0001 -8.5312E-0005

30 2.6180E+0000 5.518935298E-0001 -7.4368E-0005

31 2.7053E+0000 4.664805685E-0001 -6.2859E-0005

32 2.7925E+0000 3.775174084E-0001 -5.0871E-0005

33 2.8798E+0000 2.856811128E-0001 -3.8496E-0005

34 2.9671E+0000 1.916706114E-0001 -2.5828E-0005

35 3.0543E+0000 9.620138102E-0002 -1.2963E-0005

36 3.1416E+0000 0.000000000E+0000 -1.0402E-0011

PROGRAM 7.1.

Zagadnienie początkowo-brzegowe

dla równania parabolicznego.

Metoda niejawna pierwszego rzędu.

Liczba podprzedziałów: m = 36

Początek przedziału: a = 0.0000E+0000

Koniec przedziału: b = 3.1416E+0000

Krok czasowy: dt = 5.0000E-0003

Warunki brzegowe: ua = 0.0000E+0000

ub = 0.0000E+0000

Wartości funkcji u(x,t) dla t = 1.0000E+0000:

j x[j] u[j] błąd

0 0.0000E+0000 0.000000000E+0000 0.0000E+0000

1 8.7266E-0002 9.648928389E-0002 -3.0087E-0004

2 1.7453E-0001 1.922442261E-0001 -5.9944E-0004

3 2.6180E-0001 2.865360736E-0001 -8.9346E-0004

4 3.4907E-0001 3.786472086E-0001 -1.1806E-0003

5 4.3633E-0001 4.678766097E-0001 -1.4589E-0003

6 5.2360E-0001 5.535451873E-0001 -1.7260E-0003

7 6.1087E-0001 6.350009518E-0001 -1.9800E-0003

8 6.9813E-0001 7.116239756E-0001 -2.2189E-0003

9 7.8540E-0001 7.828311113E-0001 -2.4410E-0003

10 8.7266E-0001 8.480804295E-0001 -2.6444E-0003

11 9.5993E-0001 9.068753436E-0001 -2.8278E-0003

12 1.0472E+0000 9.587683887E-0001 -2.9896E-0003

13 1.1345E+0000 1.003364628E+0000 -3.1286E-0003

14 1.2217E+0000 1.040324656E+0000 -3.2439E-0003

15 1.3090E+0000 1.069367185E+0000 -3.3344E-0003

16 1.3963E+0000 1.090271184E+0000 -3.3996E-0003

17 1.4835E+0000 1.102877562E+0000 -3.4389E-0003

18 1.5708E+0000 1.107090375E+0000 -3.4521E-0003

19 1.6581E+0000 1.102877562E+0000 -3.4389E-0003

20 1.7453E+0000 1.090271184E+0000 -3.3996E-0003

21 1.8326E+0000 1.069367185E+0000 -3.3344E-0003

22 1.9199E+0000 1.040324656E+0000 -3.2439E-0003

23 2.0071E+0000 1.003364628E+0000 -3.1286E-0003

24 2.0944E+0000 9.587683888E-0001 -2.9896E-0003

25 2.1817E+0000 9.068753437E-0001 -2.8278E-0003

26 2.2689E+0000 8.480804296E-0001 -2.6444E-0003

27 2.3562E+0000 7.828311114E-0001 -2.4410E-0003

28 2.4435E+0000 7.116239757E-0001 -2.2189E-0003

29 2.5307E+0000 6.350009519E-0001 -1.9800E-0003

30 2.6180E+0000 5.535451874E-0001 -1.7260E-0003

31 2.7053E+0000 4.678766098E-0001 -1.4589E-0003

32 2.7925E+0000 3.786472086E-0001 -1.1806E-0003

33 2.8798E+0000 2.865360736E-0001 -8.9346E-0004

34 2.9671E+0000 1.922442261E-0001 -5.9944E-0004

35 3.0543E+0000 9.648928390E-0002 -3.0087E-0004

36 3.1416E+0000 0.000000000E+0000 -1.0402E-0011

PROGRAM 7.1.

Zagadnienie początkowo-brzegowe

dla równania parabolicznego.

Metoda Cranka-Nicolsona.

Liczba podprzedziałów: m = 36

Początek przedziału: a = 0.0000E+0000

Koniec przedziału: b = 3.1416E+0000

Krok czasowy: dt = 1.0000E-0002

Warunki brzegowe: ua = 0.0000E+0000

ub = 0.0000E+0000

Wartości funkcji u(x,t) dla t = 1.0000E+0000:

j x[j] u[j] błąd

0 0.0000E+0000 0.000000000E+0000 0.0000E+0000

1 8.7266E-0002 9.624866418E-0002 -6.0246E-0005

2 1.7453E-0001 1.917648179E-0001 -1.2003E-0004

3 2.6180E-0001 2.858215256E-0001 -1.7891E-0004

4 3.4907E-0001 3.777029589E-0001 -2.3642E-0004

5 4.3633E-0001 4.667098446E-0001 -2.9213E-0004

6 5.2360E-0001 5.521647866E-0001 -3.4562E-0004

7 6.1087E-0001 6.334174211E-0001 -3.9648E-0004

8 6.9813E-0001 7.098493666E-0001 -4.4433E-0004

9 7.8540E-0001 7.808789298E-0001 -4.8879E-0004

10 8.7266E-0001 8.459655329E-0001 -5.2953E-0004

11 9.5993E-0001 9.046138274E-0001 -5.6624E-0004

12 1.0472E+0000 9.563774645E-0001 -5.9864E-0004

13 1.1345E+0000 1.000862492E+0000 -6.2648E-0004

14 1.2217E+0000 1.037730351E+0000 -6.4956E-0004

15 1.3090E+0000 1.066700455E+0000 -6.6770E-0004

16 1.3963E+0000 1.087552326E+0000 -6.8075E-0004

17 1.4835E+0000 1.100127266E+0000 -6.8862E-0004

18 1.5708E+0000 1.104329573E+0000 -6.9125E-0004

19 1.6581E+0000 1.100127266E+0000 -6.8862E-0004

20 1.7453E+0000 1.087552326E+0000 -6.8075E-0004

21 1.8326E+0000 1.066700455E+0000 -6.6770E-0004

22 1.9199E+0000 1.037730351E+0000 -6.4956E-0004

23 2.0071E+0000 1.000862492E+0000 -6.2648E-0004

24 2.0944E+0000 9.563774645E-0001 -5.9864E-0004

25 2.1817E+0000 9.046138274E-0001 -5.6624E-0004

26 2.2689E+0000 8.459655329E-0001 -5.2953E-0004

27 2.3562E+0000 7.808789298E-0001 -4.8879E-0004

28 2.4435E+0000 7.098493666E-0001 -4.4433E-0004

29 2.5307E+0000 6.334174211E-0001 -3.9648E-0004

30 2.6180E+0000 5.521647866E-0001 -3.4562E-0004

31 2.7053E+0000 4.667098446E-0001 -2.9213E-0004

32 2.7925E+0000 3.777029589E-0001 -2.3642E-0004

33 2.8798E+0000 2.858215256E-0001 -1.7891E-0004

34 2.9671E+0000 1.917648179E-0001 -1.2003E-0004

35 3.0543E+0000 9.624866419E-0002 -6.0246E-0005

36 3.1416E+0000 0.000000000E+0000 -1.0402E-0011

PROGRAM 7.1.

Zagadnienie początkowo-brzegowe

dla równania parabolicznego.

Metoda Duforta-Frankela.

Liczba podprzedziałów: m = 36

Początek przedziału: a = 0.0000E+0000

Koniec przedziału: b = 3.1416E+0000

Krok czasowy: dt = 5.0000E-0003

Warunki brzegowe: ua = 0.0000E+0000

ub = 0.0000E+0000

Wartości funkcji u(x,t) dla t = 1.0000E+0000:

j x[j] u[j] błąd

0 0.0000E+0000 0.000000000E+0000 0.0000E+0000

1 8.7266E-0002 9.593421692E-0002 2.5420E-0004

2 1.7453E-0001 1.911383165E-0001 5.0647E-0004

3 2.6180E-0001 2.848877381E-0001 7.5488E-0004

4 3.4907E-0001 3.764689920E-0001 9.9755E-0004

5 4.3633E-0001 4.651850895E-0001 1.2326E-0003

6 5.2360E-0001 5.503608477E-0001 1.4583E-0003

7 6.1087E-0001 6.313480274E-0001 1.6729E-0003

8 6.9813E-0001 7.075302675E-0001 1.8748E-0003

9 7.8540E-0001 7.783277750E-0001 2.0624E-0003

10 8.7266E-0001 8.432017381E-0001 2.2343E-0003

11 9.5993E-0001 9.016584269E-0001 2.3892E-0003

12 1.0472E+0000 9.532529507E-0001 2.5259E-0003

13 1.1345E+0000 9.975926438E-0001 2.6434E-0003

14 1.2217E+0000 1.034340055E+0000 2.7407E-0003

15 1.3090E+0000 1.063215513E+0000 2.8172E-0003

444 7. Równania różniczkowe cząstkowe

7.3. Równania paraboliczne 445

---