# why I can not find the correct derivative with D

the coupled ODE is

f[a_?NumericQ] := NDSolve[{(x[r] (192 + (20 r^3 y''[r])/y[r]))/
r^4 == 5 x''[r]/r,
1000 (2 (x'[r])^2 + 2 x[r] x''[r]) +
250 (2 (y'[r])^2 + 2 y[r] y''[r]) == -5/(4 r),
x == 0,
y == 1/10 Sqrt[1/2 Log[5/2]],
y' == -((1 - Log[5/2])/(20 Sqrt[2 Log[5/2]])),
x' == a,
WhenEvent[y[r] < 10^(-6), end = r; "StopIntegration",
"DetectionMethod" -> "Interpolation"]},
{x, y}, {r, 1, 5/2},
Method -> {"ExplicitRungeKutta", "StiffnessTest" -> False}];


I want to plot this expression:

g[a_?NumericQ] := NIntegrate[Evaluate[{(Pi (256 x[r]^2 + 16 y[r]^2))/(2 r^3) +
5/2 Pi ((x'[r])^2 + (y'[r])^2)} /.
f[a]], {r, 1, end}][[1, 1]]


the plot of g[a] with respect to a in the region {a, 0.0026, 0.008} there is a minimum there, so I use D the find the derivative at a=0.0028

D[g[a], a] /. a -> 0.0028
0.0500445


the result is rather surprising, since the figure show the slope at a=0.0028 should be negative. what's wrong with D?

thanks a lot

one more fact, findminimum and NMinimize give me the correct answer of a where the slope is 0

• Have you looked at all the error/warning messages being printed? !Mathematica graphics it is best to take care of all these first before. – Nasser Jun 2 '14 at 2:21
• @Nasser, I set workingprecision by default, the problem still exsits. – 3c. Jun 2 '14 at 2:39
• If you are looking for a minimum, why are you plugging in a=0.0028? That minimum seems to occur past 0.004. I've punched in the following: Manipulate[ Plot[ Evaluate[ReplaceAll[{x[r], y[r]}, f[a]]], {r, 0.9, 2.6}, PlotRange -> {-0.02, 0.07}], {a, 0.0026, 0.008}] I also changed the NDSolve bounds to give me some wiggle room with the r values (notice I'm plotting just outside the range {1,2.5}. This seems to indicate that solutions go haywire at or right around 2.5 when a<0.333. Finding g[a] for those values (and D[g[a],a]) will be dodgy. – Kellen Myers Jun 2 '14 at 5:44
• See also: Manipulate[ ListPlot[ Table[{a, g[a]}, {a, 0.0026, 0.004, da}], Joined -> True, PlotRange -> {0.1, 0.3} ], {da, {0.0005, 0.0002, 0.0001, 0.00005, 0.00002, 0.00001}} ] – Kellen Myers Jun 2 '14 at 5:51
• You're forcing the system to approximate the derivative numerically. Derivative can indeed do this, but I've found it to be not so reliable in some cases. Try the ND function instead from the NumericalCalculus package. It gives -15.1987. I do not know how accurate this is and how much it's influenced by the precision used in the implementation of g. I think the precision of that NIntegrate and NDSolve may have a significant effect over this. – Szabolcs Jun 2 '14 at 14:03

This answer will be a combination of some of the comments and then a few methods to try to find the value of $a$ you seem to want.

First of all, as I note, $a=0.0028$ is clearly not the minimum (nor near it). Now, that's not a problem, you should still be able to insist that the derivative $g'(0.0028)$ be computed correctly. You are fairly sure that this is wrong:

D[g[a], a] /. a -> 0.0028
out: 0.0500449


As Szabolcs points out, you should be using a numerical derivative:

Needs["NumericalCalculus"]
ND[g[a], a, 0.0028]
out: -15.1987


Now, it was your concern that $g'(0.0028)$ should be negative, however it's interesting that this is convincing enough when we can check and see:

a0 = 0.0028;
g0 = g[a0];
m = ND[g[a], a, a0];
Show[
ListPlot[
ParallelTable[{a, g[a]}, {a, 0.0026, 0.008, 0.00001}],
Joined -> True
],
Plot[m (a - a0) + g0, {a, a0 - 0.001, a0 + 0.001}]
] The tangent line actually doesn't match -- it seems maybe to match the overall slope (a sort of average slope in a nearby interval of $a$ values, perhaps), but it's still not right.

The distraction here, which shouldn't be a problem, is the issue of some singularity or stiffness (which, I should note, you've set "StiffnessTest" -> False) at $a=0.0033$ or thereabouts. I mentioned above that you can see a "kink" in $g(a)$ just by plotting it:  And I noted that if you widen the NDSolve interval a bit and do some plotting, you can see that $x(t)$ and $y(t)$, with the substitution $f(a)$, have a blowup at $t=5/2$ when $a<0.0033$.

Manipulate[
Plot[
Evaluate[ReplaceAll[{x[r], y[r]}, f[a]]],
{r, 0.9, 2.6}, PlotRange -> {-0.02, 0.07}],
{a, 0.0026, 0.008}]  Now, the ND command has already had issues above, so it's not advisable to use it, but we could try:

FindRoot[ND[g[a], a, b], {b, 0.004}]
output: {b -> -0.00477551}


That's not right, and it's not even close. I assume because the minimum is too close to this bad value of $a$ (and with a steep region between the two values of $a$). The FindRoot command won't be very useful here, I think.

You can try NMinimize but I'm not sure that's the right command -- it's not the type of command I'd normally use in this instance because it's meant to minimize a function over its entire domain. (Worse, this is an interpolated function, so beyond the interpolation bounds, who knows what happens.) You can try it constrained if you want:

NMinimize[{ND[g[a], a, b], 0.34 <= b <= 0.5}, b]


I couldn't even get this one to work. Maybe there's a better way to use NMinimize that I'm not seeing, but I $Aborted this one. Here's what you should do: FindMinimum[g[a], {a, 0.34, 0.5}] out: {0.0841652, {a -> 0.00417839}}  This value of$a$is confirmed visually by the plot of$g(a)$you've already got. You can test this out: ND[g[a], a, 0.00417] out: 2.21822  D'oh. That's wrong. Let's try to get some precision: ND[g[a], a, 0.00417, Terms -> 30] out: -0.0396302  We can even do a quick list to see: ParallelTable[ND[g[a], a, a0], {a0, 0.0041, 0.0042, 0.00001}] ParallelTable[ND[g[a], a, a0, Terms -> 30], {a0, 0.0041, 0.0042, 0.00001}]  I've omitted output, but the first list looks useless (values like 2.20574) while the second list works properly and gives relatively accurate numbers. The smallest values of ND[g[a],a,a0] occur near the value of$a$we already found. This also let's us try to reprise our attempt to use FindRoot as follows: FindRoot[ND[g[a], a, b, Terms -> 30], {b, 0.00417}]  But even so, this takes more than 20 minutes to compute, which is longer than I'm willing to wait. There is a chance that using about 30 terms may give you accurate enough NDs that the FindRoot will work (but it might not). You can use more Terms to examine various values of the derivative. You might notice something odd happen: ParallelTable[ND[g[a], a, a0, Terms -> 30], {a0, 0.002, 0.008, 0.0001}] out: {-46.8032, -43.3339, -38.4327, -18.6678, 5.38868, -11.5114, -33.1248, 30.2456, -83.6159, -47.3033, -0.866833, -5.92381, -87.673, -558.373, 20.4022, -29.6248, -11.608, -6.2123, -3.61962, -2.10121, -1.10858, -0.412773, 0.0988319, 0.48823, 0.792357, 1.03462, 1.23061, 1.39109, 1.52376, 1.63425, 1.72682, 1.80469, 1.8704, 1.92594, 1.97291, 2.0126, 2.04607, 2.07419, 2.09769, 2.11719, 2.1332, 2.14615, 2.15644, 2.16437, 2.17023, 2.17425, 2.17666, 2.17763, 2.17732, 2.17588, 2.17343, 2.17008, 2.16594, 2.16108, 2.15559, 2.14954, 2.14299, 2.13598, 2.12859, 2.12084, 2.11278}  Notice that the first 10-15 values here are erratic. More Terms might yield some idea of whether these are "correct" in some sense (for example, fixing that tangent line we started with), but I think all this should really indicate that anything we do with values prior to$a=0.0033$is probably incorrect or at least quite strange and at least dubious. (Note: Don't use too many Terms or all the derivatives will blow up or vanish, I think due to issues with using too many terms in the Euler sums used by ND that go beyond the precision in the NDSolve and/or NIntegrate in the definitions of$f$and$g$. There may be other options in ND that can be used to get more accurate output.) But I hope this clears up the issues. In the end, the lesson seems to be: Use FindMinimum` (and a bonus lesson: Steer clear of bad values of$a\$).

• thanks very much. in the end I think all such strange behaviors originate from the stiffness problem of the ODE. now I force the function to end before the stiffness comes out, then everything is well behaved – 3c. Jun 3 '14 at 18:49