# NDSolve solution does not satisfy boundary conditions

I am trying to solve a 1st order non-linear ODE

W[y]*W'[y] + W[y]*v + Fnum == 0 /. v -> 10

Fnum = 0.05 - 1.66667 y + (270.27 y)/(1 + 270.27 y) - (660. y)/(1 + 1666.67 y)


The function Fnum has three zeros. two of them are at y=StartY and y=EndY.

In:= StartY = 0.00021578513459560855;

In:= EndY = 0.3870868735615052;

In:= Fnum = 0.05 - 1.66667 y + (270.27 y)/(1 + 270.27 y) - (660. y)/(1 +
1666.67 y)

Out= 0.05 - 1.66667 y + (270.27 y)/(1 + 270.27 y) - (660. y)/(1 + 1666.67 y)

In:= Fnum /. y -> StartY

Out= 2.62289*10^-9

In:= Fnum /. y -> EndY

Out= -5.10119*10^-7


My boundary condition is W[y=StartY] == 0. However, if it set this NDSolve gives an error

s = NDSolve[{W[y]*W'[y] + W[y]*v + Fnum == 0 /. v -> 10,
W[StartY] == 0 /. v -> 10}, W, {y, StartY, EndY}];


Power::infy: Infinite expression 1/0. encountered. >>

NDSolve::ndnum: Encountered non-numerical value for a derivative at y == 0.00021578513459560855. >>

I am guessing this happens because Mathematica cannot solve W[y]' at y=StratY

W'[y] = -v + Fum[y]/W[y] ~ -v + 0/0.


So, I am now trying to find W'[y=StartY] analytically by expanding Fnum near y=StartY and I find that W'[y=StartY]=dW. I show here only that my solution satisfies the differential equation with Fnum[StartY]=0 and boundary condition W[StartY=0]:

In:= Slope = Coefficient[Normal[Chop[Series[Fnum, {y, StartY, 1}]]], y, 1]

Out= -117.385

In:= dW  = (Sqrt[v^2 - 4*Slope] - v)/2;

In:= Solve[W*dW + W*v + Fnum == 0 /. {y -> StartY, v -> 10}, W]

Out= {{W -> -1.54903*10^-10}}


Now I go back to NDSolve and I am trying W'[StartY] = dW as boundary condition

s = NDSolve[{W[y]*W'[y] + W[y]*v + Fnum == 0 /. v -> 10, W'[StartY] == (Sqrt[v^2 - 4*Slope] - v)/2 /. v -> 10}, W, {y, StartY, EndY}]


I get the following error:

NDSolve::icordinit: The initial values for all the dependent variables are not explicitly specified. NDSolve will attempt to find consistent initial conditions for all the variables. >>

NDSolve::icres: NDSolve has computed initial values that give a zero residual for the differential system, but some components are different from those specified. If you need them to be satisfied, giving explicit initial values to all dependent variables is recommended. >>

NDSolve::ndsz: At y == 0.0854996004928462, step size is effectively zero; singularity or stiff system suspected. >>

In spite of the errors, NDsolve does manage to find a solution. However, NDSolve does find a solution which with W[StartY]=1, which is not consistent.

In:= Evaluate[W[StartY] /. s]

Out= {1.}


What am I doing wrong? What do the error mean? Why NDSolve finds a solution which is not consistent with my boundary conditions?

• Welcome to Mathematica.SE! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Read the faq! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! – user9660 Nov 24 '14 at 16:59
• I formatted your question for you this time. Please learn how to do it for your nxt question. Thanks. – Dr. belisarius Nov 24 '14 at 17:05
• The styling of your question can be improved. Please take the time to read the markdown help page. You'll better communicate your problem when you use the right formatting. – user9660 Nov 24 '14 at 17:08
• You use Slope without ever assigning it a constant numeric value. Think of trying to do a simple manual Euler's approximation of f'[x]=Slope when I won't tell you what the value of Slope is. But even when I assign some constant value to Slope, there still appear to be warning messages and those almost always must be understood and fixed before any sensible result can be had – Bill Nov 24 '14 at 17:12
• @user22375 (You need to put a @ in front of names for the person to get pinged.) According to your own code Fnum is not equal to zero: In:= Fnum /. y -> StartY, Out= 2.62289*10^-9. Execute 2.62289*10^-9 == 0 to check. (Another side point: It will encourage more users to try out your code if they can copy and paste it in Mathematica without editing it. The In and Out tags make it inconvenient. A convention many follow is to put output in comments.) – Michael E2 Nov 25 '14 at 13:13

The problem with NDSolve methods is that (as far as I know) it always check boundary points. And this a problem, as you noted. I would recommend you to look at this wonderful approach. I adapted it to your problem below

v = 10;
eq = (W[y]*W'[y] +
W[y]*v + (0.05 -
1.66667 y + (270.27 y)/(1 + 270.27 y) - (660. y)/(1 +
1666.67 y)));
Fnum[y_] :=
0.05 - 1.66667 y + (270.27 y)/(1 + 270.27 y) - (660. y)/(1 +
1666.67 y);
startY = 0.00021578513459560855;
endY = 0.3870868735615052;
pts = 500;
ygrid = Range[startY, endY, 1/(pts - 1)];
dwdy = NDSolveFiniteDifferenceDerivative[1, ygrid,
"DifferenceOrder" -> 8]["DifferentiationMatrix"] ;
varsW = Table[Unique[w\$], {Length[ygrid]}];
eqns = varsW*(dwdy.varsW) + varsW*v + Fnum[ygrid];
eqns[] = varsW[]; (* boundary condition W[StartY]\[Equal]0 *)

sol = varsW /. FindRoot[eqns, Thread[{varsW, 1/2}]];
Wsol = Interpolation[Transpose[{ygrid, sol}],
InterpolationOrder -> 2];


And here are plots. The solution

Plot[Wsol[y], {y, startY, endY}, PlotRange -> All] and the error

Plot[eq /. W -> Wsol, {y, startY, endY}, PlotRange -> All] I noticed that increasing DifferenceOrder gives better results. You can also enlarge the grid itself.

• this is a very elegant approach. Because Fnum varies rapidly near startY, W' is large there, and your error estimate also is large there, solving the equation on a very fine mesh, say 10.^-5, over the region {startY,0.02} may have merit. – bbgodfrey Nov 26 '14 at 14:20

In an earlier comment, I observed that NDSolve found the equation to be stiff, even when the calculation was begun a short distance d from StartY. Because the third root of Fnum is 0.000797339, d must be quite small, say, 0.0001. However, even with this choice of d, NDSolve promptly stops because the equation is stiff. It turns out, however, that NDSolve can solve the equation with the choice v -> 11.3 (or larger) with the corresponding negative slope value, -17.8691. (The corresponding positive slope, 6.56913, leaves the equation stiff.) The solution is Solutions for smaller values of v can be obtained by a change of variables, W -> -Sqrt[q]. (Choosing the positive root is equivalent to choosing the negative value of slope, which leaves the equation stiff.) As one would expect, the resulting solution for v -> 11.3 is the square of the curve above.

For v -> 10.0, q is complex: The negative Sqrt of this curve is the desired result.

Origin of complex-valued solution

A blow-up of the figure just above at small y helps to explain why some solutions are complex-valued. The solution vanishes at approximately y=.0052497, just as the solution does at StartY. However, here, Fnum does not vanish too. Consequently, the leading terms in q in the vicinity of this zero are 0.544893(.0052497-y) + 14.7634(.0052497-y)^(3/2), as can be demonstrated by direct substitution. Thus, although the solution is continuously differentiable at y=.0052497, it has a branch point there, giving rise to the abrupt onset of complex values. Incidentally, it is this zero, and not the one at StartY, that prevents a straightforward solution for W using NDSolve.

v=0 analytical solution

The equation can be solved analytically for v->0, yielding the following expression for q=W^2.

0.037919170215639376 + y*(-1.3080015839968322 + 1.66667*y) - 0.0004751980992057023*Log[0.0005999988000024 + y] +
0.007400007400007398*Log[0.0037000037000037003 + y]


This analytical solution may be useful in evaluating numerical solution methods for the equation. As shown in the plot below, it is negative almost everywhere, indicating that W is complex. (q is positive only in the region {StartY, 0.0015}`.) 