Splines (Nelineární interpolace)

Tomáš Kalvoda, KAM FIT, BI-ZMA 2013/2014

Motto: Jak spojit množinu zadaných bodů ,,hezkou’’ křivkou?

Propojení bodů přímkami (linární interpolace)

Jsou-li dány body (a,b) a (c,d) v rovině takové, že a≠c, pak přímka procházející těmito body je popsána rovnicí

y=((d-b)/(c-a))(x-a)+b.

Příklad

Nechť jsou zadány body:

points={{1,2},{2,-1},{3,1},{4,-2}};

Pak rovnice přímek a jejich grafy lze napočíst takto:

lineq[A_,B_]:=((B[[2]]-A[[2]])/(B[[1]]-A[[1]]))(x-A[[1]])+A[[2]];
lines=lineq@@@Partition[points,2,1];
doms=Partition[points//Transpose//First,2,1];
range={{Min[#1]-.2,Max[#1]+.2},{Min[#2]-.2,Max[#2]+.2}}&@@Transpose[points];
plots=Table[
Plot[lines[[i]],{x,doms[[i,1]],doms[[i,2]]},
PlotRange->range,
AspectRatio->1
],
{i,1,Length[lines]}
];
Show[{plots,Graphics[{Red,Disk[#,.03]&/@points}]}]

Samozřejmě, Mathematica má pro toto vlastní funkce. Stejný příklad aniž bych explicitně konstruoval rovnice přímek:

Graphics[{Blue,Line[points],Red,Disk[#,.03]&/@points}]

Propojení bodů polynomy

Místo spojování bodů přímkami (grafy polynomů prvního stupně) mužeme body spojovat polynomiálními křivkami. Často se použivá polynomů třetího stupně (kubická interpolace, splines). Důvodem je, že mají poměrně jednoduchý tvar, avšak už umožňují inflexi. Kubický polynom je funkce tvaru

f(x)=a (x^3)+b (x^2)+c x +d.

Je tedy jednoznačně zadána čtyřmi konstantami a,b,c,d. Srovnejte s linární funkcí, k jejímuž popisu stačí dvě konstanty.

Uvažujme tedy pro jednoduchost množinu čtyř bodů

{((Subscript[x, 1]),(Subscript[y, 1])),((Subscript[x, 2]),(Subscript[y, 2])),((Subscript[x, 3]),(Subscript[y, 3])),((Subscript[x, 4]),(Subscript[y, 4]))}

Předpokládejme, že jsou uspořádány podle první složky, (Subscript[x, 1])<(Subscript[x, 2])<(Subscript[x, 3])<(Subscript[x, 4]). Kdyby náhodou dvě x-ové souřadnice byly stejné, pak tyto body v ..tomto směru’’ nelze spojit a je třeba změnit souřadnice.

Hledáme tři kubické polynomy

(Subscript[f, i])(x)=(Subscript[a, i])(x^3)+(Subscript[b, i])(x^2)+(Subscript[c, i])x+(Subscript[d, i]), x∈⟨(Subscript[x, i]),(Subscript[x, i+1])⟩,  i=1,2,3.

Požadujeme, aby zadané body spojovaly, tedy máme 6 podmínek na funkční hodnoty

(Subscript[f, 1])((Subscript[x, 1]))=(Subscript[y, 1]),(Subscript[f, 1])((Subscript[x, 2]))=(Subscript[f, 2])((Subscript[x, 2]))=(Subscript[y, 2]), (Subscript[f, 2])((Subscript[x, 3]))=(Subscript[f, 3])((Subscript[x, 3]))=(Subscript[y, 3]), (Subscript[f, 3])((Subscript[x, 4]))=(Subscript[y, 4]).

Dále požadujeme, aby ve spojích měly stejnou první i druhou derivaci (stejný sklon, tečna a stejná konkávnost/konvexnost), to nám dává další čtyři podmínky

(Subscript[f, 1]^('))((Subscript[x, 2]))=(Subscript[f, 2]^('))((Subscript[x, 2])), (Subscript[f, 2]^('))((Subscript[x, 3]))=(Subscript[f, 3]^('))((Subscript[x, 3])), (Subscript[f, 1]^(''))((Subscript[x, 2]))=(Subscript[f, 2]^(''))((Subscript[x, 2])), (Subscript[f, 2]^(''))((Subscript[x, 3]))=(Subscript[f, 3]^(''))((Subscript[x, 3])).

Konečně, můžeme zadat jaký požadujeme sklon na koncích, tedy dvě podmínky

(Subscript[f, 1]^('))((Subscript[x, 1]))=A, (Subscript[f, 3]^('))((Subscript[x, 4]))=B.

Celkem máme 12 rovnic pro 12 neznámých. Rovnice jsou lineární. Řešením získáme hodnoty 12-ti konstant a pak stačí vykreslit takto vypočtene polynomy na příslušných intervalech.

Příklad

Uvažme stejné body jako v prvním příkladě.

points={{1,2},{2,-1},{3,1},{4,-2}};

Navíc požadujme sklony na krajích.

slopes={1,0};

Následující funkce počítá polynomy, vstupem je uspořádaný (podle první složky) seznam bodů, sklony na koncích a booleovská proměnná umožňující zobrazit řešenou soustavu rovnic.

Splines[points_,slopes_,showeqs_]:=Module[{cubic,Dcubic,DDcubic,n=Length[points],eqs,vars,sols,out},
cubic[ind_,x_]=a[ind](x^3)+b[ind](x^2)+c[ind]x+d[ind];
Dcubic[ind_,x_]=D[cubic[ind,x],x];
DDcubic[ind_,x_]=D[cubic[ind,x],{x,2}];
(* values *)
eqs=Table[
{
cubic[i,points[[i,1]]] == points[[i,2]],
cubic[i,points[[i+1,1]]] == points[[i+1,2]]
},{i,1,n-1}
]//Flatten;
(* add border slopes *)
eqs=Join[eqs,
{
Dcubic[1,First[points][[1]]]==slopes[[1]],
Dcubic[n-1,Last[points][[1]]]==slopes[[2]]
}
];
(* first and second derivatives *)
eqs=Join[eqs,
Table[
{
Dcubic[i,points[[i+1,1]]]==Dcubic[i+1,points[[i+1,1]]],
DDcubic[i,points[[i+1,1]]]==DDcubic[i+1,points[[i+1,1]]]
},{i,1,n-2}
]//Flatten
];
(* if required, print equations *)
If[showeqs,Print[Grid[{Range[Length[eqs]],eqs}//Transpose,Frame->All]]];
(* unknowns *)
vars=Table[{a[i],b[i],c[i],d[i]},{i,1,n-1}]//Flatten;
(* solve the system *)
sols=Solve[eqs,vars];
(* return the result *)
out=Table[cubic[i,x]/.sols,{i,1,n-1}];
Return[out];
]

V našem příkladě.

splines=Splines[points,slopes,True]

Případně bez zlomků, ale se zaokrouhlovací chybou:

%//N

Vykresleme výsledek.

splines=Splines[points,slopes,False]//N;
doms=Partition[points//Transpose//First,2,1];
range={{Min[#1]-.2,Max[#1]+.2},{Min[#2]-.2,Max[#2]+.2}}&@@Transpose[points];
color[i_]:=Which[OddQ[i],Blue,True,Red];
imgSpl=Table[
Plot[splines[[i]],{x,doms[[i,1]],doms[[i,2]]},
PlotRange->range,
PlotStyle->{color[i],Thick},
AspectRatio->1
],{i,1,doms//Length}
];
Show[{
imgSpl,Graphics[{Black,Opacity[.9],Disk[#,.03]&/@points}]
}]

Interaktivní příklad

color[i_]:=Which[OddQ[i],Blue,True,Red];
DynamicModule[{pt={{2,.5},{3,3},{5,2},{6,3}},A=0,B=0},
Column[{
Row[{"Levý sklon: ",Slider[Dynamic[A],{-10,10,.1},Appearance->"Labeled"]}],
Row[{"Pravý pravý: ",Slider[Dynamic[B],{-10,10,.1},Appearance->"Labeled"]}],
LocatorPane[Dynamic[pt],
Dynamic@Block[{spl,dom,img,aux},
aux=Sort[pt,#1[[1]]<#2[[1]]&];
spl=Splines[aux,{A,B},False];
dom=Table[{aux[[i,1]],aux[[i+1,1]]},{i,1,Length[aux]-1}];
img=Join[{Table[
Plot[spl[[i]],{x,dom[[i,1]],dom[[i,2]]},
PlotRange->{{0,8},{0,4}},
PlotStyle->{color[i]},
ImageSize->600
],{i,1,dom//Length}
]},
{Plot[{A(x-aux[[1,1]])+aux[[1,2]],B(x-Last[aux][[1]])+Last[aux][[2]]},{x,0,8},PlotRange->{{0,8},{0,4}},PlotStyle->{Gray}]}
];
Show[img]
]
]
},Center]
]