Mathematica Stack Exchange is a question and answer site for users of Mathematica. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

Can Mathematica solve $\; d^2/dx^2{}+ (1/x) d/dx\;$ form of equations?
For example, I am trying to solve the system of equations given below:

eqn1 = D[w[x], {x, 2}] + (1/x) (D[w[x], x]) - w[x];
eqn2 = D[v[x], {x, 2}] + (1/x) (D[v[x], x]) - (v[x] - u[x]);
eqn3 = D[u[x], {x, 2}] + (1/x) (D[u[x], x]) + v[x] - u[x] - w[x];
eqnSet = {eqn1 == 0, eqn2 == 0, eqn3 == 0, w[1] == 0, v[1] == 0, 
      u[1] == 0, w'[0] == 0, v'[0] == 0, u'[0] == 0};
DSolve[eqnSet, {w[x], v[x], u[x]}, x]

I am getting DSolve returning the input statement unevaluated.

share|improve this question
    
Is this 1/x term OK if you specify conditions for x=0? – Sjoerd C. de Vries Jun 25 '13 at 21:56
1  
This DSolve[{x D[w[x], {x, 2}] + (D[w[x], x]) - x w[x] == 0, w[1] == 0, w'[0] == 0}, w[x], x] works but due to the initial conditions it's a trivial solution. – b.gatessucks Jun 25 '13 at 21:57
    
I missed to see that. Thanks to both of you for pointing out my negligence. Still do we need NDSolve[] or can DSolve[] do the job? – user8250 Jun 25 '13 at 23:13
    
Just try and see. – b.gatessucks Jun 26 '13 at 5:01
    
@OleksandrR I disagree that the question arises due to a simple mistake. In fact, it arises, because DSolve is unable to solve a system of equations that you or I could solve by hand without much difficulty. – bbgodfrey Sep 27 '15 at 3:22

There is nothing wrong with the 1/x term in the ODEs. In fact, eqn1 is the standard equation for the Modified Bessel Function of order 0. However, because the three ODEs and their boundary conditions as specified in the question all are homogeneous, the solution must be identically zero. Because this result is not particularly interesting, consider the solution without boundary conditions.

DSolve[{eqn1 == 0, eqn2 == 0, eqn3 == 0}, {w[x], v[x], u[x]}, x]

DSolve also claims to be unable to solve these equations, even though they are solvable. So, we must give it some help.

DSolve[eqn1 == 0, w[x], x][[1, 1]]
(* w[x] -> BesselJ[0, I x] C[1] + BesselY[0, -I x] C[2] *)

It is very strange that DSolve does not express w in terms of Modified Bessel Functions. However, with additional help, the answer can be simplified.

cf[e_] := LeafCount[e] + 100 Count[e, _BesselY, {0, Infinity}]
s1 = FullSimplify[%%, x > 0, ComplexityFunction -> cf] /. 
    C[1] -> C[1] + I C[2] /. C[2] -> -C[2] π/2    
(* w[x] -> BesselI[0, x] C[1] + BesselK[0, x] C[2] *)

Next, we observe that eqn2 + eqn3 can be simplified greatly by substituting u[x] + v[x] -> z[x].

Simplify[(eqn2 + eqn3) /. {v'[x] -> -u'[x] + z'[x], v''[x] -> -u''[x] + z''[x]} /. s1];
s2 = FullSimplify[DSolve[% == 0, z[x], x][[1, 1]], x > 0]
(* z[x] -> (-1 + BesselI[0, x]) C[1] + BesselK[0, x] C[2] + C[4] + C[3] Log[x] *)

To obtain an expression for v[x], we eliminate u[x] in terms of z[x] from eqn2 and solve.

eqn2 /. {u[x] -> -v[x] + z[x]} /. s2
s3 = FullSimplify[DSolve[% == 0, v[x], x][[1, 1]], x > 0] /. 
    C[5] -> C[5] + I C[6] /. C[6] -> -C[6] π/2
(* v[x] -> 1/2 (-C[1] + 2 BesselI[0, x] C[1] + 2 BesselK[0, x] C[2] + C[4] + 
   2 BesselI[0, Sqrt[2] x] C[5] + 2 BesselK[0, Sqrt[2] x] C[6] + C[3] Log[x]) *)

Finally, we obtain u[x] from z[x].

Simplify[((s2 /. z[x] -> u[x] + v[x]) /. Rule -> List) - (s3 /. Rule -> List)] 
    /. List -> Rule
(* u[x] -> 1/2 (-C[1] + C[4] - 2 BesselI[0, Sqrt[2] x] C[5] - 
   2 BesselK[0, Sqrt[2] x] C[6] + C[3] Log[x]) *)

Back-substitution shows that these three expression indeed satisfy the oriinal equations. Thus, DSolve can solve this problem, but only with considerable assistance.

Alternative, more compact solution

Because the differential operators are the same in all three equations, the equations can be diagonalized, solved with DSolve, and transformed back.

var = {w[x], v[x], u[x]};
{b, m} = CoefficientArrays[{eqn1, eqn2, eqn3}, var] // Normal;
{val, vec} = Eigensystem[m]
(* {{-2, -1, 0}, {{0, -1, 1}, {1, 1, 0}, {0, 1, 1}}} *)

DSolve[Thread[(b + val var) == 0], var, x] // Flatten;
FullSimplify[%, x > 0, ComplexityFunction -> cf] /. 
  {C[1] -> C[1] + I C[2], C[3] -> C[3] + I C[4]} /. {C[2] -> -C[2] π/2, C[4] -> -C[4] π/2};
Transpose[vec].(% /. Rule[_, z_] -> z);
Thread[Rule[var, %]]
(* {w[x] -> BesselI[0, x] C[3] + BesselK[0, x] C[4], 
    v[x] -> -BesselI[0, Sqrt[2] x] C[1] - BesselK[0, Sqrt[2] x] C[2] + 
                BesselI[0, x] C[3] + BesselK[0, x] C[4] + C[6] + C[5] Log[x], 
    u[x] -> BesselI[0, Sqrt[2] x] C[1] + BesselK[0, Sqrt[2] x] C[2] + C[6] + C[5] Log[x]}*)

DSolve should, in my opinion, be able to perform the diagonalization and provide solutions in terms of Modified Bessel Functions without assistance.

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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