86104 (Чисельні методи розв’язування крайових задач для звичайних диференціальних рівнянь), страница 4

Таким чином: маємо

тобто і при

Зауваження:

Якщо яка-небудь дана нам умова апроксимується точно, то стійкість по ній можна не вимагати, тому що вона не вносить похибки у розв’язок (окрім помилок округлення, тоді стійкість по цим даним потрібна).

Для умовної апроксимації (чи стійкості) збіжність теж носить умовний характер.

Програмна реалізація(представлена на мові Delphi)

Розв’язати диференційне рівняння:

З крайовими умовами:

Розв’язання з використанням методу Гауса:

unit Unit1;

interface

uses

Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,

Dialogs, ExtCtrls, StdCtrls, Buttons;

type

TForm1 = class(TForm)

Panel1: TPanel;

Label1: TLabel;

Image1: TImage;

Image2: TImage;

Label2: TLabel;

LabeledEdit1: TLabeledEdit;

LabeledEdit2: TLabeledEdit;

LabeledEdit3: TLabeledEdit;

LabeledEdit4: TLabeledEdit;

LabeledEdit5: TLabeledEdit;

LabeledEdit6: TLabeledEdit;

LabeledEdit7: TLabeledEdit;

LabeledEdit8: TLabeledEdit;

LabeledEdit9: TLabeledEdit;

LabeledEdit10: TLabeledEdit;

LabeledEdit11: TLabeledEdit;

Label3: TLabel;

Label4: TLabel;

SpeedButton1: TSpeedButton;

LabeledEdit12: TLabeledEdit;

Label5: TLabel;

Image3: TImage;

procedure FormCreate(Sender: TObject);

procedure SpeedButton1Click(Sender: TObject);

private

{ Private declarations }

public

{ Public declarations }

end;

var

Form1: TForm1;

type Dynmas=array of array of real;

dynvec=array of real;

var a,b,pi,qi,fi,a1,a2,b1,b2,AA,BB:real;

eps,h:real;

c:dynmas;

st,m,i:integer;

x,d,y,memory:dynvec;

t_all,tx,ty,k_i:textfile;

g:boolean;

str:string;

implementation

uses Unit2;

{$R *.dfm}

function Gauss(n:Integer; a:dynmas; b:dynVec; var x:dynVec):Boolean;

Var i,j,k,l:Integer;

q,m,t:real;

Begin

for k:=0 to n-2 do

begin

l:=-1;

m:=0;

for i:=k to n-1 do

if Abs(a[i, k])>m then

begin

m:=Abs(a[i, k]);

l:=i;

end;

if l=-1 then

begin

Gauss:=false;

Exit;

end;

if l<>k then

begin

For j:=0 to n-1 do

begin

t:=a[k,j];

a[k,j]:=a[l,j];

a[l,j]:=t;

end;

t:=b[k];

b[k]:=b[l];

b[l]:=t;

end;

for i:=k+1 to n-1 do

begin

q:=a[i,k]/a[k,k];

for j:=0 to n-1 do

If j=k then

a[i,j]:= 0

else

a[i,j]:= a[i,j]-q*a[k,j];

b[i]:=b[i]-q*b[k];

end;

end;

if a[n-1,n-1]<>0 then

x[n-1]:=b[n-1]/a[n-1,n-1]

else

begin

Gauss:=false;

Exit;

end;

for i:=n-2 downto 0 do

begin

t:=0;

for j:=1 to n-i do

t:=t+a[i,i+j]*x[i+j];

x[i]:=(1/a[i,i])*(b[i]-t);

end;

Gauss := true;

end;

procedure Koef(var s:dynmas; k:integer; h:real; v:dynvec; var z:dynvec);

var i:integer;

begin

s[0,0]:=h*a1-a2; s[0,1]:=a2;

z[0]:=h*AA;

for i:=0 to 2*(k-1) do

begin

s[i+1,i]:=1-(h*pi*ln(v[i]))/2;

s[i+1,i+1]:=h*h*qi-2;

s[i+1,i+2]:=1+(h*pi*ln(v[i]))/2;

z[i+1]:=h*h*fi;

end;

s[2*k,2*k-1]:=-b2; s[2*k,2*k]:=h*b1+b2;

z[2*k]:=h*BB;

end;

procedure TForm1.FormCreate(Sender: TObject);

begin

getdir(0,str);

str:=str+'\otv\';

end;

procedure TForm1.SpeedButton1Click(Sender: TObject);

begin

if (form1.LabeledEdit1.Text='') and

(form1.LabeledEdit9.Text='') and

(form1.LabeledEdit12.Text='') then

begin

showmessage('так як ви не ввели коефіцієнти, то программа буде задіяна зі стандартним набором данних');

pi:=-1;

qi:=-2;

fi:=1;

a1:=1;

a2:=-1;

a:=0.5;

AA:=1;

b1:=1;

b2:=1;

b:=1.5;

BB:=0;

eps:=0.0001;

end

else

begin

pi:=strtofloat(form1.LabeledEdit1.Text);

qi:=strtofloat(form1.LabeledEdit2.Text);

fi:=strtofloat(form1.LabeledEdit3.Text);

a1:=strtofloat(form1.LabeledEdit4.Text);

a2:=strtofloat(form1.LabeledEdit5.Text);

a:=strtofloat(form1.LabeledEdit6.Text);

AA:=strtofloat(form1.LabeledEdit7.Text);

b1:=strtofloat(form1.LabeledEdit8.Text);

b2:=strtofloat(form1.LabeledEdit9.Text);

b:=strtofloat(form1.LabeledEdit10.Text);

BB:=strtofloat(form1.LabeledEdit11.Text);

eps:=strtofloat(form1.LabeledEdit12.Text);

end;

form2.Series1.Clear;

AssignFile(t_all,str+'otv.txt');

AssignFile(tx,str+'otv_x.txt');

AssignFile(ty,str+'otv_y.txt');

AssignFile(k_i,str+'otv_krok_vuzl.txt');

Rewrite(t_all);

m:=1;

g:=false;

While not g do

begin

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

SetLength(y,2*m+1);

SetLength(x,2*m+1);

SetLength(d,2*m+1);

for i:=0 to 2*m do

x[i]:=a+i*h;

Setlength(c,2*m+1);

for i:=0 to 2*m do

Setlength(c[i],2*m+1);

Koef(c,m,h,x,d);

if gauss(2*m+1,c,d,y)<>true then

break;

if m<>1 then

for i:=0 to m do

if abs(memory[i]-y[2*i])/15>eps then

begin

g:=false;

break;

end

else

g:=true;

SetLength(memory,2*m+1);

memory:=Copy(y);

if g then

writeln(t_all,'Крайова задача розвязана з точністю eps =',eps:0:4);

for i:=0 to 2*m do

begin

write(t_all,y[i]:0:10);

write(t_all,' ');

writeln(t_all,x[i]:0:10);

end;

Writeln(t_all,'Кількість вузлів - ',2*m+1);

Writeln(t_all,'Крок сітки - ',h:0:10);

Writeln(t_all);

st:=m;

m:=m*2;

end;

rewrite(ty);

rewrite(tx);

rewrite(k_i);

writeln(k_i,h:0:10);

writeln(k_i,2*m+1);

form2.StringGrid1.ColCount:=2*st+2;

for i:=0 to (2*st+1) do

begin

form2.StringGrid1.Cells[i+1,0]:=inttostr(i+1);

form2.StringGrid1.Cells[i+1,1]:=floattostr(x[i]);

form2.StringGrid1.Cells[i+1,2]:=floattostr(y[i]);

writeln(ty,y[i]:0:10);

writeln(tx,x[i]:0:10);

end;

for i:=0 to (2*st) do

form2.Series1.AddXY(x[i],y[i]);

form2.Label1.Caption:='Крок сітки - '+floattostr(h);

form2.Label2.Caption:='Кількість вузлів - '+floattostr(2*st+1);

CloseFile(t_all);

CloseFile(tx);

CloseFile(ty);

CloseFile(k_i);

form2.Show;

end;

end.

Результати записуємо у файл.

Графік отриманий програмою:

Розв’язання з використанням методу прогонки:

unit Unit1;

interface

uses

Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,

Dialogs, ExtCtrls, StdCtrls, Buttons;

type

TForm1 = class(TForm)

Panel1: TPanel;

Label1: TLabel;

Image1: TImage;

Image2: TImage;

Label2: TLabel;

LabeledEdit1: TLabeledEdit;

LabeledEdit2: TLabeledEdit;

LabeledEdit3: TLabeledEdit;

LabeledEdit4: TLabeledEdit;

LabeledEdit5: TLabeledEdit;

LabeledEdit6: TLabeledEdit;

LabeledEdit7: TLabeledEdit;

LabeledEdit8: TLabeledEdit;

LabeledEdit9: TLabeledEdit;

LabeledEdit10: TLabeledEdit;

LabeledEdit11: TLabeledEdit;

Label3: TLabel;

Label4: TLabel;

SpeedButton1: TSpeedButton;

LabeledEdit12: TLabeledEdit;

Label5: TLabel;

Image3: TImage;

procedure FormCreate(Sender: TObject);

procedure SpeedButton1Click(Sender: TObject);

private

{ Private declarations }

public

{ Public declarations }

end;

var

Form1: TForm1;

type Dynmas=array of array of real;

dynvec=array of real;

var a,b,pi,qi,fi,a1,a2,b1,b2,AA,BB:real;

eps,h:real;

c:dynmas;

st,m,i:integer;

w_,v_,x,d,y,memory:dynvec;

t_all,tx,ty,k_i:textfile;

g:boolean;

time1,time2,vremja:longint;

str:string;

implementation

uses Unit2;

{$R *.dfm}

Function Timer:longint;

const c60:longint=60;

var h,m,s,s100:word;

begin

decodetime(now,h,m,s,s100);

timer:=((h*c60+m)*c60+s)*100+s100;

end;

function progonka(n:Integer; a:dynmas; b:dynVec; var x:dynVec):boolean;

Var i,j,k,l:Integer;

q,m,t:real;

ls:integer;

Begin

{прямой ход}

w_[0]:=(-a[0,1]/a[0,0]);

v_[0]:=(d[0]/a[0,0]);

for i:=1 to n-1 do

begin

w_[i]:=-(a[i,i+1]/(a[i,i-1]*w_[i-1]+a[i,i]));

v_[i]:=(d[i]-a[i,i-1]*v_[i-1])/(a[i,i-1]*w_[i-1]+a[i,i]);

end;

{w_[n]:= ;

v_[n]:= ;}

for i:=0 to n-1 do

begin

x[i]:=v_[i]+w_[i]*x[i+1];

end;

x[n-1]:=v_[n-1];

{обратный ход}

x[n-1]:=v_[n-1];

for i:=n-1 downto 0 do

begin

x[i]:=w_[i]*x[i+1]+v_[i];

end;

{for k:=0 to n-2 do

begin

l:=-1;

m:=0;

for i:=k to n-1 do

if Abs(a[i, k])>m then

begin

m:=Abs(a[i, k]);

l:=i;

end;

if l=-1 then

begin

progonka:=false;

Exit;

end;

if l<>k then

begin

For j:=0 to n-1 do

begin

t:=a[k,j];

a[k,j]:=a[l,j];

a[l,j]:=t;

end;

t:=b[k];

b[k]:=b[l];

b[l]:=t;

end;

for i:=k+1 to n-1 do

begin

q:=a[i,k]/a[k,k];

for j:=0 to n-1 do

If j=k then

a[i,j]:= 0

else

a[i,j]:= a[i,j]-q*a[k,j];

b[i]:=b[i]-q*b[k];

end;

end;

if a[n-1,n-1]<>0 then

x[n-1]:=b[n-1]/a[n-1,n-1]

else

begin

progonka:=false;

Exit;

end;

for i:=n-2 downto 0 do

begin

t:=0;

for j:=1 to n-i do

t:=t+a[i,i+j]*x[i+j];

x[i]:=(1/a[i,i])*(b[i]-t);

end;}

progonka := true;

end;

procedure Koef(var s:dynmas; k:integer; h:real; v:dynvec; var z:dynvec);

var i:integer;

begin

s[0,0]:=h*a1-a2; s[0,1]:=a2;

z[0]:=h*AA;

for i:=0 to 2*(k-1) do

begin

s[i+1,i]:=1-(h*pi*ln(v[i]))/2;

s[i+1,i+1]:=h*h*qi-2;

s[i+1,i+2]:=1+(h*pi*ln(v[i]))/2;

z[i+1]:=h*h*fi;

end;

s[2*k,2*k-1]:=-b2; s[2*k,2*k]:=h*b1+b2;

z[2*k]:=h*BB;

end;

procedure TForm1.FormCreate(Sender: TObject);

begin

getdir(0,str);

str:=str+'\otv\';

vremja:=0;

end;

procedure TForm1.SpeedButton1Click(Sender: TObject);

begin

if (form1.LabeledEdit1.Text='') and

(form1.LabeledEdit9.Text='') and

(form1.LabeledEdit12.Text='') then

begin

showmessage('так як ви не ввели коефіцієнти, то программа буде задіяна зі стандартним набором данних');

pi:=-1;

qi:=-2;

fi:=1;

a1:=1;

a2:=-1;

a:=0.5;

AA:=1;

b1:=1;

b2:=1;

b:=1.5;

BB:=0;

eps:=0.0001;

end

else

begin

pi:=strtofloat(form1.LabeledEdit1.Text);

qi:=strtofloat(form1.LabeledEdit2.Text);

fi:=strtofloat(form1.LabeledEdit3.Text);

a1:=strtofloat(form1.LabeledEdit4.Text);

a2:=strtofloat(form1.LabeledEdit5.Text);

a:=strtofloat(form1.LabeledEdit6.Text);

AA:=strtofloat(form1.LabeledEdit7.Text);

b1:=strtofloat(form1.LabeledEdit8.Text);

b2:=strtofloat(form1.LabeledEdit9.Text);

b:=strtofloat(form1.LabeledEdit10.Text);

BB:=strtofloat(form1.LabeledEdit11.Text);

eps:=strtofloat(form1.LabeledEdit12.Text);

end;

time2:=timer;

form2.Series1.Clear;

AssignFile(t_all,str+'otv.txt');

AssignFile(tx,str+'otv_x.txt');

AssignFile(ty,str+'otv_y.txt');

AssignFile(k_i,str+'otv_krok_vuzl.txt');

Rewrite(t_all);

m:=1;

g:=false;

While not g do

begin

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

SetLength(y,2*m+1);

SetLength(x,2*m+1);

SetLength(d,2*m+1);

SetLength(w_,2*m+1);

SetLength(v_,2*m+1);

for i:=0 to 2*m do

x[i]:=a+i*h;

Setlength(c,2*m+1);

for i:=0 to 2*m do

Setlength(c[i],2*m+1);

Koef(c,m,h,x,d);

if progonka(2*m+1,c,d,y)<>true then