# Solving a nonlinear system of ODE with DSolve

I have been trying to solve the following system of ODE with initial values using the function DSolve. My code is :

system = {x''[t] == ((((B0 k) x[t] + B0) q)/m) y'[t],  y''[t] == -((((B0 k)
x[t] + B0) q)/m) x'[t]};
initialvalues2 = {x[0] == 0, y[0] == 0, x'[0] == 10^4, y'[0] == 0};
solution = DSolve[Join[system, initialvalues2], {x[t], y[t]}, t]


B0, k, q and m are constants defined later in the code.

However, Mathematica keeps returning the system without providing any answer. Any help would be appreciated :)

• It may be of use to note that Derivative[1][x][t]^2 + Derivative[1][y][t]^2 is a constant of the motion. – bbgodfrey Oct 26 '16 at 21:04
• This is a nonlinear system. – Michael E2 Oct 26 '16 at 23:04

To elaborate on my comment, considerable progress can be made as follows. First, the constant of the motion given in my comment is derived by

system = Map[Factor, system, Infinity] /. B0 q/m -> b
(* {Derivative[2][x][t] == b (1 + k x[t]) Derivative[1][y][t],
Derivative[2][y][t] == -b (1 + k x[t]) Derivative[1][x][t]} *)


to combine multiple constants. Then,

Simplify[Subtract @@@ system.{x'[t], y'[t]}];
Simplify[2 Integrate[%, t]];
% == (% /. t -> 0 /. {Derivative[1][x][0] -> 10000, Derivative[1][y][0] -> 0})
(* Derivative[1][x][t]^2 + Derivative[1][y][t]^2 == 100000000 *)


Then, as noted by Alexei Boulbitch, the two ODEs in system are readily combined to give

system[[1]] /. y'[t] -> Sqrt[100000000 - Derivative[1][x][t]^2]
(* Derivative[2][x][t] == b (1 + k x[t]) Sqrt[100000000 - Derivative[1][x][t]^2] *)


However, at this point we shift the zero point of x[t] by 1/k and only then solve the resulting equation. (Of course, the initial condition on x[0] also must be shifted to 1/k.)

Flatten@DSolve[Derivative[2][x][t] == b k x[t] Sqrt[100000000 - Derivative[1][x][t]^2],
x[t], t]
(* x[t] -> InverseFunction[-((2 I EllipticF[I ArcSinh[Sqrt[(b k)/(20000 + 2 C[1])] #1],
(10000 + C[1])/(-10000 + C[1])] Sqrt[2 + (b k #1^2)/(-10000 + C[1])]
Sqrt[1 + (b k #1^2)/(20000 + 2 C[1])])/
(Sqrt[(b k)/(10000 + C[1])] Sqrt[400000000 - 4 C[1]^2 - 4 b k C[1] #1^2
- b^2 k^2 #1^4])) &][t + C[2]] *)


plus a second solution with is the negative of the first. They are much simpler expressions.

I should add that the two constants C can be evaluated after inverting the solution above, applying the boundary conditions, and performing a bit of algebra.

Alternative Derivation with Constants Determined

From my perspective, it is simpler and more instructive to solve the ODE using Integrate twice, each time evaluating the resulting constant of integration. Begin with

eq = (#/y'[t]) & /@ Simplify[system[[1]] /. x[t] -> x[t] - 1/k] /.
y'[t] -> Sqrt[100000000 - Derivative[1][x][t]^2]
(* b k x[t] == Derivative[2][x][t]/Sqrt[100000000 - Derivative[1][x][t]^2] *)


Note that the zero point of x[t] is shifted by 1/k, as above. Doing so enormously simplifies subsequent results. Next Integrate eq and add a constant of integration (equivalent to C[1] above).

i1 = Integrate[# x'[t], t] & /@ eq;
i1[[1]] = i1[[1]] + c1; i1
(* c1 + 1/2 b k x[t]^2 == -Sqrt[100000000 - Derivative[1][x][t]^2] *)


Now, apply the boundary conditions to determine and eliminate c1.

Flatten@Solve[i1, c1] /. t -> 0 /. {x'[0] -> 10000, x[0] -> 1/k}
(* {c1 -> -(b/(2 k))} *)
i1 = i1 /. %
(* -(b/(2 k)) + 1/2 b k x[t]^2 == -Sqrt[100000000 - Derivative[1][x][t]^2] *)


Next, Solve this equation for x'[t], so that Integrate can again be applied. (For now, consider only the first solution from Solve.)

eq2 = Equal @@ (Solve[i1, x'[t]][[1, 1]]);
eq2 = eq2[[1]]/eq2[[2]] == 1
(* -((2 k Derivative[1][x][t])/
Sqrt[-b^2 + 400000000 k^2 + 2 b^2 k^2 x[t]^2 - b^2 k^4 x[t]^4]) == 1 *)
i2 = Reverse[Integrate[# , t] & /@ eq2];
i2[[2]] = i2[[2]] + c2; i2
(* t == c2 + (2 I k EllipticF[I ArcSinh[Sqrt[-((b k^2)/(b - 20000 k))] x[t]],
(b - 20000 k)/(b + 20000 k)] Sqrt[1 - (b k^2 x[t]^2)/(b - 20000 k)]
Sqrt[1 - (b k^2 x[t]^2)/(b + 20000 k)])/(Sqrt[-((b k^2)/(b - 20000 k))]
Sqrt[-b^2 + 400000000 k^2 + 2 b^2 k^2 x[t]^2 - b^2 k^4 x[t]^4]) *)


A substantial simplification now can be obtained by using the identity,

Sqrt[1 - (b k^2 x[t]^2)/(b - 20000 k)] Sqrt[1 - (b k^2 x[t]^2)/(b + 20000 k)] ==
Sqrt@Expand[-(b - 20000 k - b k^2 x[t]^2) (b + 20000 k -
b k^2 x[t]^2)] Sqrt[-1/(b - 20000 k)] Sqrt[1/(b + 20000 k)]


Applying it yields

(* t == c2 + (2 I EllipticF[I ArcSinh[Sqrt[-(b/(b - 20000 k))] k x[t]],
(b - 20000 k)/(b + 20000 k)])/Sqrt[b (b + 20000 k)] *)


Finally, determine and eliminate c2, and shift the zero point of x[t] back.

Simplify[Flatten@Solve[i2 /. t -> 0 /. x[0] -> 1/k, c2], k > 0 && b > 0]
(* {c2 -> -((2 I EllipticF[I ArcSinh[Sqrt[-(b/(b - 20000 k))]], (b - 20000 k)/
(b + 20000 k)])/Sqrt[b (b + 20000 k)])} *)
i2 = Simplify[i2 /. % /. x[t] -> x[t] + 1/k]
(* t == -((2 I (EllipticF[I ArcSinh[Sqrt[-(b/(b - 20000 k))]],
(b - 20000 k)/(b + 20000 k)] -
EllipticF[I ArcSinh[Sqrt[-(b/(b - 20000 k))] (1 + k x[t])],
(b - 20000 k)/(b + 20000 k)]))/Sqrt[b (b + 20000 k)]) *)


(Note that the second solution from Solve, alluded to earlier, results in the negative of the expression just obtained.) These results, although cumbersome, do provide analytical expressions for the limits on x[t] and for its oscillation period.

xlim = Solve[b/(b - 20000 k) k^2 x[t]^2 == (b + 20000 k)/(b - 20000 k), x[t]]
(* {{x[t] -> -(Sqrt[b + 20000 k]/(Sqrt[b] k))}, {x[t] -> Sqrt[b + 20000 k]/(Sqrt[b] k)}} *)

period = Simplify[2 Simplify[Subtract @@ (i2[[2]] /. x[t] -> x[t] - 1/k /. %)],
k > 0 && b > 0]
(* -((8 I EllipticF[I ArcSinh[Sqrt[-((b + 20000 k)/(b - 20000 k))]],
(b - 20000 k)/(b + 20000 k)])/Sqrt[b (b + 20000 k)]) *)


To illustrate these results graphically, consider the case {b -> 1, k -> 1}.

b1k1 = i2 /. {b -> 1, k -> 1} /. x[t] -> x;
shift = % /. {b -> 1, k -> 1};
Show[ParametricPlot[{shift # + (2 Mod[#, 2] - 1) b1k1[[1]], x},
{x, -1 - Ceiling[#, 6] Sqrt[20001], -1 + Sqrt[20001]},
AspectRatio -> 1/GoldenRatio, AxesLabel -> {t, x}] & /@
Range[0, 5], PlotRange -> {{0, .2}, Automatic}]


Numerical Solution

The corresponding numerical solution can be obtained without difficulty.

system = {x''[t] == ((((B0 k) x[t] + B0) q)/m) y'[t],
y''[t] == -((((B0 k) x[t] + B0) q)/m) x'[t]} /. {B0 -> 1, q -> 1, m -> 1, k -> 1};
initialvalues2 = {x[0] == 0, y[0] == 0, x'[0] == 10^4, y'[0] == 0};
solution = Flatten@NDSolve[{system, initialvalues2}, {x[t], y[t], x'[t], y'[t]},
{t, 0, 2}];


Plotting x[t] /. solution gives a figure identical to that obtained from the analytical solution above. For completeness, all four numerical curves are given by

Plot[Evaluate[{x[t], y[t]} /. solution], {t, 0, .2}]


Plot[Evaluate[{x'[t], y'[t]} /. solution], {t, 0, .2}]


I would close by observing that, while an analytical solution indeed can be obtained, a numerical solution is easier to derive and probably is more useful.

Making use of the note of @bbgodfrey we easily reduce this system to a single equation. I have chosen the case with the sign + in front of the radical:

eq = x''[t] == ((((B0 k) x[t] + B0) q)/m) *Sqrt[10^8 - x'[t]^2];


Then DSolve returns two solutions, which are rather heavy. I write here the first of them:

DSolve[eq, x[t], t]

(*  x[t] -> InverseFunction[-((4 I m (Sqrt[-B0 q (20000 k m + B0 q +
2 I k m C[1])] +
Sqrt[-B0 q (B0 q + 2 I k m (10000 I + C[1]))]) EllipticF[
ArcSin[\[Sqrt]((I (I B0^2 q^2 +
I Sqrt[-B0 q (20000 k m + B0 q + 2 I k m C[1])]
Sqrt[-B0 q (B0 q + 2 I k m (10000 I + C[1]))] +
B0 q (-Sqrt[-B0 q (20000 k m + B0 q +
2 I k m C[1])] +
Sqrt[-B0 q (B0 q + 2 I k m (10000 I + C[1]))] +
k (20000 I m - 2 m C[1] -
Sqrt[-B0 q (20000 k m + B0 q + 2 I k m C[1])] #1 +
Sqrt[-B0 q (B0 q +
2 I k m (10000 I +
C[1]))] #1))))/((Sqrt[-B0 q (20000 k m + B0 q +
2 I k m C[1])] +
Sqrt[-B0 q (B0 q +
2 I k m (10000 I + C[1]))]) (Sqrt[-B0 q (B0 q +
2 k m (10000 + I C[1]))] +
I B0 (q + k q #1))))], (Sqrt[-B0 q (20000 k m +
B0 q + 2 I k m C[1])] +
Sqrt[-B0 q (B0 q +
2 I k m (10000 I + C[1]))])^2/(Sqrt[-B0 q (20000 k m +
B0 q + 2 I k m C[1])] -
Sqrt[-B0 q (B0 q +
2 I k m (10000 I + C[1]))])^2] (Sqrt[-B0 q (B0 q +
2 k m (10000 + I C[1]))] +

I B0 (q +
k q #1))^2 \[Sqrt]((Sqrt[-B0 q (20000 k m + B0 q +
2 I k m C[1])] (Sqrt[-B0 q (B0 q +
2 I k m (10000 I + C[1]))] +
I B0 (q + k q #1)))/((Sqrt[-B0 q (20000 k m + B0 q +
2 I k m C[1])] +
Sqrt[-B0 q (B0 q +
2 I k m (10000 I + C[1]))]) (Sqrt[-B0 q (B0 q +
2 k m (10000 + I C[1]))] +
I B0 (q +
k q #1)))) \[Sqrt]((I Sqrt[-B0 q (20000 k m +
B0 q + 2 I k m C[
1])] (I Sqrt[-B0 q (B0 q +
2 I k m (10000 I + C[1]))] +
B0 (q + k q #1)))/((Sqrt[-B0 q (20000 k m + B0 q +
2 I k m C[1])] -
Sqrt[-B0 q (B0 q +
2 I k m (10000 I + C[1]))]) (Sqrt[-B0 q (B0 q +
2 k m (10000 + I C[1]))] +
I B0 (q + k q #1)))) \[Sqrt]((I (I B0^2 q^2 +
I Sqrt[-B0 q (20000 k m + B0 q + 2 I k m C[1])]
Sqrt[-B0 q (B0 q + 2 I k m (10000 I + C[1]))] +
B0 q (-Sqrt[-B0 q (20000 k m + B0 q + 2 I k m C[1])] +
Sqrt[-B0 q (B0 q + 2 I k m (10000 I + C[1]))] +
k (20000 I m - 2 m C[1] -
Sqrt[-B0 q (20000 k m + B0 q + 2 I k m C[1])] #1 +
Sqrt[-B0 q (B0 q +

2 I k m (10000 I +
C[1]))] #1))))/((Sqrt[-B0 q (20000 k m + B0 q +
2 I k m C[1])] +
Sqrt[-B0 q (B0 q +
2 I k m (10000 I + C[1]))]) (Sqrt[-B0 q (B0 q +
2 k m (10000 + I C[1]))] +
I B0 (q + k q #1)))))/(B0 k q Sqrt[-B0 q (B0 q +
2 k m (10000 + I C[1]))] (Sqrt[-B0 q (20000 k m + B0 q +
2 I k m C[1])] -
Sqrt[-B0 q (B0 q + 2 I k m (10000 I + C[1]))]) Sqrt[
4 m^2 (100000000 + C[1]^2) + 4 I B0 m q C[1] #1 (2 + k #1) -
B0^2 q^2 #1^2 (2 + k #1)^2])) &][t + C[2]]  *)


It is so heavy that I doubt that it can be useful.

The solution depends on 2 arbitrary constants. However, I do not think that it will be possible to solve the initial conditions with respect to C[1] and C[2]. It is for this reason probably that Mma returns no answer, if one applies the DSolve with initial conditions:

DSolve[{eq, x[0] == 0, x'[0] == 10^4}, x[t], t]


Have fun!