Problem with NDsolve for a system of equations

Question

I want to solve a system of differential equations which is not very complicated, but I cannot handle the problem with mathematica!! Please have a look at the problem and result and help me with your comments:

F=3;
ss2 = NDSolve[{3000/7*(u''[x] - 1/56*PHI''[x]*Sign[PHI'[x]]) == 0, 
               1000/7*(w''[x] + PHI'[x]) + F == 0, 
               125/2744*PHI''[x] - 1000/7*(w'[x] + PHI[x]) == 0, 
               3000/7*(u'[3/7] - 1/56*Abs[PHI'[3/7]]) == -1, 
               u[0] == 0, w[0] == 0, w[3/7] == 0, PHI'[0] == 0, PHI'[3/7] == 0}, 
         {u, w, PHI}, {x, 0, 3/7}];

And the message Mathematica generates for this:

NDSolve::ndsv: Cannot find starting value for the variable PHI'

Thank you for your attention.

Answer 1

The following gives you "reasonable" results:

DOPRIamat = {
   {1/5},
   {3/40, 9/40},
   {44/45, -56/15, 32/9},
   {19372/6561, -25360/2187, 64448/6561, -212/729},
   {9017/3168, -355/33, 46732/5247, 49/176, -5103/18656},
   {35/384, 0, 500/1113, 125/192, -2187/6784, 11/84}};
DOPRIbvec = {35/384, 0, 500/1113, 125/192, -2187/6784, 11/84, 0};
DOPRIcvec = {1/5, 3/10, 4/5, 8/9, 1, 1};
DOPRIevec = {71/57600, 0, -71/16695, 71/1920, -17253/339200, 
   22/525, -1/40};
DOPRICoefficients[5, p_] :=
  N[{DOPRIamat, DOPRIbvec, DOPRIcvec, DOPRIevec}, p];
F = 3;
ss2 = NDSolve[
  {3000/7*(u''[x] - 1/56*PHI''[x]*Sign[PHI'[x]]) == 0,
   1000/7*(w''[x] + PHI'[x]) + F == 0,
   125/2744*PHI''[x] - 1000/7*(w'[x] + PHI[x]) == 0, 
   3000/7*(u'[3/7] - 1/56*Abs[PHI'[3/7]]) == -1,
   u[0] == 0,
   w[0] == 0,
   w[3/7] == 0,
   PHI'[0] == 0,
   PHI'[3/7] == 0},
  {u, w, PHI}, {x, 0, 3/7}, 
  Method -> {"ExplicitRungeKutta", "DifferenceOrder" -> 5, 
    "Coefficients" -> DOPRICoefficients, "StiffnessTest" -> False}]

Plot[{10 u[x], 10 w[x], PHI[x]} /. ss2, {x, 0, 3/7}, PlotRange -> All, 
      AxesOrigin -> {0, 0}, Evaluated -> True]

Testing for errors in the initial conditions:

{u[0],
  w[0],
  w[3/7],
  PHI'[0],
  PHI'[3/7]} /. ss2
(*
{{-3.1407*10^-24, 
  -9.61597*10^-23, 
  -2.43547*10^-14, 
   2.29625*10^-21, 
  4 .49696*10^-13}}
 *)

Answer 2

Always adjust the parameters of complex functions, such as NDSolve, when used for more complex problems. Your problem is not so trivial. My first attempt to deal with this problem is as follows:

dEq = {3000/7*(u''[x] - 1/56*PHI''[x]*Sign[PHI'[x]]) == 0,
   1000/7*(w''[x] + PHI'[x]) + F == 0,
   125/2744*PHI''[x] - 1000/7*(w'[x] + PHI[x]) == 0
   };
bcEq = {3000/7*u'[3/7],
   u[0], w[0], w[3/7], PHI'[0], PHI'[3/7]};
bcVal = {-1, 0, 0, 0, 0, 0};

ss2Tuned = First@NDSolve[
    {dEq, Thread[bcEq == bcVal]}, {u, w, PHI}, {x, 0, 3/7},
    Method -> {"ExplicitRungeKutta",
      "DifferenceOrder" -> 5,
      "EmbeddedDifferenceOrder" -> 4,
      "StiffnessTest" -> False},
    AccuracyGoal -> 10,
    PrecisionGoal -> 10,
    MaxStepFraction -> 1/20
   ];

(You don't really need to specify RK coefficients as suggested by belisarius, it only makes the problem less readable). This code gives acceptable solution

(bcEq - bcVal) /. ss2Tuned
(* {2.90656*10^-13, -6.75703*10^-23, 3.30872*10^-24, -7.7412*10^-16, 0., 3.77476*10^-14} *)

using only 31 grid points. Since this problem is very sensible on the used method (check it!) I suggest You to reformulate it as an initial value problem, and use the shooting method to find initial values which will match the boundary values prescribed at the other end of integration interval.

solveIcProblem[a_?NumericQ, b_?NumericQ, c_?NumericQ, 
   opts : OptionsPattern[]] := Module[{dEq, ic, sol},
   dEq = {3000/7*(u''[x] - 1/56*PHI''[x]*Sign[PHI'[x]]) == 0,
     1000/7*(w''[x] + PHI'[x]) + F == 0,
     125/2744*PHI''[x] - 1000/7*(w'[x] + PHI[x]) == 0
     };
   ic = {u[0] == 0, u'[0] == a, w[0] == 0, w'[0] == b, PHI[0] == c, 
     PHI'[0] == 0};
   sol = First@NDSolve[{dEq, ic}, {u, w, PHI}, {x, 0, 3/7}, 
      FilterRules[{opts}, Options[NDSolve]]];
   {sol, {u'[3/7], w[3/7], PHI'[3/7]} - {-7/3000, 0, 0} /. sol}
  ];

rightBcDiff[a_?NumericQ, b_?NumericQ, c_?NumericQ, 
  opts : OptionsPattern[]] :=
 Last@solveIcProblem[a, b, c, opts]

FindRoot[
 rightBcDiff[a, b, c, Method -> {"ExplicitRungeKutta",
    "DifferenceOrder" -> 5,
    "EmbeddedDifferenceOrder" -> 4,
    "StiffnessTest" -> False},
  AccuracyGoal -> Infinity,
  PrecisionGoal -> 14, 
  StartingStepSize -> 10^-3], {{a, -0.0023}, {b, 0.22}, {c, -0.216}},
  AccuracyGoal -> Infinity, PrecisionGoal -> 14
 ]
solveIcProblem[a, b, c] /. % // Last

This method results with

(* {a -> -0.00233333, b -> 0.2205, c -> -0.216} *)
(* {-4.11997*10^-17, -7.77433*10^-13, 1.01738*10^-13} *)

Compared with previous approach

solveIcProblem @@ {u'[0], w'[0], PHI[0]} /. ss2Tuned // Last
(* {2.54158*10^-8, 2.12706*10^-8, -1.5725*10^-9} *)