fix Jacobian culculation of ikeda map

[hibiki-kato]

Summary

• Ikeda map's Jacobian corrected.
• Tests: randomized states and parameter grids; compare against ForwardDiff.
• Symbolic verification script included.
• No API changes; only correctness fixes.

About

By comparing with ForwardDiff, I noticed that the analytic Jacobian for the Ikeda map is not correct, and its source https://www.math.arizona.edu/~ura-reports/001/huang.pojen/2000_Report.html has expired.

This PR replaces the Jacobian with the correct closed-form expression, adds a robust test comparing it to the ForwardDiff Jacobian over randomized state and parameter grids, and includes a symbolic validation script (Mathematica/Wolfram).

Correct Jacobian (derivation sketch)

Ikeda map:

$$ \begin{aligned} t &= c - \frac{d}{1 + x_n^2 + y_n^2} \\ x_{n+1} &= a + b(x_n \cos(t) - y\sin(t)) \\ y_{n+1} &= b(x\sin(t) + y \cos(t)) \end{aligned} $$

For

$$ t=c-\frac{d}{1+x^2+y^2},\quad x’ = a+b(x\cos t-y\sin t),\quad y’ = b(x\sin t+y\cos t), $$

we have

$$ \partial_x t = \frac{2dx}{(1+x^2+y^2)^2},\qquad \partial_y t = \frac{2dy}{(1+x^2+y^2)^2}. $$

Hence

$$ \begin{aligned} J_{11}&=b\left(\cos t-\frac{2d}{(1+r^2)^2}[x^2\sin t+xy\cos t]\right),\\ J_{12}&=b\left(-\sin t-\frac{2d}{(1+r^2)^2}[xy\sin t+y^2\cos t]\right),\\ J_{21}&=b\left(\ \sin t+\frac{2d}{(1+r^2)^2}[x^2\cos t-xy\sin t]\right),\\ J_{22}&=b\left(\cos t+\frac{2d}{(1+r^2)^2}[xy\cos t-y^2\sin t]\right), \end{aligned} $$

where $\quad r^2=x^2+y^2$.

Julia Test

For various parameters and random points, the analytical Jacobian and the ForwardDiff Jacobian are verified to be approximately the same.

 using Random, Test
using StaticArrays
using DynamicalSystems
using PredefinedDynamicalSystems

# --- Correct analytic Jacobian (same as above) ---
@inbounds function ikedamap_jacob(u, p, n)
    a, b, c, d = p
    x, y = u
    t  = c - d/(1 + x^2 + y^2)
    s, cc = sin(t), cos(t)
    r2 = 1 + x^2 + y^2
    aux = 2d/(r2^2)

    J11 = b*( cc - aux*( x^2*s + x*y*cc) )
    J12 = b*( -s - aux*( x*y*s + y^2*cc) )
    J21 = b*(  s + aux*( x^2*cc - x*y*s) )
    J22 = b*( cc + aux*( x*y*cc - y^2*s) )

    return @SMatrix [J11 J12; J21 J22]
end

# ---------------- Tests vs ForwardDiff ----------------
Random.seed!(1234)

u0s = [rand(2) for _ in 1:100]

as = 0.9:0.01:1.1
bs = 0.85:0.01:0.95
cs = 0.35:0.01:0.45
ds = 5.5:0.1:6.5

@testset "Ikeda Jacobian ≈ ForwardDiff" begin
    for a in as, b in bs, c in cs, d in ds
        # PredefinedDynamicalSystems.ikedamap expects named params;
        # our J expects a positional container (tuple or vector). We’ll pass both.
        p_named = (a=a, b=b, c=c, d=d)
        p_vec   = [a, b, c, d]

        for u0 in u0s
            ds_ = PredefinedDynamicalSystems.ikedamap(u0; p_named...)
            tands    = TangentDynamicalSystem(ds_; J=ikedamap_jacob)  # analytic J
            tands_FD = TangentDynamicalSystem(ds_)                    # ForwardDiff J

            J_analytic = tands.J(u0, p_vec, 0)
            J_fd       = tands_FD.J(u0, p_vec, 0)

            @test isapprox(J_analytic, J_fd; atol=1e-10, rtol=1e-10)
        end
    end
end

Output

Test Summary:                |    Pass    Total   Time
Ikeda Jacobian ≈ ForwardDiff | 2795100  2795100  45.7s

Symbolic verification with Wolfram Mathematica

Also, I verified it's correct by comparing it with Mathematica's result.

ClearAll[x, y, a, b, c, d, t, f, Jmat, Jhand, aux];

t[x_, y_, c_, d_] := c - d/(1 + x^2 + y^2);
f[x_, y_, a_, b_, c_, d_] := {
  a + b*(x*Cos[t[x, y, c, d]] - y*Sin[t[x, y, c, d]]),
  b*(x*Sin[t[x, y, c, d]] + y*Cos[t[x, y, c, d]])
};

Jmat = D[f[x, y, a, b, c, d], {{x, y}}] // FullSimplify; 

aux = 2 d/(1 + x^2 + y^2)^2;
Jhand = {
  { b*(Cos[t[x,y,c,d]] - aux*(x^2*Sin[t[x,y,c,d]] + x*y*Cos[t[x,y,c,d]])),
    b*(-Sin[t[x,y,c,d]] - aux*(x*y*Sin[t[x,y,c,d]] + y^2*Cos[t[x,y,c,d]])) },
  { b*( Sin[t[x,y,c,d]] + aux*(x^2*Cos[t[x,y,c,d]] - x*y*Sin[t[x,y,c,d]])),
    b*( Cos[t[x,y,c,d]] + aux*(x*y*Cos[t[x,y,c,d]] - y^2*Sin[t[x,y,c,d]])) }
};

(* Equality should return True *)
Simplify[Jmat == Jhand,
  Assumptions -> Element[{x,y,a,b,c,d}, Reals] && (1 + x^2 + y^2) > 0
]

(* Zero matrix confirmation *)
Simplify[Jmat - Jhand,
  Assumptions -> Element[{x,y,a,b,c,d}, Reals] && (1 + x^2 + y^2) > 0
] // MatrixForm

Output

True

$$ \begin{pmatrix} 0 & 0\\ 0 & 0 \end{pmatrix} $$

[Datseris]

thank you!