# Jacobi Elliptic Function Solution using DSolve

I am happy that I have a solution to a particular differential equation by hand:

eq = 2 * D[psi[x], {x,1}]^2 == (psi[x]^2 + f0^2)*(2(1-f0^2) - psi[x]^2) + (2*ju^2/f0^2)


We can write this as

eq = 2 * D[psi[x], {x,1}]^2 == (an^2 + psi[x]^2)(bn^2 - psi[x]^2)


where

an = Sqrt[- ( 1 - (3/2) f0^2) + Sqrt[(1 - (f0^2)/2)^2 + 2 ju^2/f0^2] ];
bn = Sqrt[+ ( 1 - (3/2) f0^2) + Sqrt[(1 - (f0^2)/2)^2 + 2 ju^2/f0^2] ];

mn = (bn^2) / (an^2 + bn^2);
m1n = 1 - mn //FullSimplify;

solpsi[x_] = Sqrt[an^2 * mn] * JacobiSD[x*Sqrt[(an^2 + bn^2)/2], mn ];


This solution can be verified:

JEFrules = {JacobiCD[x_,mn_] -> Sqrt[1 - (1-mn)*JacobiSD[x, mn]^2]
,JacobiND[x_,mn_] -> Sqrt[1 + mn*JacobiSD[x, mn]^2] }
lhs = 2 * D[solpsi[x], {x,1}]^2/.JEFrules;
rhs = (an^2 + solpsi[x]^2)(bn^2 -solpsi[x]^2);
lhs - rhs //Simplify


with output equal to zero as expected.

If I now pretend that I don't know the solution, I can't obtain this solution solpsi using DSolve. How can I solve my original differential equation, without needing to do the subsequent hand manipulations with MMA? Do I need to use something other than DSolve?

For ease: my attempt to solve directly using either

eq = 2 * D[psi[x], {x,1}]^2 == (psi[x]^2 + f0^2)*(2(1-f0^2) - psi[x]^2) + (2*ju^2/f0^2);
sol = DSolve[{eq, psi[0] == 0}, psi[x], x]


or

eq = 2 * D[psi[x], {x,1}]^2 == (an^2 + psi[x]^2)(bn^2 - psi[x]^2)
sol = DSolve[{eq, psi[0] == 0}, psi[x], x]


both provide very weird solutions, which certainly aren't as pleasant as the above solpsi; I'm not even sure they are equivalent when they clearly should be!

A final example, with output included:

eq = Sqrt[2]*D[psi[x], {x,1}] == Sqrt[(an^2 + psi[x]^2)(bn^2 - psi[x]^2)]
sol = DSolve[{eq, psi[0] == 0}, psi[x], x]


Output:

{{psi[x] -> InverseFunction[(I*EllipticF[I*ArcSinh[Sqrt[2]*Sqrt[f0/(-2*f0 + 3*f0^3 + Sqrt[4*f0^2 - 4*f0^4 + f0^6 + 8*ju^2])]*#1], -((-2*f0 + 3*f0^3 + Sqrt[4*f0^2 - 4*f0^4 + f0^6 + 8*ju^2])/(2*f0 - 3*f0^3 + Sqrt[4*f0^2 - 4*f0^4 + f0^6 + 8*ju^2]))]*Sqrt[1 - (2*f0*#1^2)/(2*f0 - 3*f0^3 + Sqrt[4*f0^2 - 4*f0^4 + f0^6 + 8*ju^2])]*Sqrt[1 + (2*f0*#1^2)/(-2*f0 + 3*f0^3 + Sqrt[4*f0^2 - 4*f0^4 + f0^6 + 8*ju^2])])/(Sqrt[2]*Sqrt[f0/(-2*f0 + 3*f0^3 + Sqrt[4*f0^2 - 4*f0^4 + f0^6 + 8*ju^2])]*Sqrt[-2*f0^4 + (2*ju^2)/f0^2 + f0^2*(2 - 3*#1^2) - #1^2*(-2 + #1^2)]) & ][x/Sqrt[2]]}}


Clearly not solpsi! Any advice appreciated.

• Spell psi/Psi consistently? ;-) Apr 23, 2021 at 12:10
• @MichaelE2 Gah! Curse my attempt to create a MWE... I'll correct it, but the problem still stands I think! Thanks for your sharp eyes
Apr 23, 2021 at 12:20
• I get bn JacobiSN[(an x)/Sqrt[2], -(bn^2/an^2)], which is certainly simple. I didn't check equivalence though. Apr 23, 2021 at 12:39
• That's very odd... May I ask for you to provide a snippet of your code? Or have you just copied and pasted my lines above? I have added my typical outputs above, in terms of InverseFunctions, EllipticF and ArcSinh... oh, I think you've done DSolve without using the an or bn expressions right?
Apr 23, 2021 at 12:42
• eq = 2*D[psi[x], {x, 1}]^2 == (an^2 + psi[x]^2) (bn^2 - psi[x]^2); sol = DSolve[{eq, psi[0] == 0}, psi, x] -- it was one of the codes you posted, but with the spelling fix. Apr 23, 2021 at 12:46

Thanks to Michael for starting this approach. Without specifying an or bn, we obtain a solution like:

bn JacobiSN[(an x)/Sqrt[2], -(bn^2/an^2)]


This is secretly equivalent to the psisol form as desired above, by using an RIR transform, described in the wiki page and some textbooks - relevant section of wiki page linked here..

It is then possible to manipulate the SN function into the SD form required. A priori, I am not sure if MMA can find this form, so I will leave this question open for the time being.

EDIT; from the solution above, it's actually quite easy to algebraically solve for the required forms of mn and m1n which exactly match our expected solution. Problem solved with a few equations solving each of the arguments!

• So what is the exact transformation? Apr 24, 2021 at 12:07
• It is the transform linked in the table; pq(u,m) = gamma_{pq} (pq)' (k' u, -m/m'). One can identify the relevant result from this (where we relate sd to sn using this). However, obtaining this directly from MMA is still challenging.
• Of course you are correct, one could use a manual transformation. The real challenge I think is identifying what m and m' are here though from this form. They should be the same as the mn and m1n as defined in my original statement - in traditional JEF notation, here k^2 = m and k'^2 = m'. It is trivial on paper when ones knows what the forms are; it is not so clear to me using MMA directly!