Derivative of an interpol function, what to do with charfun2

On 2013-11-23, ???????? ?????? <vonavy15 at gmail.com> wrote:

> f:log(%e^(-x)+sin(y))$
> y0:1$
> x1:0$
> x2:10$
> n:40$
> s:float((x2-x1)/n);
> results: rk(f,y,y0,[x,x1,x2,s])$
> load(interpol)$
> g:cspline(results)$
> d:diff(g,x);
> e:d-f;
>
> I want to get some values of e between x1 and x2 and output a descrete plot
> of this values.

Hmm, that's an interesting problem. I think I've solved it. The main
problem is that we have to tell Maxima what's the derivative of a
characteristic function (it doesn't know it by default). There are
other minor speedbumps along the way.

    /* (0) */
    fpprintprec : 4 $
    display2d : false $
    ratprint : false $
    
    f:log(%e^(-x)+sin(y))$
    y0:1$
    x1:0$
    x2:10$
    /* (1) */ n : 8 $
    s:float((x2-x1)/n);
    results: rk(f,y,y0,[x,x1,x2,s])$
    load(interpol)$
    g:cspline(results)$
    /* (2) */ g : ''g;
    d:diff(g,x);
    
    /* (3) */
    matchdeclare ([aa, bb], all, xx, symbolp);
    tellsimp ('diff (charfun (aa <= xx), xx), delta (xx - aa));
    tellsimp ('diff (charfun (aa <= xx and xx < bb), xx), delta (xx - aa) - delta (xx - bb));
    tellsimp ('diff (charfun (xx < bb), xx), - delta (xx - bb));
    
    /* (4) */
    matchdeclare (nn, notequal (0));
    tellsimp (delta (nn), 0);
    
    /* (5) */ d : ''d;
    
    /* (6) */
    FOO (p) := ev (d, x=p[1]);
    BAR (p) := ev (f, x=p[1], y=p[2]);
    e : map (lambda ([p], FOO (p) - BAR (p)), results);
    
    /* (7) */
    subst (delta (0.0) = 0, e);

(0) incidental
(1) make problem a little smaller for convenience
(2) cspline returns unevaluated charfun2 expressions; evaluate them
(3) tell Maxima what's the derivative of charfun
(4) tell Maxima what's the value of delta away from 0
(5) apply rules for charfun and delta
(6) spline derivative (FOO) minus exact derivative (BAR)
(7) At this point I see that e contains terms like <eps>*delta(0.0)
where <eps> is on the order of 1e-16. These must be due only to
round-off error (there can't really be any delta terms) so I'll clobber
them.

For the final result I get

    [-.1316,.03553,-.009746,.002185,-7.5764E-4,1.5791E-4,-7.8254E-5,
     9.9061E-5,-3.326E-4]

I guess the spline derivative is a poor approximation near x = 0,
but it is more accurate farther away.

Hope this helps, and thanks for the interesting problem.

Robert Dodier