Re: Fun math puzzle: cin(X)

[Dirk Laurie <dirk.laurie@...>]

Op Di. 9 Apr. 2019 om 00:09 het Albert Chan <albertmcchan@yahoo.com> geskryf:
>
>
> > On Apr 4, 2019, at 5:25 PM, Egor Skriptunoff <egor.skriptunoff@gmail.com> wrote:
> >
> > The most straightforward approach is
> > to build Taylor series of cin(x) one term at a time.
> >
> >    function maclaurin_of_cin(k)
> >       for n = #c + 1, k do
> >          a = {}
> >          local e, h = f(c), f(d)
> >          s, a = -s/(2*n)/(2*n+1)
> >          local t = (s-h-e)/3
> >          assert(math.abs(t) < 0.056)
> >          c[n], d[n] = t, e + 2*t
> >       end
> >       return c[k]
> >    end
>
> Say, I wanted a function, din(x), such that din(din(din(din(x)))) = sin(x)
>
> What is required to modify above to do maclaurin_of_din(k) ?

I can't speak for the code, but the usual method is equivalent to
Newton's method for roots, applied to functions. It will work for any
function s whose Maclaurin series starts at x.

The well-known algorithm for square roots of real numbers is
   x ← x + (a/x-x)/2.
For n-th roots, it becomes
   x ← x + (a/x^(n-1)-x)/n
For n-th root of a function, it becomes:
   f ← f + (g(g(...(s)))-f)/n    -- n-1 applications of g
where g is the inverse function of f.

You can do this analytically, but usually for one step only; the
inverse function soon starts to become intractable.

Examples for din:
   f(x)=x, g(x)=x, next approximation is x + (sin(x)-x)/3
   f(x) = 2*sin(x/2), g(x)=2*asin(x/2), next approximation is
2*sin(x/2) + (2*asin(2*asin(2*asin(sin(x)/2)/2)/2)-2*sin(x/2))/4

This second approximation is not bad at all, applied four times
starting at 1 gives 0.84465, should be 0.84147. You can draw a
reasonable graph with that. Applied four times at pi/2, gives 1.058,
still not outrageously bad.

What we do instead is to work with a power series for f; then the
reversal of the series is a truncated power series for g. We add just
one term to both series at a time. You already know up to x^k the
identity holds, so you compare coefficients of x^(k+1). I think that
is the basis for Egor's code, but my letters f and g probably mean
different things to his f and g.