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
|
(* -*-Mathematica-*- script to determine the polynomial coefficients
to approximate functions
2013-02-12 Erik Schnetter <eschnetter@perimeterinstitute.ca>
Based on the "minimax" algorithm, but using least squares instead *)
funcs =
{
{Sin[2 Pi #]&, 0, 1/4, 1, 2, 10, 1, 1, 1, {0, 1/4}},
{Cos[2 Pi #]&, 0, 1/4, 0, 2, 10, 1, 1, 1, {0, 1/4}},
(* {Tan[Pi #]&, 0, 1/2, 1, 2, 3, 2, 2, 3, {0, 1/4}}, *)
{Log[2, (1+#)/(1-#)]&, -3 + 2 Sqrt[2], 3 - 2 Sqrt[2],
0, 1, 15, 0, 1, 0, {0}},
{2^#&, -1/2, +1/2, 0, 1, 15, 0, 1, 1, {}},
{ArcTan, 0, 1/2, 1, 2, 15, 1, 2, 1, {0, 1}}
};
findcoeffs[func_, xmin_, xmax_, dmin_, dstep_, ndegrees_, denndegrees_,
conslocs_] :=
Module[{prec, npts,
degree,
qs, q2x, x2q,
funcq,
funcqpade, polydenom,
coeffs, powers, poly,
norm,
A, r, b,
approx,
error1, error2, error3, error},
Print[{"function ", func, ndegrees, denndegrees}];
(* Working precision *)
prec = 30;
(* Number of test points *)
npts = 10000;
degree = dmin + ndegrees dstep;
(* Test points in normalized interval [0,1] *)
qs = Table[n/(npts-1), {n, 0, npts-1}];
(* A (discrete) L2-norm based on the test points *)
norm[f_] = Sqrt[Total[N[(Abs[f[#1]]^2 &) /@ qs, prec]] / Length[qs]];
(* Transform q to x coordinate *)
q2x[q_] = xmin + (xmax - xmin) q;
x2q[x_] = (x - xmin) / (xmax - xmin);
(* Function with rescaled input *)
funcq[q_] = func[q2x[q]];
(* Use denominator of Pade approximant as denominator of our
approximant *)
funcqpade[q_] = PadeApproximant[funcq[q], {q, 0, denndegrees}];
polydenom[q_] = Denominator[funcqpade[q]];
(* List of expansion coefficients *)
coeffs = Table[c[i], {i, dmin, degree, dstep}];
(* Corresponding list of powers of q *)
powers[q_] = Table[If[i==0, 1, q^i], {i, dmin, degree, dstep}];
(* Polynomial *)
poly[q_] = coeffs . powers[q];
(* We determine the expansion coefficients via a
least-squares method *)
A = N[Table[powers[q], {q, qs}], prec];
r = N[(funcq[#1] polydenom[#1] &) /@ qs, prec];
(* r = N[(Limit[funcq[q] polydenom[q], q->#1] &) /@ qs, prec]; *)
b = LeastSquares[A, r];
(* Define approximating polynomial using this solution *)
approx[q_] = ((poly[q] /. MapThread[#1 -> #2 &, {coeffs, b}]) /
N[polydenom[q], prec]);
(* Calculate three kinds of errors to check solution: *)
(* (1) the (discrete) norm form above: *)
error1 = norm[approx[#1] - funcq[#1] &];
(* (2) A non-discrete L2-norm using an integral: *)
error2 = Sqrt[NIntegrate[Abs[approx[q] - funcq[q]]^2, {q, 0, 1},
WorkingPrecision -> prec]];
(* (3) The maximum of the error: *)
error3 = NMaxValue[{Abs[approx[q] - funcq[q]], q >= 0 && q <= 1}, {q},
WorkingPrecision -> prec];
error = Max[error1, error2, error3];
(* Evaluate at constraint locations to check solution: *)
consvals = Map[Abs[approx[x2q[#]] - func[#]]&, conslocs];
Print[{"error ", CForm[error], " constraints ", CForm[consvals]}];
Write[outfile, {func, ndegrees, denndegrees, CForm[error],
CForm[HornerForm[approx[x2q[x]]]]}]];
findcoeffs2[func_, xmin_, xmax_,
dmin_, dstep_, maxndegrees_,
dendmin_, dendstep_, maxdenndegrees_,
conslocs_] :=
Do[findcoeffs[func, xmin, xmax, dmin, dstep, ndegrees, denndegrees, conslocs],
{ndegrees, maxndegrees},
{denndegrees, dendmin, dendstep maxdenndegrees, dendstep}];
outfile = OpenWrite["/tmp/tmp"];
(* findcoeffs[Tan[Pi #]&, 0, 1/2, 1, 2, 3, 4, {0, 1/4}]; *)
(* findcoeffs[Tan[Pi #]&, 0, 1/2, 1, 2, 3, 4, {0, 1/4}]; *)
(* findcoeffs2[Tan[Pi #]&, 0, 1/2, 1, 2, 5, 2, 2, 5, {0, 1/4}]; *)
(*
PADE DOES NOT REMOVE SINGULARITIES
PROBABLY NEED TO DO THAT MANUALLY.
*)
(* findcoeffs[ArcTan[#]&, 0, 1, 1, 2, 10, 1, {0, 1}]; *)
(* findcoeffs2[ArcTan, 0, 1, 1, 2, 10, 1, 2, 10, {0, 1}]; *)
(* findcoeffs2[ArcTan[#]&, 0, 1/2, 1, 2, 15, 1, 2, 1, {0, 1/2}]; *)
(* findcoeffs2[Log[2, (1+#)/(1-#)]&,
-3 + 2 Sqrt[2], 3 - 2 Sqrt[2], 0, 1, 15, 0, 1, 0, {0}]; *)
outfile = OpenWrite["coeffs.out"];
Write[outfile, "(* Coefficients for function approximations *)"];
Map[findcoeffs2@@# &, funcs];
Close[outfile];
|
