summaryrefslogtreecommitdiff
path: root/graphics/asymptote/examples/interpolate1.asy
diff options
context:
space:
mode:
authorNorbert Preining <norbert@preining.info>2019-09-02 13:46:59 +0900
committerNorbert Preining <norbert@preining.info>2019-09-02 13:46:59 +0900
commite0c6872cf40896c7be36b11dcc744620f10adf1d (patch)
tree60335e10d2f4354b0674ec22d7b53f0f8abee672 /graphics/asymptote/examples/interpolate1.asy
Initial commit
Diffstat (limited to 'graphics/asymptote/examples/interpolate1.asy')
-rw-r--r--graphics/asymptote/examples/interpolate1.asy226
1 files changed, 226 insertions, 0 deletions
diff --git a/graphics/asymptote/examples/interpolate1.asy b/graphics/asymptote/examples/interpolate1.asy
new file mode 100644
index 0000000000..8845c8dd2f
--- /dev/null
+++ b/graphics/asymptote/examples/interpolate1.asy
@@ -0,0 +1,226 @@
+// 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");