Num. Math.—9

Warm-up: Lagrange interpoláció és lineáris algebra

Lineáris algebra.
Mutassuk meg a példánkon (X={1,2,3,4,5}, Y={0,6,24,60,120}), hogy a Lagrange interpoláció kezelhető a határozatlan együtthatók módszerével.

p=α0+α1 x+α2 (x^2)+α3 (x^3)+α4 (x^4);

Az (Subscript[α, j])együtthatókat az interpolációs krit alapján számoljuk ki.

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

es=Table[(p/.x->X[[j]])==Y[[j]],{j,1,5}]

Solve[es]

InterpolatingPolynomial[Transpose[{X,Y}],x]//Expand

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

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

InterpolatingPolynomial[Transpose[{{1,2,3},{2,5,10}}],x]

A Mathematica InterpolatingPolynomial függvény is így dolgozik.

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

Illusztráljuk a példánkon (X={1,2,3,4,5}, Y={0,6,24,60,120}), hogy az osztott differenciák tkp. egy alkalmas másik bázisban vett koordinátasorral egyeznek meg.
Megoldás. Legf. 4-edfokú polinomok lin tere (Subscript[P, 4]). Az intterpolációs feltϵtelrendszer ekvivalens egy Vandermonde együtthatómátrixú lineáris egyenletrendszerrel.

V[l_]:=Table[((l[[i]])^j),{i,Length[l]},{j,0,Length[l]-1}]

V[{1,2,3,4,5}]//MatrixForm

U={u0,u1,u2,u3,u4};

V[{1,2,3,4,5}][[2]].U

Inverse[V[{1,2,3,4,5}]].{0,6,24,60,120}

LinearSolve[V[{1,2,3,4,5}],{0,6,24,60,120}]

%.Table[(x^(j-1)),{j,5}]

A másik bázis elemei, bázisátmenetmátrix (BTM)

B ={1,(x-1),(x-1)(x-2),(x-1)(x-2)(x-3),(x-1)(x-2)(x-3)(x-4)}

(BTM=Transpose[Table[PadRight[CoefficientList[Expand[B[[j]]],x],5],{j,5}]])//TableForm

Inverse[BTM].Inverse[V[{1,2,3,4,5}]].{0,6,24,60,120}

Tehát:

%.B

Ekvidisztáns alappontrendszer, progresszív differenciák

Clear[x0,x1,y0,y1]

OD[{x0,x0+h},{f0,f1}]

OD[{x0,x0+h,x0+2h},{f0,f1,f2}]//Simplify

OD[{x0,x0+h,x0+2h,x0+3h},{f0,f1,f2,f3}]//Simplify

Progresszív differenciák

Clear[PD];

PD[{y0_}]:=y0;
PD[y_List]:=PD[Drop[y,1]-Drop[y,-1]]

F0={f0,f1,f2,f3};

Table[PD[Take[F0,j]],{j,Length[F0]}]

Nem kell (h^(n ))-nel osztani. Indoklás:

Legyen  a vizsgált pont x:=x0+th és nézzük pl. a kvadratikus tagot:

(x-x0)(x-x1)=t(t-1)(h^2) ⟹ f[x0;x1;x2](x-x0)(x-x1) (f0-2f1+f2)/(2(h^2)) (h^2)t(t-1)=PD[f0,f1,f2](({{t}, {2}}))

Csak az ált. bin eüh és a progreszív differencia kell.

Sin[π / 9] közelítése interpolációval. Lépésköx: π/12. x0=0

F0=N[Sin[{0,Pi/12,Pi/6,Pi/4,Pi/3}],4];

F0

Table[PD[Take[F0,j]],{j,1,Length[F0]}]

Interpolációs formula progresszív differenciákkal és ált. binomiális együtthatókkal

Newt2Interp[x0_,t_,DL_List]:=(∑_(j=0)^(Length[DL]-1))Binomial[t,j]DL[[j+1]]

Newt2Interp[0,4/3,Table[PD[Take[F0,j]],{j,1,Length[F0]}]]

Sin[Pi/9]//N

Binomiális együttható általánosítása

Binomial[4/3,3]

(4/3)(1/3)(-2/3)/6

A formula az eredeti általános interpolációs formulának valóban speciális esete

Clear[x0]

NewtInterp[{x0,x0+h,x0+2h,x0+3h,x0+4h},{f0,f1,f2,f3,f4},x]/.{x-> x0+t h}/.{t->4/3}//Simplify

(∑_(j=0)^4)Binomial[4/3,j]PD[Take[{f0,f1,f2,f3,f4},j+1]]//Simplify

*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]