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
|
#include <2geom/chebyshev.h>
#include <2geom/sbasis.h>
#include <2geom/sbasis-poly.h>
#include <vector>
using std::vector;
#include <gsl/gsl_math.h>
#include <gsl/gsl_chebyshev.h>
namespace Geom{
SBasis cheb(unsigned n) {
static std::vector<SBasis> basis;
if(basis.empty()) {
basis.push_back(Linear(1,1));
basis.push_back(Linear(0,1));
}
for(unsigned i = basis.size(); i <= n; i++) {
basis.push_back(Linear(0,2)*basis[i-1] - basis[i-2]);
}
return basis[n];
}
SBasis cheb_series(unsigned n, double* cheb_coeff) {
SBasis r;
for(unsigned i = 0; i < n; i++) {
double cof = cheb_coeff[i];
//if(i == 0)
//cof /= 2;
r += cheb(i)*cof;
}
return r;
}
SBasis clenshaw_series(unsigned m, double* cheb_coeff) {
/** b_n = a_n
b_n-1 = 2*x*b_n + a_n-1
b_n-k = 2*x*b_{n-k+1} + a_{n-k} - b_{n - k + 2}
b_0 = x*b_1 + a_0 - b_2
*/
double a = -1, b = 1;
SBasis d, dd;
SBasis y = (Linear(0, 2) - (a+b)) / (b-a);
SBasis y2 = 2*y;
for(int j = m - 1; j >= 1; j--) {
SBasis sv = d;
d = y2*d - dd + cheb_coeff[j];
dd = sv;
}
return y*d - dd + 0.5*cheb_coeff[0];
}
SBasis chebyshev_approximant (double (*f)(double,void*), int order, Interval in, void* p) {
gsl_cheb_series *cs = gsl_cheb_alloc (order+2);
gsl_function F;
F.function = f;
F.params = p;
gsl_cheb_init (cs, &F, in[0], in[1]);
SBasis r = compose(clenshaw_series(order, cs->c), Linear(-1,1));
gsl_cheb_free (cs);
return r;
}
struct wrap {
double (*f)(double,void*);
void* pp;
double fa, fb;
Interval in;
};
double f_interp(double x, void* p) {
struct wrap *wr = (struct wrap *)p;
double z = (x - wr->in[0]) / (wr->in[1] - wr->in[0]);
return (wr->f)(x, wr->pp) - ((1 - z)*wr->fa + z*wr->fb);
}
SBasis chebyshev_approximant_interpolating (double (*f)(double,void*),
int order, Interval in, void* p) {
double fa = f(in[0], p);
double fb = f(in[1], p);
struct wrap wr;
wr.fa = fa;
wr.fb = fb;
wr.in = in;
//printf("%f %f\n", fa, fb);
wr.f = f;
wr.pp = p;
return compose(Linear(in[0], in[1]), Linear(fa, fb)) + chebyshev_approximant(f_interp, order, in, &wr) + Linear(fa, fb);
}
SBasis chebyshev(unsigned n) {
static std::vector<SBasis> basis;
if(basis.empty()) {
basis.push_back(Linear(1,1));
basis.push_back(Linear(0,1));
}
for(unsigned i = basis.size(); i <= n; i++) {
basis.push_back(Linear(0,2)*basis[i-1] - basis[i-2]);
}
return basis[n];
}
};
/*
Local Variables:
mode:c++
c-file-style:"stroustrup"
c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
indent-tabs-mode:nil
fill-column:99
End:
*/
// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :
|