10. Григоренко Я. М., Панкратова Н. Д. Обчислювальні методи в задачах прикладної математики. – К.: Либідь, 1995. – 280 с.
Додаток.
Інструкція користувача та
тексти програм.
Для проведення обчислювальних експериментів по інтерполюванню функцій двох змінних було складено програму, яка будує двовимірний інтерполяційний многочлен (у формі Лагранжа) і двовимірний інтерполяційний ланцюговий дріб з подальшою перевіркою на точність наближення. Проміжки інтерполювання і кількість точок розбиття проміжку по х і по у а також кількість контрольних точок розбиття по кожній змінній (для оцінки похибки) задаються в програмі. На виході програма генерує текстовий файл з максимальними абсолютними і відносними похибками наближення. Функція двох дійсних змінних, яку потрібно інтерполювати, задається безпосередньо в текстах програм в функції Func(). Константи MaxXiMaxY визначають максимальну кількість точок розбиття по відповідних змінних.
Текст програми :
{$A+,B+,D+,E+,F-,G-,I+,L+,N-,O-,P-,Q+,R+,S+,T-,V+,X+}
{$M 65520,0,655360}
Uses Crt;
Const MaxX=50;
MaxY=50;
Type MyArr=Array[0..MaxX,0..MaxY] Of Real;
Var Nx,Ny,Cx,Cy:Integer;
X:Array[0..MaxX] Of Real;
Y:Array[0..MaxY] Of Real;
B:MyArr;
Xa,Xb,Ya,Yb:Real;
D1,D2:^MyArr;
cc,cc1:Integer;
Function Func(x,y:Real):Real;
Begin
Func:=1/(x*x+y*y+x*y);
End;
Procedure DataInput;
Var i,j:Integer;
Begin
{ Write('Input Xa : '); ReadLn(Xa);
Write('Input Xb : '); ReadLn(Xb);
Write('Input Ya : '); ReadLn(Ya);
Write('Input Yb : '); ReadLn(Yb);}
Xa:=1; Xb:=2; Ya:=1; Yb:=2;
{ Write('Input Nx : '); ReadLn(Nx);
Write('Input Ny : '); ReadLn(Ny);}
nx:=cc; ny:=cc1*2-1;
{ For i:=0 To Nx Do X[i]:=(Xa+Xb)/2+(Xb-Xa)*Cos(Pi*i/Nx)/2;
For i:=0 To Ny Do Y[i]:=(Ya+Yb)/2+(Yb-Ya)*Cos(Pi*i/Ny)/2;}
For i:=0 To Nx Do X[i]:=Xa+(Xb-Xa)*i/Nx;
For i:=0 To Ny Do Y[i]:=Ya+(Yb-Ya)*i/Ny;
End;
Procedure BuildCoefTable;
Function Xij(i,j:Integer):Real;
Begin
If i>j Then Xij:=X[i]-X[j] Else Xij:=1;
End;
Function Yij(i,j:Integer):Real;
Begin
If i>j Then Yij:=Y[i]-Y[j] Else Yij:=1;
End;
Function Teta(t,s:Integer):Integer;
Begin
If s>t Then Teta:=-1 Else Teta:=0;
End;
Function Delta(k,i,j:Integer):Real;
Begin
Delta:=Xij(i,k)*Yij(j,k)/
( D1^[i,j]+
Teta(k,j)*D1^[i,k]+
Teta(k,i)*D1^[k,j]+
Teta(k,i)*Teta(k,j)*D1^[k,k]
);
End;
Var i,j,s,k,Mx:Integer;
Begin
For i:=0 To Nx Do
For j:=0 To Ny Do
Begin
D1^[i,j]:=Func(X[i],Y[j]);
End;
k:=0;
D2^:=D1^;
If Nx>Ny Then Mx:=Nx Else Mx:=Ny;
While k<Mx+1 Do
Begin
For i:=0 To Nx Do
For j:=0 To Ny Do
Begin
If i>j Then s:=i Else s:=j;
If s=k Then B[i,j]:=D2^[i,j];
End;
For i:=0 To Nx Do
For j:=0 To Ny Do
Begin
D2^[i,j]:=Delta(k,i,j);
End;
D1^:=D2^;
k:=k+1;
End;
End;
Function Drib(xx,yy:Real):Real;
Var n:Integer;
Function GetH(m,k:Integer):Real;
Begin
If m=n+1 Then GetH:=0
Else
Begin
GetH:=(xx-X[m-1])/(B[m,k]+GetH(m+1,k));
End;
End;
Function GetL(m,k:Integer):Real;
Begin
If m=n+1 Then GetL:=0
Else
GetL:=(yy-Y[m-1])/(B[k,m]+GetL(m+1,k));
End;
Function GetG(k:Integer):Real;
Begin
If k=n+1 Then GetG:=0
Else
GetG:=(xx-X[k-1])*(yy-Y[k-1])/
(B[k,k]+GetH(k+1,k)+GetL(k+1,k)+GetG(k+1));
End;
Begin
If Nx<Ny Then n:=Nx Else n:=Ny;
Drib:=B[0,0]+GetH(1,0)+GetL(1,0)+GetG(1);
End;
Function Polinom(xx,yy:Real):Real;
Var p,q,s,s1,p1,q1:Real; i,j,k:Integer;
Begin
s:=0;
For i:=0 To Nx Do
For j:=0 To Ny Do
Begin
p:=1; q:=1;
For k:=0 To Nx Do If k<>i Then p:=p*(xx-X[k])/(X[i]-X[k]);
For k:=0 To Ny Do If k<>j Then q:=q*(yy-Y[k])/(Y[j]-Y[k]);
s1:=p*q*Func(X[i],Y[j]);
s:=s+s1;
End;
Polinom:=s;
End;
Procedure GetMaxError;
Var i,j:Integer; dx,dy,MaxErr1,p1,p2,p3,VidnErr1,MaxErr2,VidnErr2:Real; F:Text;
Begin
MaxErr1:=0; VidnErr1:=0; MaxErr2:=0; VidnErr2:=0;
dx:=(Xb-Xa)/Cx; dy:=(Yb-Ya)/Cy;
For i:=0 To Cx Do
For j:=0 To Cy Do
Begin
p1:=Func(Xa+i*dx,Ya+j*dy);
p2:=Drib(Xa+i*dx,Ya+j*dy);
p3:=Polinom(Xa+i*dx,Ya+j*dy);
If Abs(p1-p3)>MaxErr1 Then
Begin
MaxErr1:=Abs(p1-p3); VidnErr1:=Abs((p1-p3)/p1);
end;
If Abs(p1-p2)>MaxErr2 Then
Begin
MaxErr2:=Abs(p1-p2); VidnErr2:=Abs((p1-p2)/p1);
End;
End;
Assign(f,'mix.txt'); Append(f);
WriteLn(f,nx:4,ny:4,MaxErr2:19:12,VidnErr2:19:12,MaxErr1:19:12,VidnErr1:19:12);
Close(f);
End;
Begin
For cc:=1 To 10 Do For cc1:=1 To 5 Do
Begin
DataInput; cx:=33; cy:=33;
WriteLn('Nx=',nx,' Ny=',ny);
New(D1); New(D2); BuildCoefTable; Dispose(D1); Dispose(D2);
GetMaxError;
End;
WriteLn('Press <ENTER>'); ReadLn;
End.