# Problem with NDsolve for a system of equations

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'


-
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

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}}
*)

-
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 – Dr. 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 ... – Dr. belisarius Oct 29 '12 at 3:29
I can't agree more… – xzczd Oct 29 '12 at 3:33

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} *)

-