Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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.

share|improve this question
1  
You're sure that's all the boundary conditions you need? –  J. M. Oct 28 '12 at 15:10
    
@Nasser, w[0] == 0 and w[3/7] == 0 would look to me that OP wants to solve a boundary-value problem (which NDSolve[] has no trouble handling). Still, it looks as he doesn't have all the conditions he needs... –  J. M. Oct 29 '12 at 1:54
add comment

2 Answers

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]

Mathematica graphics

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}}
 *)
share|improve this answer
    
Er…why not add some explain for the select of the coefficients? –  xzczd Oct 29 '12 at 3:22
    
@xzczd The explanation is on the help system –  belisarius Oct 29 '12 at 3:23
    
Oh, that's true…I should have read the help more carefully. –  xzczd Oct 29 '12 at 3:28
    
@xzczd If only the help was friendlier ... –  belisarius Oct 29 '12 at 3:29
    
I can't agree more… –  xzczd Oct 29 '12 at 3:33
add comment

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} *)
share|improve this answer
add comment

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.