# Two first order non-linear differential equation: Using DSolve to obtain the solutions analytically

I have two simple first order (coupled) differential equations which I denote as system together with its initial conditions

system = {A1'[x] == (j1/2)*A1[x]^2*Conjugate[A1[x]] + j1*A1[x]*A2[x]*Conjugate[A2[x]],
A2'[x] == (j2/2)*A2[x]^2*Conjugate[A2[x]] + j2*A1[x]*A2[x]*Conjugate[A1[x]],
A1[0] == √[P10], A2[0] == √[P20];


where j1, j2, P10, P20 are all constants and independent of x. I then do

DSolve[system,{A1,A2},x]


And I'm returned with: There are fewer dependent variables than equations, so the system is overdetermined. But I have exactly 2 equations with 2 unknowns, so I'm not sure why Mathematica is returning this error. What am I doing wrong here?

Thank you.

• Is \Sqrt a typo or you really write this in your code? Commented Sep 7, 2020 at 2:06
• Which variables and constants, if any, are complex? Do you have reason to suppose that analytical solutions actually exist? DSolve is unable to find a solution, although it does not generate an error message for me. Commented Sep 7, 2020 at 3:56
• Try quitting the kernel and then retry your code. Commented Sep 7, 2020 at 5:38
• \Sqrt is a typo, but I didn't type that wrong in my code @xzczd Commented Sep 7, 2020 at 6:29
• I mentioned j1, j2, P10 and P20 are constants @bbgodfrey An analytical solution should exist given that it represents the pump amplitude of a physical system. Even if it doesn't it should return something that says unable to find solution. Not what im getting. Also quitting the kernel didn't change anything Commented Sep 7, 2020 at 6:32

## 1 Answer

Amazingly, this system can be solved symbolically, at least for j2 = j1. How useful the answer may be is a different issue. If all symbols are real, then system (without initial conditions) can be rewritten as

system = {p1'[x] == j1 p1[x]^2 + 2 j1 p1[x] p2[x],
p2'[x] == 2 j2 p1[x] p2[x] + j2 p2[x]^2}


where p1[x] = A1[x]^2 and p2[x] = A2[x]^2. (If the symbols are complex, then p1[x] = A1[x]*Conjugate[A1[x]], p2[x] = A2[x]*Conjugate[A2[x]], and j1 and j2 replaced by their real parts.) The obvious next step, DSolve[system /. j2 -> j1, {p1, p2}, x] returns unevaluated with unexpected error messages. So, instead, eliminate p2, which is straightforward.

D[Simplify[#/p1[x]], x] & /@ First[system];
% /. Rule @@ Last[system];
Solve[First[system], p2[x]] // Flatten;
eq = Simplify[%% /. %]
(* 3 j1 j2 p1[x]^3 + 2 p1''[x] ==
2 (j1 + j2) p1[x] p1'[x] + ((2 j1 + j2) p1'[x]^2)/(j1 p1[x]) *)


which can be solved.

s = DSolve[eq /. j2 -> j1, b1[x], x][[1, 1]]
(* p1[x] -> InverseFunction[Inactive[Integrate][4/(3 (4 j1 K[1]^2 +
(E^((3 C[1])/8) j1^2 K[1]^3)/(4 E^((3 C[1])/8) j1^3 K[1]^5 +
Sqrt[-E^(((9 C[1])/8)) j1^6 K[1]^9 + 16 E^((3 C[1])/4) j1^6 K[1]^10])^(1/3) +
(4 E^((3 C[1])/8) j1^3 K[1]^5 + Sqrt[-E^(((9 C[1])/8)) j1^6 K[1]^9 +
16 E^((3 C[1])/4) j1^6 K[1]^10])^(1/3))), {K[1], 1, #1}] &][x + C[2]] *)


Perhaps, more convenient is

s[[2, 1]] - C[2] -> s[[2, 0, 1]][s[[1]]] - C[2];
% /. {K[1], 1, p1[x]} -> {K[1], p10, p1[x]} /. C[2] -> 0
(* x -> Inactive[Integrate][4/(3 (4 j1 K[1]^2 +
(E^((3 C[1])/8) j1^2 K[1]^3)/(4 E^((3 C[1])/8) j1^3 K[1]^5 +
Sqrt[-E^(((9 C[1])/8)) j1^6 K[1]^9 + 16 E^((3 C[1])/4) j1^6 K[1]^10])^(1/3) +
(4 E^((3 C[1])/8) j1^3 K[1]^5 + Sqrt[-E^(((9 C[1])/8)) j1^6 K[1]^9 +
16 E^((3 C[1])/4) j1^6 K[1]^10])^(1/3))), {K[1], p10, p1[x]}] *)


p2 then can be determined by inserting the p1solution into system[[1]], and the initial condition p20 applied to obtain C[1]. If j2!= j1, DSolve returns unevaluated.

Addendum: Obtaining first integral

C[1] can be derived more simply by obtaining a first integral of eq by the standard transformation valid for autonomous second order ODEs. Doing so provides some useful insight as well.

eqv = eq /. {p1''[x] -> v'[z] v[z], p1'[x] -> v[z], p1[x] -> z}
(* 3 j1 j2 z^3 + 2 v[z] v'[z] == (j1 + j2) z v[z] + ((2 j1 + j2) v[z]^2)/(j1 z) *)

eqc1 = (DSolve[eqv, v[z], z] // First // Simplify) /. {v[z] -> p1'[x], z -> p1[x]}
(* (1/(3 j2)) 2 (j1 + j2)^2 (Log[p1[x]]/j1 +
(Log[(2 (j1 + j2) (j1 p1[x]^2 - p1'[x]))/(j1 p1[x]^2)] - (3 j2
Log[-((2 (j1 + j2) (3 j1 j2 p1[x]^2 + (-2 j1 + j2) p1'[x]))/(3 j1 j2 p1[x]^2))])
/(2 j1 - j2))/(j1 - 2 j2)) == C[1] *)


from which C[1] is obtained by setting x == 0 and applying initial conditions. To obtain the complete solution for p1[x] requires solving eqc1 for p1'[x], which can be done for j2 = j1 and, perhaps, a few other cases. (This explains why DSolve[eq, b1[x], x] is unable to obtain solutions for arbitrary choices of j2/j1.) So, consider j2 = j1.

Simplify[eqc1 /. {p1'[x] -> j1 p10^2 + 2 j1 p10  p20, p1[x] -> p10}]
Simplify[% /. j2 -> j1]
Exp[3/8 #] & /@ %
(* -((64 j1^2 (-p10 + p20)^3)/(27 p10 p20)) == E^((3 C[1])/8) *)


This procedure also provides a simple evaluation of p2[x]. Interchange p1[x] and p2[x], and p10 and p20 in the code block just before this addendum and also in the expression for C[1] just above.