Ykhjuy

I have a question to a physical task in Mathematica.

We have this equation of motion:

$$mcdotddot{x} = -m(omega_0^2cdot x+(epsilon x^3))=-frac{d}{dx}V(x)$$

For energy of masspoint there is the condition :

$$epsilon Ell momega_0^4$$

I have to write a procedure that uses the law of the conservation of energy for the potential $V(x)$ to calculate $t(x_1) - t(x_0)$ when there are given two points $x_0$ and $x_1.

How could I do this in Mathematica?

In a numerical model, energy is conserved with some accuracy; in this example, the deviation from the initial value is about 1.5*10^-12

m = 1; omega0 = 1; eps = 1/100; v0 = 1;

eq = m*x''[t] == -m*(omega0^2*x[t] + eps*x[t]^3);

ic = {x[0] == 0, x'[0] == v0};



X = NDSolveValue[{eq, ic}, x, {t, 0, 10}, WorkingPrecision -> 30];



Plot[m/2*X'[t]^2 + m/2*omega0^2*X[t]^2 + m/4*eps*X[t]^4 - 

  m/2*v0^2, {t, 0, 10},AxesLabel -> {"t", "E-E0"}]

Using the law of conservation of energy, we express $x'(t) $ and then time as a function of $x$

t=Integrate[1/Sqrt[v0^2 - omega0^2*x^2 - eps/2*x^4], x]

(*-((I Sqrt[2 + (2 eps x^2)/(omega0^2 - Sqrt[omega0^4 + 2 eps v0^2])]

   Sqrt[1 + (eps x^2)/(omega0^2 + Sqrt[omega0^4 + 2 eps v0^2])]

   EllipticF[

   I ArcSinh[Sqrt[eps/(omega0^2 - Sqrt[omega0^4 + 2 eps v0^2])] x], (

   omega0^2 - Sqrt[omega0^4 + 2 eps v0^2])/(

   omega0^2 + Sqrt[omega0^4 + 2 eps v0^2])])/(

 Sqrt[eps/(omega0^2 - Sqrt[omega0^4 + 2 eps v0^2])] Sqrt[

  2 v0^2 - 2 omega0^2 x^2 - eps x^4]))*)

This problem can be solved symbolically as follows. Multiply the expression (m (omega0^2 x[t] + eps x[t]^3) + m x''[t]) by x'[t] and integrate to obtain an expression for the energy of this nonlinear oscillator.

eq = Integrate[(m (omega0^2 x[t] + eps x[t]^3) + m x''[t]) x'[t], t]

(* 1/2 m omega0^2 x[t]^2 + 1/4 eps m x[t]^4 + 1/2 m x'[t]^2 *)

with constant of integration v0, the conserved energy. Then, apply DSolve.

s = DSolve[eq == v0, x[t], t] // Last

(* {x[t] -> InverseFunction[-((I EllipticF[I ArcSinh[Sqrt[(eps m)/(m omega0^2 -

    Sqrt[m (m omega0^4 + 4 eps v0)])] #1], (m omega0^2 - Sqrt[m (m omega0^4 + 

    4 eps v0)])/(m omega0^2 + Sqrt[m (m omega0^4 + 4 eps v0)])] Sqrt[1 + (eps m #1^2)

    /(m omega0^2 - Sqrt[m (m omega0^4 + 4 eps v0)])] Sqrt[1 + (eps m #1^2)

    /(m omega0^2 + Sqrt[m (m omega0^4 + 4 eps v0)])])/(Sqrt[(eps m)/(m omega0^2 - 

    Sqrt[m (m omega0^4 + 4 eps v0)])] Sqrt[4 v0 - m #1^2 (2 omega0^2 + eps #1^2)])) &]

    [t/(Sqrt[2] Sqrt[m]) + C[1]]} *)

(The other solution is the negative of the first.) Since the question requests t as a function of x, s must be inverted. In the absence of a Mathematica command to accomplish this, we use the following ungainly expression.

st = Rule[(s[[1, 2, 1]] /. C[1] -> 0) Sqrt[2] Sqrt[m], 

    Head[s[[1, 2]]][[1]][x[t]] Sqrt[2] Sqrt[m]]

(* t -> -((I Sqrt[2] Sqrt[m] EllipticF[I ArcSinh[

       Sqrt[(eps m)/(m omega0^2 - Sqrt[m (m omega0^4 + 4 eps v0)])] x[t]], 

       (m omega0^2 - Sqrt[m (m omega0^4 + 4 eps v0)])/

       (m omega0^2 + Sqrt[m (m omega0^4 + 4 eps v0)])] 

       Sqrt[1 + (eps m x[t]^2)/(m omega0^2 - Sqrt[m (m omega0^4 + 4 eps v0)])] 

       Sqrt[1 + (eps m x[t]^2)/(m omega0^2 + Sqrt[m (m omega0^4 + 4 eps v0)])])/(

       Sqrt[(eps m)/(m omega0^2 - Sqrt[m (m omega0^4 + 4 eps v0)])] 

       Sqrt[4 v0 - m x[t]^2 (2 omega0^2 + eps x[t]^2)])) *)  

This result for various values of eps can be plotted as

st /. {m -> 1, omega0 -> 1, v0 -> 1};

Plot[Evaluate@Table[Last[%], {eps, {10^1, 1, 10^-1, 10^-10}}], {x[t], -2, 2}, 

    AxesLabel -> {x, t}, AspectRatio -> 1, ImageSize -> Large, 

    LabelStyle -> {Bold, Black, 15}]

Decreasing eps corresponds to increasing values of x and t at the turning points.