Чисельні методи розвязування крайових задач для звичайних диференціальних рівнянь (стр. 7 из 8)
("Апроксимація + Стійкість =>Збіжність").
Доведення: Запишемо нев’язку різницевої схеми (3)-(4).
(*)Функція u(x) задовольняє задачі (*) — збуреній задачі (3)-(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
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);