1% ---------------------------------------------------------------------------
2% Driver (owned by the app, not editable in the browser).
3%
4% It reads params.json, calls into the method script, and writes out.json.
5% The editable method script is appended to the end of this file before the
6% run, so the functions it defines -- bary_weights, bary_eval, and optionally
7% local_blend -- are local functions of this script and are visible here.
8% Every driver-owned helper is prefixed drv_ so that it cannot collide with
9% anything the method script defines.
10% ---------------------------------------------------------------------------
12drv_p = jsondecode(fileread('params.json'));
13switch drv_p.mode
14 case 'explore'
15 drv_out = drv_explore(drv_p);
16 case 'converge'
17 drv_out = drv_converge(drv_p);
18 otherwise
19 error('driver: unknown mode ''%s''', drv_p.mode);
20end
21drv_fid = fopen('out.json', 'w');
22fprintf(drv_fid, '%s', jsonencode(drv_out));
23fclose(drv_fid);
25% ── one set of nodes, one degree d: everything the first three tabs draw ────
26function out = drv_explore(p)
27x = nodes_of(p.nodes, p.a, p.b, p.n, p.seed);
28n = numel(x) - 1;
29d = min(max(round(p.d), 0), n);
30y = testfun(p.f, x, p.fexpr);
31t = linspace(p.a, p.b, p.ngrid);
32ft = testfun(p.f, t, p.fexpr);
34w = bary_weights(x, d);
35w = reshape(w, 1, []);
36r = bary_eval(x, y, w, t);
37r = reshape(r, 1, []);
39out = struct();
40out.n = n;
41out.d = d;
42out.x = x;
43out.y = y;
44out.t = t;
45out.ft = ft;
46out.r = r;
47out.err = r - ft;
48out.maxerr = drv_maxabs(r - ft);
50% Barycentric weights, and the integer form they take on a uniform mesh.
51% Section 4 lists these: 1,1,...,1,1 for d = 0, then 1,2,2,...,2,2,1 for
52% d = 1, 1,3,4,...,4,3,1 for d = 2, and so on.
53out.w = w;
54aw = abs(w);
55pos = aw(aw > 0);
56if isempty(pos)
57 base = 1;
58else
59 base = min(pos);
60end
61out.wscaled = aw / base;
62out.wsign = sign(w);
63out.wIsInteger = all(abs(out.wscaled - round(out.wscaled)) < 1e-7);
64out.wAlternates = all(w(1:end - 1) .* w(2:end) < 0);
66% the degree-n polynomial interpolant, in barycentric form (equation 2)
67if p.want.poly
68 out.rpoly = reshape(bary_eval(x, y, drv_lagrange_weights(x), t), 1, []);
69end
71% the clamped C^2 cubic spline of Tables 3 and 4
72if p.want.spline
73 [~, dya] = testfun(p.f, p.a, p.fexpr);
74 [~, dyb] = testfun(p.f, p.b, p.fexpr);
75 out.rspline = reshape(cubic_spline(x, y, dya, dyb, t), 1, []);
76end
78% the blend of local polynomials, equations (4) and (5)
79if p.want.blend
80 out.hasBlend = false;
81 try
82 [P, L] = local_blend(x, y, d, t);
83 out.P = drv_rows(P);
84 out.L = drv_rows(L);
85 out.wlo = x(1:(n - d + 1));
86 out.whi = x((d + 1):(n + 1));
87 out.hasBlend = true;
88 catch blenderr
89 out.blendError = blenderr.message;
90 end
91end
93if p.want.poles
94 out.poles = drv_poles(x, y, w, p);
95end
96end
98% ── denominator of r, its sign on the real line, and its roots ─────────────
99function s = drv_poles(x, y, w, p)
100n = numel(x) - 1;
101mg = 0.35 * (p.b - p.a);
102tw = linspace(p.a - mg, p.b + mg, p.ngridwide);
104% Work in a variable rescaled to [-1, 1]. This only multiplies the
105% denominator by a positive constant, so signs and roots are untouched, but it
106% keeps the products of n factors from over- or underflowing.
107c0 = (p.a + p.b) / 2;
108sc = (p.b - p.a) / 2;
109xs = (x - c0) / sc;
110ts = (tw - c0) / sc;
112q = drv_denom_eval(xs, w, ts);
114s = struct();
115s.t = tw;
116% A plain plot of q is useless: it spans many orders of magnitude. The signed
117% n-th root keeps every sign and every zero exactly where it was and brings the
118% magnitudes into a range that can be drawn.
119qmax = max(abs(q));
120if qmax == 0
121 qmax = 1;
122end
123s.u = sign(q) .* (abs(q) / qmax).^(1 / max(1, n));
125% Real poles: sign changes of the denominator, located by linear interpolation
126% of the crossing. This works for any n, which the root-finding below does not.
127sg = sign(q);
128idx = find(sg(1:end - 1) .* sg(2:end) < 0);
129rp = zeros(1, numel(idx));
130for m = 1:numel(idx)
131 j = idx(m);
132 rp(m) = tw(j) + (tw(j + 1) - tw(j)) * abs(q(j)) / (abs(q(j)) + abs(q(j + 1)));
133end
134s.realPoles = rp;
136% All roots in the complex plane. Theorem 1 says none of them are real, and
137% that is the picture: every root sits off the real axis. We only do this for
138% moderate n, since finding roots from monomial coefficients is not reliable
139% for large degree.
140s.rootsRe = [];
141s.rootsIm = [];
142s.rootsShown = false;
143if n <= p.rootsMaxN
144 [cc, floorc] = drv_denom_coeffs(xs, w);
145 % The top d coefficients of q vanish identically -- that is exactly the
146 % statement that r reproduces polynomials of degree d, i.e. that the first
147 % d moments of the weights are zero -- but they come out of the sum as
148 % roundoff rather than as zero. Left in place they contribute spurious
149 % roots, some of them real, which would be a lie in a picture whose whole
150 % point is that there are no real roots. So drop every leading
151 % coefficient that is below the rounding-error floor of its own sum.
152 k0 = find(abs(cc) > floorc, 1);
153 if ~isempty(k0) && numel(cc) - k0 >= 1
154 z = reshape(roots(cc(k0:end)), 1, []);
155 s.rootsRe = real(z) * sc + c0;
156 s.rootsIm = imag(z) * sc;
157 s.rootsShown = true;
158 end
159end
161% the classical alternative, with the poles it puts in the interval
162if p.want.classical
163 [rc, cp] = classical_rational(x, y, tw);
164 s.classical = reshape(rc, 1, []);
165 s.classicalPoles = cp;
166end
167end
169% q(t) = sum_k w_k prod_{j ~= k} (t - x_j), the denominator of r written as a
170% polynomial. Its zeros are exactly the poles of r.
171function q = drv_denom_eval(xs, w, ts)
172D = ts(:) - xs(:).';
173n1 = numel(xs);
174q = zeros(numel(ts), 1);
175for k = 1:n1
176 pr = ones(numel(ts), 1);
177 for j = 1:n1
178 if j ~= k
179 pr = pr .* D(:, j);
180 end
181 end
182 q = q + w(k) * pr;
183end
184q = q.';
185end
187function [c, floorc] = drv_denom_coeffs(xs, w)
188n1 = numel(xs);
189n = n1 - 1;
190c = zeros(1, n1);
191for k = 1:n1
192 ck = 1;
193 for j = 1:n1
194 if j ~= k
195 ck = conv(ck, [1, -xs(j)]);
196 end
197 end
198 c = c + w(k) * ck;
199end
200% Rounding-error floor, coefficient by coefficient. The nodes have been
201% rescaled so that every |x_j| <= 1, so the coefficient of t^(n-m) in each of
202% the n+1 products is at most binomial(n, m); summing them accumulates at most
203% n+1 roundings of terms of size |w_k| times that bound.
204floorc = zeros(1, n1);
205sw = sum(abs(w));
206for m = 0:n
207 floorc(m + 1) = 8 * eps * n1 * sw * nchoosek(n, m);
208end
209end
211% the weights of equation (2), which put the degree-n polynomial interpolant
212% into barycentric form
213function w = drv_lagrange_weights(x)
214n1 = numel(x);
215w = zeros(1, n1);
216for k = 1:n1
217 pr = 1;
218 for j = 1:n1
219 if j ~= k
220 pr = pr / (x(k) - x(j));
221 end
222 end
223 w(k) = pr;
224end
225end
227% ── error against n, for a range of d: Tables 1 to 4 ───────────────────────
228function out = drv_converge(p)
229ns = reshape(p.ns, 1, []);
230ds = reshape(p.ds, 1, []);
231t = linspace(p.a, p.b, p.ngrid);
232ft = testfun(p.f, t, p.fexpr);
233[~, dya] = testfun(p.f, p.a, p.fexpr);
234[~, dyb] = testfun(p.f, p.b, p.fexpr);
236E = zeros(numel(ds), numel(ns));
237esp = zeros(1, numel(ns));
238epo = zeros(1, numel(ns));
240for ib = 1:numel(ns)
241 n = ns(ib);
242 x = nodes_of(p.nodes, p.a, p.b, n, p.seed);
243 y = testfun(p.f, x, p.fexpr);
244 for ia = 1:numel(ds)
245 d = min(ds(ia), n);
246 w = reshape(bary_weights(x, d), 1, []);
247 E(ia, ib) = drv_maxabs(reshape(bary_eval(x, y, w, t), 1, []) - ft);
248 end
249 if p.want.spline
250 esp(ib) = drv_maxabs(reshape(cubic_spline(x, y, dya, dyb, t), 1, []) - ft);
251 end
252 if p.want.poly
253 wl = drv_lagrange_weights(x);
254 epo(ib) = drv_maxabs(reshape(bary_eval(x, y, wl, t), 1, []) - ft);
255 end
256end
258out = struct();
259out.ns = ns;
260out.ds = ds;
261out.E = drv_rows(E);
262out.orders = drv_rows(drv_orders(ns, E));
263if p.want.spline
264 out.splineErr = esp;
265 out.splineOrders = drv_orders(ns, esp);
266end
267if p.want.poly
268 out.polyErr = epo;
269end
270end
272function o = drv_orders(ns, E)
273o = zeros(size(E));
274o(:, 1) = NaN;
275for j = 2:numel(ns)
276 o(:, j) = log(E(:, j - 1) ./ E(:, j)) / log(ns(j) / ns(j - 1));
277end
278end
280% jsonencode flattens a matrix with a single row, so hand matrices over as a
281% cell array of rows to keep the shape on the JavaScript side predictable.
282function c = drv_rows(A)
283c = cell(1, size(A, 1));
284for i = 1:size(A, 1)
285 c{i} = A(i, :);
286end
287end
289function m = drv_maxabs(v)
290v = abs(v(:));
291v = v(isfinite(v));
292if isempty(v)
293 m = Inf;
294else
295 m = max(v);
296end
297end