> #==================================

> # Uebung 9 Aufgabe 2: Interpolation

> #=================================

> #

> restart: with(linalg):

Warning, new definition for norm

Warning, new definition for trace

> nm:=5; x:=vector([0,1,2,4, -1, 3]); f:=vector([-3,1,2,7,1,6]);

> # Lagrange

> n:=3; p(x):=sum(f[jj]*( product((x-x[ii])/(x[jj]-x[ii]),ii=1..jj-1) * product((x-x[ii])/(x[jj]-x[ii]),ii=jj+1..n+1) ), jj=1..n+1);

> simplify(p(x)); # umwandeln in Normalform

> # Ueber lin Gleichungssytem

> A:=matrix(n+1,n+1): b:=vector(n+1):

> for id from 1 to n+1 do A[id,1]:=1; for jd from 2 to n+1 do A[id,jd]:=x[id]^(jd-1): od: b[id]:=f[id]: od: A:=evalm(A); b:=evalm(b);

> a:=augment(A,b); a:=gausselim(a); koeff:=backsub(a);

> lgsp(x):=sum(koeff[ii]*x^(ii-1), ii=1..n+1);

> # Newton'sche Berechnung

> NT:=matrix(nm+1,nm+1): for id from 1 to n+1 do NT[id,1]:=f[id]: od:

> for jd from 1 to n do for id from 1 to n+1 - jd do NT[id,jd+1]:=(NT[id+1,jd]-NT[id,jd])/(x[id+jd]-x[id]): od:od:

> for id from 1 to n+1 do printf(`%8.5f | `,x[id]); for jd from 1 to n+2-id do printf(`%8.5f `, NT[id,jd]); od; printf(`\n`); od;

0.00000 | -3.00000 4.00000 -1.50000 .50000

1.00000 | 1.00000 1.00000 .50000

2.00000 | 2.00000 2.50000

4.00000 | 7.00000

> np(x):=NT[1,n+1]: for id from n to 1 by -1 do np(x):=np(x)*(x-x[id])+NT[1,id]: od:

> np(x);

> simplify(np(x)); # Normalform zum Vergleich

> # Plots

> ll:=[[x[ii],f[ii]] $ii=1..n+1];

> plot({ll,[x,np(x), x=0..4]});

> # nun eine weitere Stuetzstelle

> # Lagrange: im Prinzip Alles neu

> n:=4; p(x):=sum(f[jj]*( product((x-x[ii])/(x[jj]-x[ii]),ii=1..jj-1) * product((x-x[ii])/(x[jj]-x[ii]),ii=jj+1..n+1) ), jj=1..n+1);

> simplify(p(x));

> # Ueber lin Gleichungssytem - auch Alles neu

> A:=matrix(n+1,n+1): b:=vector(n+1):

> for id from 1 to n+1 do A[id,1]:=1; for jd from 2 to n+1 do A[id,jd]:=x[id]^(jd-1): od: b[id]:=f[id]: od: A:=evalm(A); b:=evalm(b);

> a:=augment(A,b); a:=gausselim(a); koeff:=backsub(a);

> lgsp(x):=sum(koeff[ii]*x^(ii-1), ii=1..n);

> # Newton: Nur eine weitere Schraegzeile

> NT[n+1,1]:=f[n+1]: for jd from 1 to n do NT[n+1-jd,jd+1]:=(NT[n+2-jd,jd]-NT[n+1-jd,jd])/(x[n+1]-x[n+1-jd]): od:

> for id from 1 to n+1 do printf(`%8.5f | `,x[id]); for jd from 1 to n+2-id do printf(`%8.5f `, NT[id,jd]); od; printf(`\n`); od;

0.00000 | -3.00000 4.00000 -1.50000 .50000 .46667

1.00000 | 1.00000 1.00000 .50000 .03333

2.00000 | 2.00000 2.50000 .43333

4.00000 | 7.00000 1.20000

-1.00000 | 1.00000

> np(x):=NT[1,n+1]: for id from n to 1 by -1 do np(x):=np(x)*(x-x[id])+NT[1,id]: od:

> np(x); simplify(np(x));

> xh:=vector([-1,0,1,2,4]): fh:=vector([1,-3,1,2,7]): ll:=[[xh[ii],fh[ii]] $ii=1..n+1];

> plot({ll,[x,np(x), x=-1..4]});

> # eine weitere Stuetzstelle - noch eine Schraegzeile

> n:=5;NT[n+1,1]:=f[n+1]: for jd from 1 to n do NT[n+1-jd,jd+1]:=(NT[n+2-jd,jd]-NT[n+1-jd,jd])/(x[n+1]-x[n+1-jd]): od:

> for id from 1 to n+1 do printf(`%8.5f | `,x[id]); for jd from 1 to n+2-id do printf(`%8.5f `, NT[id,jd]); od; printf(`\n`); od;

0.00000 | -3.00000 4.00000 -1.50000 .50000 .46667 -.24167

1.00000 | 1.00000 1.00000 .50000 .03333 -.25833

2.00000 | 2.00000 2.50000 .43333 -.48333

4.00000 | 7.00000 1.20000 -.05000

-1.00000 | 1.00000 1.25000

3.00000 | 6.00000

> np(x):=NT[1,n+1]: for id from n to 1 by -1 do np(x):=np(x)*(x-x[id])+NT[1,id]: od:

> np(x); simplify(np(x));

> xh:=vector([-1,0,1,2,3,4]): fh:=vector([1,-3,1,2,6,7]): ll:=[[xh[ii],fh[ii]] $ii=1..n+1];

> plot({ll,[x,np(x), x=-1..4]});

>