2
$\begingroup$

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.

Improve this question
$\endgroup$
6
  • 1
    $\begingroup$ Is \Sqrt a typo or you really write this in your code? $\endgroup$
    – xzczd ♦
    Commented Sep 7, 2020 at 2:06
  • $\begingroup$ 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. $\endgroup$
    – bbgodfrey
    Commented Sep 7, 2020 at 3:56
  • $\begingroup$ Try quitting the kernel and then retry your code. $\endgroup$
    – Michael E2
    Commented Sep 7, 2020 at 5:38
  • $\begingroup$ \Sqrt is a typo, but I didn't type that wrong in my code @xzczd $\endgroup$
    – kowalski
    Commented Sep 7, 2020 at 6:29
  • $\begingroup$ 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 $\endgroup$
    – kowalski
    Commented Sep 7, 2020 at 6:32

1 Answer 1

Reset to default
4
$\begingroup$

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.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.