// 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");