How to solve these two differential equations

Question

I have two differential equations and I try to use function DSolve to solve them together.

DSolve[{r^2*A''[r] + 2*r*A'[r] - 2 A[r] + lambda^2*r^2*A[r] + 
2 lambda*r^2*dd'[r] == 0, 
r*dd''[r] + 2 dd'[r] + 3 lambda^2*r*dd[r] + 2 lambda*r*A'[r] + 
4 lambda*A[r] == 0}, {A, dd}, r]

There are two functions as A and dd. But it failed.

I try to use Maple software to solve it and it successes and needs little time.

The Maple code:

sol3 := r^2*diff(A(r), r, r) + 2*r*diff(A(r), r) - 2*A(r) + lambda^2*r^2*A(r) +2*lambda*r^2*diff(dd(r), r) = 0

sol4 := r*diff(dd(r), r, r) + 2*diff(dd(r), r) + 3*lambda^2*r*dd(r) + 2*lambda*r*diff(A(r), r) + 4*lambda*A(r) = 0

dsolve({sol3, sol4})

Therefore, the equations has solutions. I want to know whether Mathematica can solve it and how to do it. Thank you very much.

Answer 1

Edited to derive simpler expression for dd

The symbolic solution requested in the Question is obtained as follows.

eq1 = Expand[(r^2*A''[r] + 2*r*A'[r] - 2 A[r] + lambda^2*r^2*A[r] + 
    2 lambda*r^2*dd'[r])/r^2];
(* lambda^2 A[r] - 2 A[r]/r^2 + 2 A'[r]/r + 2 lambda dd'[r] + A''[r] *)
eq2 = Expand[(r*dd''[r] + 2 dd'[r] + 3 lambda^2*r*dd[r] + 
    2 lambda*r*A'[r] + 4 lambda*A[r])/r];
(* (4 lambda A[r])/r + 3 lambda^2 dd[r] + 2 lambda A'[r] + 2 dd'[r]/r + dd''[r] *)

Combine the two expressions to obtain an expression for dd in terms of A

eqdd = -Subtract @@ First[Eliminate[{eq1 == 0, eq2 == 0, D[eq1, r] == 0}, 
    {dd'[r], dd''[r]}]]
(* -6*lambda^2*A[r]/r - 6*lambda^3*dd[r] - 3*lambda^2*A'[r] + 4*A''[r]/r + A'''[r] *)

and an expression for A alone and solve it.

eqA = Subtract @@ First[Eliminate[{D[eqdd, r] == 0, eq1 == 0}, dd'[r]]]
(* 3*lambda^4*r*A[r] - 4*A''[r]/r + 4*A'''[r] + r*A''''[r] *)

sA = DSolveValue[eqA == 0, A[r], r]
(*   (E^((-1)^(3/4) 3^(1/4) lambda r) (I + (-3)^(1/4) lambda r) C[1])/r^2 
   + (E^((-3)^(1/4) lambda r) (-I + (-1)^(3/4) 3^(1/4) lambda r) C[2])/r^2 
   + (E^(-(-1)^(3/4) 3^(1/4) lambda r) (-I + (-3)^(1/4) lambda r) C[3])/r^2
   + (E^(-(-3)^(1/4) lambda r) (I + (-1)^(3/4) 3^(1/4) lambda r) C[4])/r^2 *)

Finally, substitute A into eqdd and solve it.

sdd = Solve[Simplify[eqdd /. A -> Function[{r}, sA]] == 0, dd[r]][[1, 1]] // Values
(* -(1/(lambda r^3)) (E^((-1)^(3/4) 3^(1/4) lambda r) (I + (-3)^(1/4) lambda r) C[1] + 
    E^((-3)^(1/4) lambda r) (-I + (-1)^(3/4) 3^(1/4) lambda r) C[2] + 
    E^(-(-1)^(3/4) 3^(1/4) lambda r) (-I + (-3)^(1/4) lambda r) C[3] +
    E^(-(-3)^(1/4) lambda r) (I + (-1)^(3/4) 3^(1/4) lambda r) C[4]) *)

as requested.

Addendum: Eliminating complex numbers from the solution.

It often is desirable not to have complex numbers in the solutions of real ODEs. This can be accomplished by converting exponentials to trigonometric and hyperbolic functions, performing a complex expansion, and replacing the existing constants of integration by constants that absorb the complex numbers.

sA1 = Collect[ComplexExpand@ExpToTrig@sA, {r, Cos[_] Cosh[_], Cos[_] Sinh[_], 
    Sin[_] Cosh[_], Sin[_] Sinh[_]}, Simplify];
Cases[sA1, Cos[_] Cosh[_] z_ | Sin[_] Cosh[_] z_ | Cos[_] Sinh[_] z_ | 
    Sin[_] Sinh[_] z_ -> z, Infinity, 4];
cf = Solve[Thread[{c1, c2, c3, c4} == %], {C[1], C[2], C[3], C[4]}] // Flatten;
sA2 = Collect[sA1 /. cf, {r, Cos[_] Cosh[_], Cos[_] Sinh[_], Sin[_] Cosh[_], 
    Sin[_] Sinh[_]}, Simplify]
(* (1/(r^2))(c1 Cos[(3^(1/4) lambda r)/Sqrt[2]] Cosh[(3^(1/4) lambda r)/ Sqrt[2]] + 
             c2 Cosh[(3^(1/4) lambda r)/Sqrt[2]] Sin[(3^(1/4) lambda r)/Sqrt[2]] + 
             c3 Cos[(3^(1/4) lambda r)/Sqrt[2]] Sinh[(3^(1/4) lambda r)/ Sqrt[2]] + 
             c4 Sin[(3^(1/4) lambda r)/Sqrt[2]] Sinh[(3^(1/4) lambda r)/Sqrt[2]]) + 
   (1/r)(-((3^(1/4) (c2 + c3) lambda Cos[(3^(1/4) lambda r)/Sqrt[2]] Cosh[(
  3^(1/4) lambda r)/Sqrt[2]])/Sqrt[2]) + 
          (3^(1/4) (c1 - c4) lambda Cosh[(3^(1/4) lambda r)/Sqrt[2]] Sin[(
  3^(1/4) lambda r)/Sqrt[2]])/Sqrt[2] - 
          (3^(1/4) (c1 + c4) lambda Cos[(3^(1/4) lambda r)/Sqrt[2]] Sinh[(
  3^(1/4) lambda r)/Sqrt[2]])/Sqrt[2] - 
          (3^(1/4) (c2 - c3) lambda Sin[(3^(1/4) lambda r)/Sqrt[2]] Sinh[(
  3^(1/4) lambda r)/Sqrt[2]])/Sqrt[2]) *)

sdd2 = Collect[ComplexExpand@ExpToTrig@sdd /. cf, {r, Cos[_] Cosh[_], 
    Cos[_] Sinh[_], Sin[_] Cosh[_], Sin[_] Sinh[_]}, Simplify]
(* (1/(r^2))((3^(1/4) (c2 + c3) Cos[(3^(1/4) lambda r)/Sqrt[2]] Cosh[(
   3^(1/4) lambda r)/Sqrt[2]])/Sqrt[2] + 
            (3^(1/4) (-c1 + c4) Cosh[(3^(1/4) lambda r)/Sqrt[2]] Sin[(
   3^(1/4) lambda r)/Sqrt[2]])/Sqrt[2] + 
            (3^(1/4) (c1 + c4) Cos[(3^(1/4) lambda r)/Sqrt[2]] Sinh[(
   3^(1/4) lambda r)/Sqrt[2]])/Sqrt[2] + 
            (3^(1/4) (c2 - c3) Sin[(3^(1/4) lambda r)/Sqrt[2]] Sinh[(
   3^(1/4) lambda r)/Sqrt[2]])/Sqrt[2]) + 
   (1/(r^3))(-((c1 Cos[(3^(1/4) lambda r)/Sqrt[2]] Cosh[(3^(1/4) lambda 
   r)/Sqrt[2]])/lambda) - 
               (c2 Cosh[(3^(1/4) lambda r)/Sqrt[2]] Sin[(3^(1/4) lambda 
   r)/Sqrt[2]])/lambda - 
               (c3 Cos[(3^(1/4) lambda r)/Sqrt[2]] Sinh[(3^(1/4) lambda 
   r)/Sqrt[2]])/lambda - 
               (c4 Sin[(3^(1/4) lambda r)/Sqrt[2]] Sinh[(3^(1/4) lambda 
   r)/Sqrt[2]])/lambda) *)

For convenience, the solutions can be plotted for one non-zero c at a time.

Plot[Evaluate[Replace[sA2, {{c1 -> 1}, {c2 -> 1}, {c3 -> 1}, {c4 -> 1}}, Infinity] 
    /. {c1 -> 0, c2 -> 0, c3 -> 0, c4 -> 0, lambda -> 1}], {r, 0, 5}, 
    ImageSize -> Large, AxesLabel -> {r, A}, LabelStyle -> {15, Bold, Black}, 
    PlotLegends -> Placed[{c1, c2, c3, c4}, {.3, .8}]]

Plot[Evaluate[Replace[sdd2, {{c1 -> 1}, {c2 -> 1}, {c3 -> 1}, {c4 -> 1}}, Infinity] 
    /. {c1 -> 0, c2 -> 0, c3 -> 0, c4 -> 0, lambda -> 1}], {r, 0, 5}, 
    ImageSize -> Large, AxesLabel -> {r, dd}, LabelStyle -> {15, Bold, Black}, 
    PlotLegends -> Placed[{c1, c2, c3, c4}, {.6, .8}]]

Answer 2

ClearAll[A, r, dd, sol3, sol4, lambda]
sol3 = r^2*D[A[r], r, r] + 2*r*D[A[r], r] - 2*A[r] + 
  lambda^2*r^2*A[r] + 2*lambda*r^2*D[dd[r], r]

sol4 = r*D[dd[r], r, r] + 2*D[dd[r], r] + 3*lambda^2*r*dd[r] + 
  2*lambda*r*D[A[r], r] + 4*lambda*A[r]
lambda = 0;

s = NDSolve[{sol3 == 0, sol4 == 0, A[1] == 1, A'[1] == 1, dd[1] == 1, 
   dd'[1] == 1}, {A, dd}, {r, 1, 2}]

Plot[{A[r], dd[r]} /. s, {r, 1, 2}, PlotRange -> All]

For

lambda = 0.5;

If the integration interval is extented towards zero the singularity at r==0 get important and influences, dominates the solution:

lambda = 5.0; s2 = NDSolve[{sol3 == 0, sol4 == 0, A1 == 1, A'1 == 1, dd1 == 1, dd'1 == 1}, {A, dd}, {r, MachinePrecision, 5}] Plot[{A[r], dd[r]} /. s2, {r, 0.2, 2}, PlotRange -> All, PlotStyle -> Automatic]

Mathematica has no DSolve solution since the details are not given in the question. The boundary conditions for this solution are taken by choice. They can be selected from a wide range. Same is true for lambda.

lambda = 5.0;
s2 = DSolve[{sol3 == 0, sol4 == 0, A[1] == 0, A'[1] == .1, dd[1] == 2,
    dd'[1] == 1}, {A, dd}, {r, 0.1, 5}]
(*{{A -> Function[{r}, -(0.0333333/r^2) + 0.0333333 r], 
  dd -> Function[{r}, 3 - 1/r]}}*)

Plot[{A[r], dd[r]} /. s2, {r, 0.1, 5}, PlotRange -> All, 
 PlotStyle -> Automatic]