diff options
author | Karl Berry <karl@freefriends.org> | 2013-04-07 18:19:31 +0000 |
---|---|---|
committer | Karl Berry <karl@freefriends.org> | 2013-04-07 18:19:31 +0000 |
commit | 752012c605d34cd943795527a9738475a6958fcc (patch) | |
tree | 4ee06acdd8333a662c2d6f6ef716235053468f55 /Master/texmf/asymptote/interpolate.asy | |
parent | 9789d09132f18a838e84f041b4b3aff28d3426ec (diff) |
texmf -> texmf-dist: start with unique dirs from texmf
git-svn-id: svn://tug.org/texlive/trunk@29712 c570f23f-e606-0410-a88d-b1316a301751
Diffstat (limited to 'Master/texmf/asymptote/interpolate.asy')
-rw-r--r-- | Master/texmf/asymptote/interpolate.asy | 140 |
1 files changed, 0 insertions, 140 deletions
diff --git a/Master/texmf/asymptote/interpolate.asy b/Master/texmf/asymptote/interpolate.asy deleted file mode 100644 index a600c94c125..00000000000 --- a/Master/texmf/asymptote/interpolate.asy +++ /dev/null @@ -1,140 +0,0 @@ -// 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); - }; -} |