summaryrefslogtreecommitdiff
path: root/Master/texmf-dist/asymptote/interpolate.asy
diff options
context:
space:
mode:
Diffstat (limited to 'Master/texmf-dist/asymptote/interpolate.asy')
-rw-r--r--Master/texmf-dist/asymptote/interpolate.asy140
1 files changed, 140 insertions, 0 deletions
diff --git a/Master/texmf-dist/asymptote/interpolate.asy b/Master/texmf-dist/asymptote/interpolate.asy
new file mode 100644
index 00000000000..a600c94c125
--- /dev/null
+++ b/Master/texmf-dist/asymptote/interpolate.asy
@@ -0,0 +1,140 @@
+// 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);
+ };
+}