Integrating and plotting solutions of a nonlinear dynamical system

Question

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?

Answer 1

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<a<b$ then $\int_0^z \frac{1}{\sqrt{G(u)}} du\;$ where $0<z<a\;$ can be expressed in terms of logarithm:

Integrate[ (Sqrt[3] 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]] == 4 && cfl[[3]] == 0, {a, b}] // Flatten

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

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

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

{{w[θ] -> WeierstrassP[ θ - C[1], {1/12 (-2B + ω^2/k^6), 
                       -((3A k^11 - 3B k^6 ω + ω^3)/(216k^9))}]},
 {w[θ] -> WeierstrassP[ θ + C[1], {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)[[1]] // 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[[1]] //. {C[1] -> 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$.

Answer 2

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]

This solution, and your attempt separation of variables too, need additional knowledge about the parameters!