Interpolaatio
Luento ke 31.1.2001
# read("/p/edu/mat-1.414/maple/v201.mpl"); # tai
# Ota www-sivulta ...maple/v201.mpl ja sitten
# read("v201.mpl"):
read("/mathwork/apiola/opetus/peruskurssi/v2-3/201/maple/v201.mpl");
M\357\277\275\357\277\275ritell\357\277\275\357\277\275n kuitenkin my\357\277\275s t\357\277\275ss\357\277\275
L:=proc(j,xd,x)
local oso,nimi,i,j1;
j1:=j+1;
oso:=product((x-xd[i]),i=1..nops(xd))/(x-xd[j1]);
nimi:=subs(x=xd[j1],oso);
oso/nimi;
end:
L(0,[a,b,c],x);
seq(L(j,xd,x),j=0..nops(xd)-1);
xd:=[seq(X[j],j=0..3)];
seq(L(j,xd,x),j=0..nops(xd)-1);
add(y[j]*L(j,xd,x),j=0..nops(xd)-1);
Nyt m\357\277\275\357\277\275rittelemme interpolaatiofunktion:
lagint:=(xd,yd,x)->add(yd[j+1]*L(j,xd,x),j=0..nops(xd)-1);
lagint([a,b,c],[ya,yb,yc],x); xd:=[a,b,c]: lagint(xd,map(f,xd),x);
p:=unapply(%,x);
p(a),p(b),p(c);
Maplen oma interp:
interp(xd,map(f,xd),x);
interR:=(f,xd,x)->(D@@nops(xd))(f)(xi)/(nops(xd))!*product((x-xd[i]),i=1..\134
nops(xd)):
interR(f,[a,b,c],x);
<Text-field style="Heading 1" layout="Heading 1">Rungen koe</Text-field>
with(plots):
f:=x->1/(1+25*x^2);x:='x':n:=5: xd:=linspace(-1,1,n):yd:=map(f,xd);
p:=interp(xd,yd,x);zip((x,y)->[x,y],xd,yd);
datakuva:=plot(zip((x,y)->[x,y],xd,yd),style=point,symbol=circle):fjap:=plot([f(x),p],x=-1..1,color=[brown,blue]):
display([datakuva,fjap]);
n:=10: xd:=linspace(-1,1,n):yd:=map(f,xd);
p:=interp(xd,yd,x);zip((x,y)->[x,y],xd,yd);
datakuva:=plot(zip((x,y)->[x,y],xd,yd),style=point,symbol=circle):fjap:=plot([f(x),p],x=-1..1,color=[brown,blue]):display([datakuva,fjap]);
n:=15: xd:=linspace(-1,1,n):yd:=map(f,xd): p:=interp(xd,yd,x):zip((x,y)->[x,y],xd,yd):datakuva:=plot(zip((x,y)->[x,y],xd,yd),style=point,symbol=circle):fjap:=plot([f(x),p],x=-1..1,color=[brown,blue]):display([datakuva,fjap]);
n:=15: xd:=linspace(-1,1,n):yd:=map(f,xd): p:=lagint(xd,yd,x):zip((x,y)->[x,y],xd,yd):datakuva:=plot(zip((x,y)->[x,y],xd,yd),style=point,symbol=circle):fjap:=plot([f(x),p],x=-1..1,color=[brown,blue]):display([datakuva,fjap]);
<Text-field style="Heading 1" layout="Heading 1">Splini</Text-field>
restart: read("/home/apiola/opetus/peruskurssi/v2-3/201/maple/v201.mpl");
s1:=x->1+a1*(x-1)+b1*(x-1)^2+c1*(x-1)^3;s2:=x->1/2+a2*(x-2)+b2*(x-2)^2+c2*(x-2)^3;s3:=x->1/3+a3*(x-3)+b3*(x-3)^2+c3*(x-3)^3;
Interpolaatioehdot oikeissa p\357\277\275\357\277\275tepisteiss\357\277\275 (samalla jatkuvuus sis\357\277\275solmuissa)
ie:=s1(2)=1/2,s2(3)=1/3,s3(4)=1/4;
1. derivaatat jatkuvat sis\357\277\275solmuissa:
d1:=D(s1)(2)=D(s2)(2),D(s2)(3)=D(s3)(3);
2. derivaatat jatkuvat sis\357\277\275solmuissa:
d2:=(D@D)(s1)(2)=(D@D)(s2)(2),(D@D)(s2)(3)=(D@D)(s3)(3);
Reunaehdot: Toiset derivaatat nollia reunoilla. Luonnolliset eli vapaat reunaehdot. (Pitk\357\277\275 ohut taipuisa vapa tai sauva, joka pakotetaan kulkemaan datapisteiden kautta, asettuu n\357\277\275in.
T\357\277\275ss\357\277\275 voi vaihtaa my\357\277\275s muunlaiset reunaehdot.
re:=(D@D)(s1)(1)=0,(D@D)(s3)(4)=0;
solve({ie,d1,d2,re},{a1,b1,c1,a2,b2,c2,a3,b3,c3});
assign(%);
s:=piecewise(x>=1 and x<2, s1(x),x>=2 and x<3, s2(x),x>=3 and x<=4, s3(x));
skuva:=plot(s,x=1..4):
Maplen oma spline:
ss:=spline([1,2,3,4],[1,1/2,1/3,1/4],x,cubic);
sskuva:=plot(.01+ss,x=1..4,color=blue): # Nostetaan sinist\357\277\275 m\357\277\275\357\277\275r\357\277\275ll\357\277\275 0.01, jotta erottuvat.
plots[display]([skuva,sskuva]);
plot(s-ss,x=1..4);
Taas t\357\277\275m\357\277\275 skaalausilmi\357\277\275, 1:1 antaa pelk\357\277\275n viivan.
Muuten, olipa nautinnollista tehd\357\277\275 Maplella!