# Integrating and plotting solutions of a nonlinear dynamical system

So I have an ODE

$$k^2\; u'(\theta)^2 = -\frac{1}{3}u^3 + \frac{\omega}{k}u^2-2B\;u+2A$$

where $$B,A$$ are constants.

I attempted to directly integrate in Mathematica with

G = k^2 D[u,theta]^2 = -1/3 u^3 + omega/k  u^2 - 2B u + 2A
I = Integrate[k^2 D[u,theta]^2 = -1/3 u^3 + omega/k u^2 - 2B u + 2A,u]


But got no luck. The official solution says that I should get

G = 1/3 ( (u1 - u) (u2 - u) (u3 - u) )


As my solution. I'm them asked to integrate $$1/\sqrt{G(u)}$$ as

Integrate[ Sqrt[G], {u, u1, u2}]


But again, got nowhere. How might I also plot G? say for u1 = 1, u2 = 2, u3 = 1?

• You can get e.g. Integrate[1/Sqrt[(1 - u)^2 (2 - u)], u] yielding -2 I ArcTan[Sqrt[-2 + u]] if u>2, however in general if u1 != u3 i.e. (u1 not equal to u3 ) this will be an elliptic integral, then u[theta] will be an elliftic function expressible in terms of e.g. the Weierstrass elliptic function $\wp$ see e.g. this post 1. If you expect more detailed help you should ask your question more precisely. Apr 1, 2022 at 13:26
• If constants k, A, B and $\omega$ are symbolic you cannot prescribe u1, u2, u3 but they can be expressed in terms of k, A, B and $\omega$. One can plot a general solution only if an initial condition is given. Apr 1, 2022 at 13:29
• I've explained how you can solve your equation here. I've mentioned it in comments to your question. Apr 1, 2022 at 14:42
• There is a difference betwen = and ==. The double == is for representing an equation. Apr 1, 2022 at 15:27
• @Roberto_1986 Explain meaning of symbols $ω,k,A,B$. What does this ode describe? What are interesting ranges of values of the constants? What is domain of $\theta$? Apr 14, 2022 at 1:54

Special case

There is a double root of $$G$$, i.e. $$G(a)=0$$. $$\frac{1}{\sqrt{G(u)}}=\frac{\sqrt{3} k}{\sqrt{(a-u)(a-u)(b-u)}}$$ Then the integral $$\int_a^b \frac{1}{\sqrt{G(u)}} du$$ does not converge. If $$\;0 then $$\int_0^z \frac{1}{\sqrt{G(u)}} du\;$$ where $$0 can be expressed in terms of logarithm:

Integrate[ (Sqrt k)/Sqrt[(a - u) (b - u) (a - u)], {u, 0, z},
Assumptions -> 0 < z < a < b]//TraditionalForm General case - solution in terms of special functions

There are three different roots of $$G(u)$$

Since the original question haven't mentioned anything regarding a physical model, interesting ranges of parameters, domain of variables etc. we are to guess what we could expect from a mathematical analysis. We can observe that our differential equation describes an elliptic curve in the phase space $$(u(\theta),\frac{d u(\theta)}{d \theta})$$. Ususally parameters $$k$$ and $$\omega$$ denote a spring (or elasticity) constant and an angular velocity (or frequency) respectively. The ODE may describe a nonlinear (anharmonic) oscillator and so we are rather interested in periodic solutions. A convenient tool for visualisation of the system behaviour is Manipulate and we can now demonstrate what happens is the phase space, here up stands for u'[θ] and u for u[θ]:

Manipulate[
ContourPlot[
k^2 up^2 == 2A - 2B u + ω u^2/k - u^3/3 , {u, -10, 10}, {up, -10, 10},
ContourStyle -> Thick, Axes -> True, AxesLabel -> TraditionalForm /@ {u[θ], u'[θ]}],
{{A, 1/4}, -4, 1}, {{B, 1}, -6, 2}, {{k, 2/3}, 0, 2}, {{ω, 5/4}, 1/2, 2}] Here the unbounded curve is related to singular solutions while the closed curve is related to periodic solutions. Now we would like to parametrize these curves solving the ODE.

In order to solve this differential equation, let's recast it into the canonical Weierstrass form $$w'(\theta)^2=4w(\theta)^3-g_2 w(\theta)-g_3\quad$$ with help of a linear transformation $$u(\theta)=a\; w(\theta)+b$$. Solving appropriate algebraic equations we will get coefficients $$a$$ and $$b$$, while $$g_2$$ and $$g_3$$ can be easily expressed in terms of the given constants. This procedure had been applied e.g. to solve geodesics equations in the Schwarzschild space-time and it was clearly expressed in comments under earlier question. Dividing the original ode by $$k^2 a^2$$ we construct an appropriate coefficient list and solving imposed conditions we get $$a$$ and $$b$$:

cfl = CoefficientList[1/a^2 (2A/k^2 - 2B u[θ]/k^2 + ω u[θ]^2/k^3 - u[θ]^3/3k^2
) /. u[θ] -> a w[θ] + b, w[θ]];
sol = Solve[ cfl[] == 4 && cfl[] == 0, {a, b}] // Flatten

{{a -> -12 k^2, b -> ω/k}}

{g2, g3} = {-cfl[], -cfl[]};

dsol = DSolve[ w'[θ]^2 == 4w[θ]^3 - g2 w[θ] - g3, w[θ], t]/.sol //Simplify

{{w[θ] -> WeierstrassP[ θ - C, {1/12 (-2B + ω^2/k^6),
-((3A k^11 - 3B k^6 ω + ω^3)/(216k^9))}]},
{w[θ] -> WeierstrassP[ θ + C, {1/12 (-2B + ω^2/k^6),
-((3A k^11 - 3B k^6 ω + ω^3)/(216 k^9))}]}}


and now we obtain the solution (the both solutions differ by the sign of the integration constant i.e. $$\theta \pm c_1$$ and so we write down only the first one)

(usol = a w[θ] + b /. sol /. dsol)[] // TraditionalForm Since the solutions are expressed by the Weierstrass elliptic function $$\wp$$ they involve singularities in a finite independent variable θ (assuming that $$A, B, k,\omega\;$$ are real) unless we impose certain complex parameter $$c_1$$.

We define half-periods:

{WHP1, WHP2, WHP3} = (Through @ { WeierstrassHalfPeriodW1, WeierstrassHalfPeriodW2,
WeierstrassHalfPeriodW3}@{g2, g3}) /. sol;


and roots $$e_1,e_2,e_3$$ of the Weierstrass polynomial:

{WE1, WE2, WE3} =
Through @ { WeierstrassE1, WeierstrassE2, WeierstrassE3}@{g2, g3} /. sol;


now we can factorize the Weierstrass polynomial $$w'(\theta)^2=4w(\theta)^3-g_2 w(\theta)-g_3=4(w(\theta)-e_1)(w(\theta)-e_2)(w(\theta)-e_3)$$

FullSimplify[ WeierstrassP[{WHP1, WHP2, WHP3}, {g2, g3}] == {WE1, WE2, WE3} /. sol]

True

FullSimplify[ 4w[θ]^3 - g2 w[θ] - g3 == 4(w[θ] - WE1)(w[θ] - WE2)(w[θ] - WE3) /. sol /.
{A -> 1/4, B -> 1, k -> 2/3, ω -> 5/4} // N]

True


Let's plot an example:

Column[{{WHP1, WHP2, WHP3}, {WE1, WE2, WE3}} /.
{A -> 1/4, B -> 1, k -> 2/3, ω -> 5/4} // N // Chop]

 {0. - 1.57585 I, -3.62049 + 1.57585 I, 3.62049}
{-0.662408, 0.325374, 0.337034}

Plot[ReIm @ usol[] //. {C -> 1.575 I, A -> 1/4, B -> 1,
k -> 2/3, ω -> 5/4},
{ θ, -12, 12}, Evaluated -> True, PlotRange -> All] The solution takes the minimum and maximum values at $$\theta = m WHP3$$, where $$m$$ is integer, and its period is $$2 WHP$$.

To solve differential equation you need Mathematica functions DSolve(symbolical) or NDSolve(numerical)

The corrected form of your ode is

ode = k^2 u'[theta]^2 == -1/3*u[theta]^3 + omega/k*u[theta]^2 -2 B u[theta] + 2 A


general solution

U=DSolve[ode, u, theta]