// Lagrange and Hermite interpolation in Asymptote // Author: Olivier Guibé // Acknowledgements: Philippe Ivaldi // diffdiv(x,y) computes Newton's Divided Difference for // Lagrange interpolation with distinct values {x_0,..,x_n} in the array x // and values y_0,...,y_n in the array y, // hdiffdiv(x,y,dyp) computes Newton's Divided Difference for // Hermite interpolation where dyp={dy_0,...,dy_n}. // // fhorner(x,coeff) uses Horner's rule to compute the polynomial // a_0+a_1(x-x_0)+a_2(x-x_0)(x-x_1)+...+a_n(x-x_0)..(x-x_{n-1}), // where coeff={a_0,a_1,...,a_n}. // fspline does standard cubic spline interpolation of a function f // on the interval [a,b]. // The points a=x_1 < x_2 < .. < x_n=b form the array x; // the points y_1=f(x_1),....,y_n=f(x_n) form the array y // We use the Hermite form for the spline. // The syntax is: // s=fspline(x,y); default not_a_knot condition // s=fspline(x,y,natural); natural spline // s=fspline(x,y,periodic); periodic spline // s=fspline(x,y,clamped(1,1)); clamped spline // s=fspline(x,y,monotonic); piecewise monotonic spline // Here s is a real function that is constant on (-infinity,a] and [b,infinity). private import math; import graph_splinetype; typedef real fhorner(real); struct horner { // x={x0,..,xn}(not necessarily distinct) // a={a0,..,an} corresponds to the polyonmial // a_0+a_1(x-x_0)+a_2(x-x_0)(x-x_1)+...+a_n(x-x_0)..(x-x_{n-1}), real[] x; real[] a; } // Evaluate p(x)=d0+(x-x0)(d1+(x-x1)+...+(d(n-1)+(x-x(n-1))*dn))) // via Horner's rule: n-1 multiplications, 2n-2 additions. fhorner fhorner(horner sh) { int n=sh.x.length; checklengths(n,sh.a.length); return new real(real x) { real s=sh.a[n-1]; for(int k=n-2; k >= 0; --k) s=sh.a[k]+(x-sh.x[k])*s; return s; }; } // Newton's Divided Difference method: n(n-1)/2 divisions, n(n-1) additions. horner diffdiv(real[] x, real[] y) { int n=x.length; horner s; checklengths(n,y.length); for(int i=0; i < n; ++i) s.a[i]=y[i]; for(int k=0; k < n-1; ++k) { for(int i=n-1; i > k; --i) { s.a[i]=(s.a[i]-s.a[i-1])/(x[i]-x[i-k-1]); } } s.x=x; return s; } // Newton's Divided Difference for simple Hermite interpolation, // where one specifies both p(x_i) and p'(x_i). horner hdiffdiv(real[] x, real[] y, real[] dy) { int n=x.length; horner s; checklengths(n,y.length); checklengths(n,dy.length); for(int i=0; i < n; ++i) { s.a[2*i]=y[i]; s.a[2*i+1]=dy[i]; s.x[2*i]=x[i]; s.x[2*i+1]=x[i]; } for(int i=n-1; i > 0; --i) s.a[2*i]=(s.a[2*i]-s.a[2*i-2])/(x[i]-x[i-1]); int stop=2*n-1; for(int k=1; k < stop; ++k) { for(int i=stop; i > k; --i) { s.a[i]=(s.a[i]-s.a[i-1])/(s.x[i]-s.x[i-k-1]); } } return s; } typedef real realfunction(real); // piecewise Hermite interpolation: // return the piecewise polynomial p(x), where on [x_i,x_i+1], deg(p) <= 3, // p(x_i)=y_i, p(x_{i+1})=y_i+1, p'(x_i)=dy_i, and p'(x_{i+1})=dy_i+1. // Outside [x_1,x_n] the returned function is constant: y_1 on (infinity,x_1] // and y_n on [x_n,infinity). realfunction pwhermite(real[] x, real[] y, real[] dy) { int n=x.length; checklengths(n,y.length); checklengths(n,dy.length); if(n < 2) abort(morepoints); if(!increasing(x,strict=true)) abort("array x is not strictly increasing"); return new real(real t) { int i=search(x,t); if(i == n-1) { i=n-2; t=x[n-1]; } else if(i == -1) { i=0; t=x[0]; } real h=x[i+1]-x[i]; real delta=(y[i+1]-y[i])/h; real e=(3*delta-2*dy[i]-dy[i+1])/h; real f=(dy[i]-2*delta+dy[i+1])/h^2; real s=t-x[i]; return y[i]+s*(dy[i]+s*(e+s*f)); }; } realfunction fspline(real[] x, real[] y, splinetype splinetype=notaknot) { real[] dy=splinetype(x,y); return new real(real t) { return pwhermite(x,y,dy)(t); }; }