// Lagrange and Hermite interpolation in Asymptote
// Author: Olivier Guibé

import interpolate;
import graph;

// Test 1: The Runge effect in the Lagrange interpolation of 1/(x^2+1).

unitsize(2cm);

real f(real x) {return(1/(x^2+1));}
real df(real x) {return(-2*x/(x^2+1)^2);}

real a=-5, b=5;
int n=15;
real[] x,y,dy;
x=a+(b-a)*sequence(n+1)/n;
y=map(f,x);
dy=map(df,x);
for(int i=0; i <= n; ++i)
 dot((x[i],y[i]),5bp+blue);
horner h=diffdiv(x,y);
fhorner p=fhorner(h);
draw(graph(p,a,b,n=500),"$x\longmapsto{}L_{"+string(n)+"}$");
draw(graph(f,a,b),red,"$x\longmapsto{}\frac{1}{x^2+1}$");

xlimits(-5,5);
ylimits(-1,1,Crop);

xaxis("$x$",BottomTop,LeftTicks);
yaxis("$y$",LeftRight,RightTicks);

attach(legend(),point(10S),30S);

shipout("runge1");

erase();

// Test 2: The Runge effect in the Hermite interpolation of 1/(x^2+1).

real f(real x) {return(1/(x^2+1));}
real df(real x) {return(-2*x/(x^2+1)^2);}

real a=-5, b=5;
int n=16;
real[] x,y,dy;
x=a+(b-a)*sequence(n+1)/n;
y=map(f,x);
dy=map(df,x);
for(int i=0; i <= n; ++i)
 dot((x[i],y[i]),5bp+blue);
horner h=hdiffdiv(x,y,dy);
fhorner ph=fhorner(h);
draw(graph(p,a,b,n=500),"$x\longmapsto{}H_{"+string(n)+"}$");
draw(graph(f,a,b),red,"$x\longmapsto{}\frac{1}{x^2+1}$");

unitsize(2cm);

xlimits(-5,5);
ylimits(-1,5,Crop);

xaxis("$x$",BottomTop,LeftTicks);
yaxis("$y$",LeftRight,RightTicks);

attach(legend(),point(10S),30S);

shipout("runge2");

erase();

// Test 3: The Runge effect does not occur for all functions:
// Lagrange interpolation of a function whose successive derivatives
// are bounded by a constant M (here M=1) is shown here to converge.

real f(real x) {return(sin(x));}
real df(real x) {return(cos(x));}

real a=-5, b=5;
int n=16;
real[] x,y,dy;
x=a+(b-a)*sequence(n+1)/n;
y=map(f,x);
dy=map(df,x);
for(int i=0; i <= n; ++i)
 dot((x[i],y[i]),5bp+blue);
horner h=diffdiv(x,y);
fhorner p=fhorner(h);

draw(graph(p,a,b,n=500),"$x\longmapsto{}L_{"+string(n)+"}$");
draw(graph(f,a,b),red,"$x\longmapsto{}\cos(x)$");

xaxis("$x$",BottomTop,LeftTicks);
yaxis("$y$",LeftRight,RightTicks);

attach(legend(),point(10S),30S);

shipout("runge3");

erase();

// Test 4: However, one notes here that numerical artifacts may arise
// from limit precision (typically 1e-16).

real f(real x) {return(sin(x));}
real df(real x) {return(cos(x));}

real a=-5, b=5;
int n=72;
real[] x,y,dy;
x=a+(b-a)*sequence(n+1)/n;
y=map(f,x);
dy=map(df,x);
for(int i=0; i <= n; ++i)
 dot((x[i],y[i]),5bp+blue);
horner h=diffdiv(x,y);
fhorner p=fhorner(h);

draw(graph(p,a,b,n=500),"$x\longmapsto{}L_{"+string(n)+"}$");
draw(graph(f,a,b),red,"$x\longmapsto{}\cos(x)$");

ylimits(-1,5,Crop);

xaxis("$x$",BottomTop,LeftTicks);
yaxis("$y$",LeftRight,RightTicks);

attach(legend(),point(10S),30S);

shipout("runge4");

erase();

// Test 5: The situation is much better using Tchebychev points.

unitsize(2cm);

real f(real x) {return(1/(x^2+1));}
real df(real x) {return(-2*x/(x^2+1)^2);}

real a=-5, b=5;
int n=16;
real[] x,y,dy;
fhorner p,ph,ph1;
for(int i=0; i <= n; ++i)
 x[i]=(a+b)/2+(b-a)/2*cos((2*i+1)/(2*n+2)*pi);
y=map(f,x);
dy=map(df,x);
for(int i=0; i <= n; ++i)
 dot((x[i],y[i]),5bp+blue);
horner h=diffdiv(x,y);
fhorner p=fhorner(h);

draw(graph(p,a,b,n=500),"$x\longmapsto{}T_{"+string(n)+"}$");
draw(graph(f,a,b),red,"$x\longmapsto{}\frac{1}{x^2+1}$");

xlimits(-5,5);
ylimits(-1,2,Crop);

xaxis("$x$",BottomTop,LeftTicks);
yaxis("$y$",LeftRight,RightTicks);
attach(legend(),point(10S),30S);

shipout("runge5");

erase();

// Test 6: Adding a few more Tchebychev points yields a very good result.

unitsize(2cm);

real f(real x) {return(1/(x^2+1));}
real df(real x) {return(-2*x/(x^2+1)^2);}

real a=-5, b=5;
int n=26;
real[] x,y,dy;
for(int i=0; i <= n; ++i)
 x[i]=(a+b)/2+(b-a)/2*cos((2*i+1)/(2*n+2)*pi);
y=map(f,x);
dy=map(df,x);
for(int i=0; i <= n; ++i)
 dot((x[i],y[i]),5bp+blue);
horner h=diffdiv(x,y);
fhorner p=fhorner(h);
draw(graph(p,a,b,n=500),"$x\longmapsto{}T_{"+string(n)+"}$");
draw(graph(f,a,b),red,"$x\longmapsto{}\frac{1}{x^2+1}$");

xlimits(-5,5);
ylimits(-1,2,Crop);

xaxis("$x$",BottomTop,LeftTicks);
yaxis("$y$",LeftRight,RightTicks);
attach(legend(),point(10S),30S);


shipout("runge6");

erase();

// Test 7: Another Tchebychev example.

unitsize(2cm);

real f(real x) {return(sqrt(abs(x-1)));}

real a=-2, b=2;
int n=30;
real[] x,y,dy;
for(int i=0; i <= n; ++i)
 x[i]=(a+b)/2+(b-a)/2*cos((2*i+1)/(2*n+2)*pi);
y=map(f,x);
dy=map(df,x);
for(int i=0; i <= n; ++i)
 dot((x[i],y[i]),5bp+blue);
horner h=diffdiv(x,y);
fhorner p=fhorner(h);
draw(graph(p,a,b,n=500),"$x\longmapsto{}T_{"+string(n)+"}$");
draw(graph(f,a,b),red,"$x\longmapsto{}\sqrt{|x-1|}$");

xlimits(-2,2);
ylimits(-0.5,2,Crop);

xaxis("$x$",BottomTop,LeftTicks);
yaxis("$y$",LeftRight,RightTicks);
attach(legend(),point(10S),30S);

shipout("runge7");

You are viewing proxied material from ftp.icm.edu.pl. The copyright of proxied material belongs to its original authors. Any comments or complaints in relation to proxied material should be directed to the original authors of the content concerned. Please see the disclaimer for more details.