1. Interpolācija

Report
Interpolēšanas metožu izmantošana
Delphi un Matlab vidē.
Praksē nereti rodas nepieciešamība noskaidrot
sakarības dažādos procesos un parādībās un
izteikt šīs sakarības matemātiskas izteiksmes
veidā.
Ir zināmi trīs funkciju uzdošanas veidi:
analītiskais, grafiskais un tabulārais.
Interpolācija visbiežāk jāizmanto tad, ja
funkcija ir uzdota tabulas veidā.
Funkciju interpolācija
y
x
1
2
…
i
i+1
…
n–1
n
Lineārā interpolācija
y  a*xb
- Taisnes vienādojums
yi  a * xi  b
y i 1  a * x i 1  b
Funkcijas vērtības
kaimiņ mezglu punktos
Atrodam a un b vērtības:
b  yi  a * xi
y i 1  a * x i 1  y i  a * x i
a 
y i 1  y i
x i 1  x i
b  yi  xi *
Nosacījums
y i 1  y i
x i 1  x i
xi < Xt < xi+1
Piemērs.
xt := strtofloat(edit1.text);
For i := 1 to n-1 do
begin
if (xt < x[i+1]) and (xt > x[i]) then
begin
a:= (y[i+1] – y[i])/(x[i+1] – x[i]);
b:= y[i] – x[i] * a;
yt:= a*xt + b;
end
end;
Kvadrātiskā interpolācija
a*x2 + b*x + c = 0;
Parabolas vienādojums
y[i-1] = x[i-1]*x[i-1]*a + x[i-1]*b + c
y[i] = x[i]*x[i]*a + x[i]*b + c
y[i+1] = x[i+1]*x[i+1]*a + x[i+1]*b + c
No vienādojumu sistēmas
atrodam a, b un c vērtības:
a = ∆a /∆; b = ∆b /∆; c = ∆c /∆;
kur ∆, ∆a, ∆b, ∆c – trešās kārtas
determinanti.
xi-1 <= x <= xi+1
Splain-interpolācija
m kārtības splain dēvējas funkcija,
kura ir m pakāpes polinoms uz katras
atstarpes starp mezglu punktiem un
visos iekšējos mezglos apmierina
funkcijas un atvasinājumu
kontinuitātes nosacījumiem.
Vislielāko izplatīšanu saņēma kuba
splains
y=ax3+bx2+cx+d.
Atradīsim atkarības viņa a, b, c, d
koeficentu aprēķināšanai pa funkcijas
un tās pirmo atvasinājumu vērtībām
kaimiņ mezglu punktos.
Izkravāšanu vienkāršojumam pieņemsim,
ka atstarpes kreisa leņķiska punkta xi,xi+1
argumenta vērtība, kurai meklējam splainfunkciju, vienāds nullei. To vienmēr var
izdarīt pārveidi x=xi+z, kur z jauns
arguments.
y=az3+bz2+cz+d
xi+1
z = x - xi
xi
z1 = 0
z2 = xi+1 - xi
Tad uzdevums ir noformulējams sekojoši:
atrast nezināmus a, b, c, d splaina koeficentus,
apmierinoši mezglu punktos nosacījumam
z1=0; y(z1)=f1; y'(z1)=f1'; y(z2)=f2; y' (z2)=f2'.
Šeit f1,f2,f1',f2' - attiecīgi interpolēšanas
funkcijas un tās atvasinājuma vērtības mezglu
punktos z1, z2.
Lai atrastu a,b,c,d saliksim vienādojumu
sistēmu, izmantojot nosacījumus mezglu
punktos
f1 = d;
f1′ = c;
f2 = az23 + bz22 + cz2 +d;
f2′ = 3az22 + 2bz +c; tātad
d = f1;
c = f1′
Paliekot atrastās d, c vērtības divos
pēdējos vienādojumos un atļaujot tos
relatīvi a, b saņemsim
a = (z2(f2′ - f1′ ) – 2(f2 - f1′ * z2 – f1))/z23
b = (3(f2 - f1′ * z2 – f1) – z2(f2′ - f1′ ))/z22
Paliekot atrastos koeficentus kuba
splainu un uzdodot argumenta vērtību,
aprēķināsim funkcijas vērtību starp
mezglu punktiem.
y = az3 + bz2 + cz + d
Bieži pie interpolācijas ir esamas tikai
funkcijas vērtības mezglu punktos un nav
atvasinājumu vērtību. Tad pirmie atvasinājumi
iekšējos mezglu punktos aprēķina pa
formulām
fi′ = (fi+1 – fi-1)/(xi+1 – xi-1),
fi′ = (fi+1 – fi-1)/(2h), i = 2, .. N-1
h- attālums starp mezglu punktiem.
N – mezglu punktu skaits
Atvasinājuma vērtību pirmajā un pēdējā
mezglu punktos aprēķina pa formulām:
f1′ = f2′ - f2″ h; fn′ = fn-1′ – fn-1″ h;
f2″ = (f1 + f3 – 2f2)/h2;
fn-1″ = (fn-2 + fn – 2fn-1)/h2;
f2″ , fn-1″ - otro atvasinājumu pietuvinātās
vērtības otrajā un n - 1 punktos.
Programmas fragmenti Delphi vidē
Funkcijas kubisks splains
function k_spl(xx:Real):Real;
var
d:Real; // koeficient d
{pirmā atvasinājuma aprēķins mezgla punktos
(i>= 2) and (i <= n-1)}
function fi_1(i:Byte):Real;
begin
fi_1 := (y[i+1] - y[i-1])/(x[i+1]-x[i-1]);
end;
{otrā atvasinājuma aprēķins punktā i = 2}
function f22:single;
begin
f22 := (y[1] + y[3] - 2*y[2])/sqr(x[2]-x[1]);
end;
{pirmā atvasinājuma aprēķins punktā i = 1}
function f12:single;
begin
f12 := fi_1(2) - f22*(x[2]-x[1]);
end;
{otrā atvasinājuma aprēķins punktā i = n-1}
function fn2:single;
begin
fn2 := (y[n-2] + y[n] - 2*y[n-1])/sqr(x[n]-x[n-1]);
end;
{pirmā atvasinājuma aprēķins punktā i = n}
function f1n:single;
begin
f1n := fi_1(n-1) - fn2*(x[n]-x[n-1]);
end;
{koeficienta c aprēķināšana}
function c(i:byte):single;
begin
if (i>= 2) and (i < n-1) then c := fi_1(i);
if i = 1 then c := f12;
if i = n-1 then c := f1n;
end;
{koeficienta a aprēķinšāna}
function a(i:Byte):Real;
begin
if (i>= 2) and (i < n-1) then
begin
a := ((x[i+1]-x[i])*(fi_1(i+1)-fi_1(i)) - 2*(y[i+1] - fi_1(i)*(x[i+1]-x[i])
- y[i]))/power(x[i+1]-x[i],3);
end;
if i= 1 then
begin
a := ((x[i+1]-x[i])*(fi_1(i+1)-fi_1(i)) - 2*(y[i+1] - f12*(x[i+1]-x[i])
- y[i]))/power(x[i+1]-x[i],3);
end;
if i= n-1 then
begin
a := ((x[i+1]-x[i])*(f1n-fi_1(i)) - 2*(y[i+1] - fi_1(i)*(x[i+1]-x[i])
- y[i]))/power(x[i+1]-x[i],3);
end;
end;
{koeficienta b aprēķināšana}
function b(i:Byte):Real;
begin
if (i>= 2) and (i < n-1) then
begin
b := (3*(y[i+1] - fi_1(i)*(x[i+1]-x[i]) - y[i])-(x[i+1]-x[i])*
(fi_1(i+1)-fi_1(i)))/sqr(x[i+1]-x[i]);
end;
if i = 1 then
begin
b := (3*(y[i+1] - f12*(x[i+1]-x[i]) - y[i])- (x[i+1]-x[i])*(fi_1(i+1)fi_1(i)))/sqr(x[i+1]-x[i]);
end;
if i= n-1 then
begin
b := (3*(y[i+1] - fi_1(i)*(x[i+1]-x[i]) - y[i])(x[i+1]-x[i])*(f1n-fi_1(i)))/sqr(x[i+1]-x[i]);
end;
end;
{Funkcijas ķermenis}
begin
for j := 1 to n-1 do
if(xx > x[j]) and (xx < x[j+1]) then it := j;
{it – kreisa mezgla numurs
x[it] - kreisa mezgla koordināta
xx – argumenta vērtība}
d := y[it];
k_spl := a(it)*power(xx-x[it],3) + b(it)*sqr(xx-x[it])+
c(it)*(xx-x[it]) + d;
end;
{Funkcijas y(xt) vērtības aprēķināšana}
procedure TForm1.Button1Click(Sender: TObject);
begin
xt := strtofloat(edit1.Text);
{ lasīt vektorus x un y no stringgrid komponenta}
for i := 1 to n do
begin
x[i] := strtofloat(stringgrid1.Cells[i-1,0]);
y[i] := strtofloat(stringgrid1.Cells[i-1,1]);
end;
yt := k_spl(xt) ;
{vērtības y(xt) izvade}
edit2.Text := floattostrf(yt,fffixed,10,3);
end;
{Grafika konstruēšana}
procedure TForm1.Button2Click(Sender: TObject);
begin
Series1.Clear;
Series2.Clear;
for i := 1 to n do
begin
x[i] := strtofloat(stringgrid1.Cells[i-1,0]);
y[i] := strtofloat(stringgrid1.Cells[i-1,1]);
Series1.AddXY(x[i],y[i]);
end;
xg := x[1];
yg := y[1];
Series2.AddXY(xg,yg);
repeat
xg := xg + 0.05; yg := k_spl(xg) ;
Series2.AddXY(xg,yf);
until xg >= x[n];
end;
{Lasīt vektorus x un y no faila, izmantojot OpenDialog komponentu}
procedure TForm1.Open1Click(Sender: TObject);
begin
if OpenDialog1.Execute then
begin
AssignFile(ff,OpenDialog1.FileName);
Reset(ff);
i:=0;
while not Eof(ff) do
begin
i := i+1;
readln(ff,x[i],y[i]);
end;
end;
CloseFile(ff);
end;
{Fails dati.txt Notepadā}
-5.0
-4.2
-3.4
-2.6
-1.8
-1.0
-0.2
0.6
1.4
2.2
9.13
1.26
-0.7
1.03
-2.05
-7.23
-6.15
-2.02
-3.02
-5.43

similar documents