Num. Math.—10

Lagrange Interpoláció, polinomfejlesztési algoritmus

A polinomot meghatárohatjuk a következô eljárással (minden egyes lineáris egyenletben pontosan egy ismeretlen (A)):

Az apriori alak: (Subscript[p, j+1])=(Subscript[p, j])+A ω

X = {1, 2, 3, 4, 5}; Y = {0, 6, 24, 60, 120};

p0=0;

p1=p0+A(x-1)

p1/.x->2

Solve[(p1/.x->2)==6]

p1=p0+6(x-1)

p2=p1+A(x-1)(x-2)

p2/.x->3

Solve[(p2/.x->3)==24]

p2=p1+6(x-1)(x-2)

p3=p2+A(x-1)(x-2)(x-3)

Solve[(p3/.x->4)==60]

p3=p2+1 (x-1)(x-2)(x-3)

p4=p3+A(x-1)(x-2)(x-3)(x-4)

Solve[(p4/.x->5)==120,A]

p4=p3+0(x-1)(x-2)(x-3)(x-4)

Expand[%]

X={1,2,3};
F={2,5,10};

Clear[LP];
LP[0]:=(P=F[[1]]);

LP[n_]:=Module[{pl},pl=P+A Product[x-X[[j]],{j,1,n}];
P=pl/.Solve[(pl/.x->X[[n+1]])==F[[n+1]],A][[1]];{pl,(pl/.x->X[[n+1]])==F[[n+1]],Solve[(pl/.x->X[[n+1]])==F[[n+1]],A][[1]],P,Expand[P]}]

Manipulate[TableForm[{{LP[n]}}],{{n,0},0,2,1}]

Newton előállítása osztott differenciákkal (rekurzív)

A Lagrange alaknál új adat hozzávételére nincs lehetôség ⟶ Newton féle elôállítás osztott differenciák segítségével

Clear[OD,x0,y0,x1,y1]

NewtInterp[x_List,y_List,var_]:=(∑_(j=1)^(Length[x]))OD[Take[x,j],Take[y,j]](∏_(i=1)^(j-1))(var-x[[i]])

NewtInterp[{x0,x1,x2},{y0,y1,y2},x]

OD[{x_},{y_}]:=y;
OD[x_List,y_List]:=(OD[Drop[x,1],Drop[y,1]]-OD[Drop[x,-1],Drop[y,-1]])/(Last[x]-First[x])

OD[{x0,x1},{y0,y1}]

NewtInterp[{x0,x1,x2},{y0,y1,y2},x]

Konkrét Lagrange polinom előállításhoz:

NewtInterp[{1,2,3,4,5},{0,6,24,60,120},x]

Expand[%]

Megjegyzés. Séma az f[x0], f[x0;x1], f[x0;x1;x2], f[x0;x1;x2;x3],f[x0;x1;x2;x3;x4] osztott differencia sorozat kiszámításához.

(Subscript[L, 4])[x]=((((0)+6 (-1+x))+6 (-2+x) (-1+x))+1(-3+x) (-2+x) (-1+x))+0(-1+x)(-2+x)(-3+x)(-4+x)(-5+x)


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

1    0  
               6   
2    6               6
             18            1
3   24             9            0
             36            1
4   60             12
             60
5  120

Hermite: Black-Box Feladatmegoldás

Mi a következõ Hermite interpolációs probléma megoldása?

x0={1,2,3,4}
y0={1,2,4,3}
yd0={1,0,0,-1}

Ábrázoljuk az "adatokkal" a polinomot!
Oldjuk meg a feladatot a Mathematica függvénnyel!

Built-in Mathematica függvény

?InterpolatingPolynomial

InterpolatingPolynomial[{{1,{-1,1}},{3,{21,25}}},x]//Expand

x0={1,2,3,4};
y0={1,2,4,3};
yd0={1,0,0,-1};

hp=Expand[InterpolatingPolynomial[Transpose[{x0,Transpose[{y0,yd0}]}],x]]

Grafikon

Plot[hp,{x,0,5},PlotRange->{-1,5},PlotStyle->Thickness[.005]]

Grafikon+pontok

L[{xx_,yy_,ydd_},var_]:=ydd(var-xx)+yy

L[{1,4,2},x]//Expand

ls=Map[L[#,x]&,Transpose[{x0,y0,yd0}]]

grafikon+pontok+érintőegyenesek

Plot[Evaluate[Prepend[ls,hp]],{x,0,5},PlotRange->{-1,5},PlotStyle->Prepend[Table[{Blue},{4}],{Thickness[.005],Black}],Epilog->{Red,PointSize[.02],Map[Point[#]&,Transpose[{x0,y0}]]}]

grafikon+pontok+érintőegyenes-szakaszok

Plot[hp,{x,0,5},PlotStyle->Black,Epilog->{RGBColor[1,0,0],PointSize[.02],Map[Point[#]&,Transpose[{x0,y0}]],RGBColor[0,0,1],Map[Line[{{#[[1]]-.5,#[[2]]-#[[3]].5},{#[[1]]+.5,#[[2]]+#[[3]].5}}]&,Transpose[{x0,y0,yd0}]]}]

HDEMO

Manipulate[Grid[{{Panel[Grid[{Prepend[DF,"DF:"],Prepend[F,"F:"]}]],GraphPlot[MM[n],PlotRange->{{1-.2,Length[X]+.2},{-.3,1.3}},VertexCoordinateRules->Join[Table[j->{j,0},{j,Length[X]}],Table[Length[X]+j->{j,1},{j,Length[X]}]],ImageSize->{200,100},VertexRenderingFunction-> ({GM[#2],If[ChFLists[IntegerPart[#2]],If[#2==GInd,Circle[#1,.15],Circle[#1,.1]],If[#2==GInd,Disk[#1,.15],Disk[#1,.1]]]}&)],Panel[Grid[{MapIndexed[Button[#,n=
f[
#2[[1]],n]]&,ll]}]]},{x,TableForm[{X}],GHP}},Frame->All],{{n,1},1,100,.1,Enabled->False},ControlPlacement->Right, FrameLabel->Style["Hermite Interpolation 2012",Red],SaveDefinitions->True]

Hermite Demo aux

GM[j_]:=If[j>Length[X],Red,Blue]

ChFLists[j_]:=If[j>Length[X],If[MemberQ[MDF[[All,1]],j],True,False],If[MemberQ[MF[[All,1]],j],True,False]]

ll={"⇑","⇓","⇐","⇒","+","I"};

f[6,s_]:=(X={1,3};F={2,0};DF={0,1};GInd=1;PGInd=-1;GHP=0;
M=IdentityMatrix[2Length[X]];
MF=Transpose[{Range[Length[X]],F}];CMF={};
MDF=Transpose[{Range[Length[X]]+Length[X],DF}]);

f[i_,s_]:=Module[{},If[i==
5,If[1<=GInd≤Length[X],GHP=GHP+A Product[If[j==GInd,1,(x-X[[j]])^(2-Count[MF[[All,1]],j]-Count[MDF[[All,1]],j+Length[X]])],{j,Length[X]}];GHP=GHP/.(Solve[(GHP/.x->X[[GInd]])==F[[GInd]],A][[1]]);MF=DeleteCases[MF,{GInd,_?NumberQ}];If[PGInd>0,M[[PGInd,GInd]]=1];PGInd=GInd,If[Length[X]+1<=GInd≤2Length[X]∧Cases[MF,{GInd-Length[X],_}]=={},GHP=GHP+A Product[If[j==GInd-Length[X],x-X[[j]],(x-X[[j]])^(2-Count[MF[[All,1]],j]-Count[MDF[[All,1]],j+Length[X]])],{j,Length[X]}];GHP=GHP/.(Solve[(D[GHP,x]/.x->X[[GInd-Length[X]]])==DF[[GInd-Length[X]]],A][[1]]);MDF=DeleteCases[MDF,{GInd,_?NumberQ}];If[PGInd>0,M[[PGInd,GInd]]=1];PGInd=GInd]],If[i==4∧(1<=GInd<Length[X]∨Length[X]<GInd<2Length[X]),GInd++,If[i==3∧(1<GInd≤Length[X]∨Length[X]+1<GInd<=2Length[X]),GInd--,If[i==1∧1≤GInd≤Length[X],GInd=GInd+Length[X],If[i==2∧Length[X]<GInd≤2Length[X],GInd=GInd-Length[X]]]]]];s+1]

MM[n_]:=M

Példa

A Hermite interpolációs polinomot ismét szukcesszív lineáris egyenletmegoldással keressük meg, lásd a Lagrange-problémánál lévð algoritmust!

1. Példa p(1)=2, p'(1)=0, p(3)=0, p'(3)=1

azaz X={1,3}
          Y={2,0}
       YD={0,1}

InterpolatingPolynomial[{{1,{2,0}},{3,{0,1}}},x]

InterpolatingPolynomial[{{1,{2,0}},{3,{0,1}}},x]//Expand

Poliomfejlesztés és a 
A ND-s háromszögtáblázat  értelemszerű módosítása

Clear[p0,p1,p2,p3]

p0=2;

p1=p0+A(x-1)

Solve[(D[p1,x]/.x->1)==0,A]

p1=2+0(x-1)//Expand

p2=p1+A(x-1)(x-1)

Solve[(p2/.x->3)==0,A]

p2=p1+(-1/2)(x-1)^2

p3=p2+A(x-1)^2(x-3)

Solve[(D[p3,x]/.x->3)==1,A]

p3=p2+3/4(x-1)^2(x-3)

(Subscript[x, i]) (Subscript[f, i  ])          1 .      2.     3.

1    2  
               0   
1    2              -1/2
             -1                 3/4
3   0              1           
              1          
3   0

Hermite alappolinomok (a Lagrange alappolinomok felhasználásával)

Az interpolációs probléma megoldása egy legfeljebb (2n-1)-edfokú polinom. Pl. n=3, deg P≤5
Az elôállítás alapgondolata: Kétféle bázispolinom: (Subscript[h, j])ill. (Subscript[h, n+j]) alakjuk:  ((Subscript[a, j])x+(Subscript[b, j]))(((Subscript[l, j])[x])^( 2)) ill. (x-(Subscript[x, j]))(((Subscript[l, j])[x])^( 2))

Mindkettô támaszkodik a Lagrange-féle (Subscript[l, j]) bázispolinomokra. A feltételek (Subscript[h, j])[(Subscript[x, k])]=(Subscript[δ, j k]), (Subscript[h, j])'[(Subscript[x, k])]=0, (Subscript[h, n+j])[(Subscript[x, k])]=0,(Subscript[h, n+j])' [(Subscript[x, k])]=(Subscript[δ, j k])

6 feltétel automatikusan teljesül, 2 pedig meghatározza a lineáris faktor együtthatóit (1-2(x-(Subscript[x, j]))(Subscript[l, j])' [x])

LagrBase[j_,x_List,var_]:=((∏_(k=1)^j)((var-x[[k]])/(x[[j+1]]-x[[k]])))((∏_(k=j+2)^(Length[x]))((var-x[[k]])/(x[[j+1]]-x[[k]])))

Clear[HermBase]

HermBase[j_,x_List,var_]/;j<Length[x]:=(1-2(var-x[[j+1]])((D[LagrBase[j,x,var],var])/.var->x[[j+1]]))LagrBase[j,x,var]^2;

HermBase[j_,x_List,var_]:=(var-x[[j-Length[x]+1]])LagrBase[j-Length[x],x,var]^2;

HermInterp[x_List,y_List,yd_List,var_]:=(∑_(j=0)^(Length[x]-1))y[[j+1]]HermBase[j,x,var]+(∑_(j=0)^(Length[x]-1))yd[[j+1]]HermBase[Length[x]+j,x,var]

Feladat

x={x0,x1,x2}={1,2,3}
y={y0,y1,y2}={2,32,242}
yd={yd0,yd1,yd2}={4,79,404}

Adjuk meg a Hermite alappolinomokat és a Hermite interpolációs polinomot!
Ellenõrizzük a nevezetes tulajdonságokat!

Javaslat

x0={1,2,3};
y0={2,32,242};
yd={4,79,404};

HermBase[3,x0,x]

LagrBase[0,x0,x]

HermBase[0,x0,x]

q=(a x+b)LagrBase[0,x0,x]^2

q/.x->x0[[1]]

D[q,x]/.x->x0[[1]]

q

Solve[{(q/.x->x0[[1]])==1,(D[q,x]/.x->x0[[1]])==0},{a,b}]

q/.Solve[{(q/.x->x0[[1]])==1,(D[q,x]/.x->x0[[1]])==0},{a,b}][[1]]//Expand

HermBase[0,x0,x]//Expand

H=Table[HermBase[j,x0,x],{j,0,2 Length[x0]-1}]

Length[%]

A nulladik és a harmadik bázispolinomra kirótt feltételek vizsgálata

h0[x_]=Expand[HermBase[0,x0,x]];

h3[x_]=Expand[HermBase[3,x0,x]];

h3[x0[[3]]]

hd0[x_]=D[HermBase[0,x0,x],x];

hd3[x_]=D[HermBase[3,x0,x],x];

{{h0[1],h0[2],h0[3]},{hd0[1],hd0[2],hd0[3]}}

{{h3[1],h3[2],h3[3]},{hd3[1],hd3[2],hd3[3]}}

Az elsô és a negyedik bázispolinomra kirótt feltételek vizsgálata

h1[x_]=Expand[HermBase[1,x0,x]];

h4[x_]=Expand[HermBase[4,x0,x]];

hd1[x_]=D[HermBase[1,x0,x],x];

hd4[x_]=D[HermBase[4,x0,x],x];

{{h3[1],h3[2],h3[3]},{hd3[1],hd3[2],hd3[3]}}

{{h1[1],h1[2],h1[3]},{hd1[1],hd1[2],hd1[3]}}

{{h4[1],h4[2],h4[3]},{hd4[1],hd4[2],hd4[3]}}

h2[x_]=Expand[HermBase[2,x0,x]];

h5[x_]=Expand[HermBase[5,x0,x]];

hd2[x_]=D[HermBase[2,x0,x],x];

hd5[x_]=D[HermBase[5,x0,x],x];

{{h2[1],h2[2],h2[3]},{hd2[1],hd2[2],hd2[3]}}

Hermite Interpolációs Polinom

HermInterp[x0,y0,yd,x]

Fokszám: 5

Expand[HermInterp[x0,y0,yd,x]]

p1=2-x+(x^5);

D[p1,x]/.x->3

x0

y0

Clear[p2]

p2[x_]:=2-x+(x^5);

p2[3]

*CRM

Base Algorithm: (Newton)

Consider only the first two and solve it

r=r1’+σm1 (m1m2)
σ=(r2-r1’) c (m2)
c=m1^(-1) (m2)

X={1,2,3};Y={2,5,10};

PolynomialExtendedGCD[(x-1),x-2,x]

c=1;r1=2;r2=5;σp=(r2-(r1))c

σ=PolynomialRemainder[%,x-2,x]

r1+σ (x-1)//Expand

PolynomialExtendedGCD[(x-1)(x-2),x-3,x]

{(1/2),-(x/2)}.{(x-1)(x-2),x-3}//Expand

c=1/2;r1=3x-1;r2=5;σp=(10-(r1))c

σ=PolynomialRemainder[%,x-3,x]

r1+σ (x-1)(x-2)//Expand

ALSO works with Hermite!

InterpolatingPolynomial[{{1,{2,0}},{3,{0,1}}},x]//Expand

PolynomialExtendedGCD[(x-1)^2,(x-3)^2,x]

(1-x/4) (x-1)^2+x/4 (x-3)^2

Expand[%]

p=2 (x/4) (x-3)^2+(x-3) (1-x/4) (x-1)^2//Expand

pp=PolynomialRemainder[p,(x-1)^2 (x-3)^2,x]

PolynomialRemainder[pp,(x-1)^2,x]

PolynomialRemainder[pp,(x-3)^2,x]

InterpolatingPolynomial[{{1,{2,4}},{2,{32,79}},{3,{242,404}}},x]//Expand

PolynomialExtendedGCD[(x-1)^2,(x-2)^2,x]

p1=(2+4(x-1)) (2x-1) (x-2)^2+(32+79(x-2)) (5-2x) (x-1)^2//Expand

pp1=PolynomialRemainder[p1,(x-1)^2(x-2)^2,x]

PolynomialExtendedGCD[(x-1)^2(x-2)^2,(x-3)^2,x]

p2=pp1 1/16(-16+48 x-40 (x^2)+12 (x^3)) (x-3)^2+(242+404(x-3)) 1/16(40-12x) (x-1)^2(x-2)^2//Expand

PolynomialRemainder[p2,(x-1)^2(x-2)^2(x-3)^2,x]