1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
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);
};
}
|
