Re: Mathematica takes wrong derivative
- To: mathgroup at smc.vnet.net
- Subject: [mg59703] Re: Mathematica takes wrong derivative
- From: Maxim <ab_def at prontomail.com>
- Date: Thu, 18 Aug 2005 00:16:35 -0400 (EDT)
- References: <dds971$8vd$1@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
On Tue, 16 Aug 2005 08:47:29 +0000 (UTC), Mukhtar Bekkali
<mbekkali at gmail.com> wrote:
> I believe Mathematica 5.0.1.0 takes wrong derivative. See my code below
> where I plotted the function D2 and its derivative. The derivative is
> obviously wrong.
>
Consider a simpler example:
y = x;
f[x_?NumericQ] := x*y
f[1] will substitute 1 for x on the right-hand side, giving 1*y, which in
turn will evaluate to x. Inside of Plot we're effectively evaluating
Block[{x = 1}, f[x]], so this time the evaluation chain will look like
Block[{x = 1}, f[x]] --> Block[{x = 1}, f[1]] --> Block[{x = 1}, x] --> 1
(Block won't be reevaluated on each step, I just want to stress the fact
that all computations are performed in a context where the value of x is
1). So this kind of definition will work only inside Block-like
constructs, such as Plot or numerical root finding/optimization functions.
Furthermore, even there it may conflict with the setting Compiled -> True.
Next, what happens when we try to find the derivative of f numerically? In
the simplest case we're using the approximation (f[x + eps] - f[x])/eps,
which becomes ((x + eps)*y - x*y)/eps --> ((x + eps)*x - x*x)/eps --> x.
So the 'derivative' of f[x] will be equal to x, not 2*x. Indeed,
Plot[f'[x], {x, -1, 1}] will draw the line y = x.
So it is a good idea to make the dependencies on \[Delta] (and probably on
\[Beta] too) explicit:
{r, u, \[Lambda], \[Phi], \[Gamma][\[Beta]_]} = {2, 4, (9/10)*(u/(1 + u)),
1 - (\[Lambda]/(1 - \[Lambda]))*(1/u), (1 - \[Beta])^r/(\[Beta]^\[Beta]*(1
- \[Beta])^(1 - \[Beta]))};
z1[\[Beta]_] = (\[Beta]^\[Beta]*(1 - \[Beta])^(1 -
\[Beta])*\[Gamma][\[Beta]]*(u/(1 + u)))^(-1);
z2[\[Beta]_, \[Delta]_] = (\[Delta]^\[Beta]*(1 - \[Delta])^(1 -
\[Beta])*\[Gamma][\[Beta]]*((u*(1 - \[Beta]))/(u*(1 - \[Beta]) + 1))^(1 -
\[Beta]))^(-1);
V21[\[Delta]_] := NIntegrate[z2[\[Beta], \[Delta]]^(-u), {\[Beta], 0,
(\[Delta]*(1 + u))/u}]
V22[\[Delta]_] := NIntegrate[z1[\[Beta]]^(-u), {\[Beta], (\[Delta]*(1
+ u))/u, 1}]
D2[(\[Delta]_)?NumericQ] := (\[Delta]*V21[\[Delta]]
+ NIntegrate[(\[Beta]*u)/(1 + u)/z1[\[Beta]]^u, {\[Beta], (\[Delta]*(1
+ u))/u, 1}])/(V21[\[Delta]] + V22[\[Delta]])^\[Phi]
Plot[D2[\[Delta]] // Evaluate, {\[Delta], 0, u/(1 + u)}];
Plot[D2'[\[Delta]] // Re // Evaluate, {\[Delta], 0, u/(1 + u)}];
Now the derivative will be evaluated correctly, but not very precisely --
we had to chop off the imaginary part and there are some irregularities in
the left part of the plot. It is better to allow D2 to be evaluated for
symbolic arguments as well, because D has rules for symbolic
differentiation of NIntegrate:
Clear[D2]
D2[\[Delta]_] := (\[Delta]*V21[\[Delta]] + NIntegrate[(\[Beta]*u)/(1
+ u)/z1[\[Beta]]^u, {\[Beta], (\[Delta]*(1 + u))/u, 1}])/(V21[\[Delta]]
+ V22[\[Delta]])^\[Phi]
Plot[D2[\[Delta]] // Evaluate, {\[Delta], 0, u/(1 + u)}];
Plot[D2'[\[Delta]] // Evaluate, {\[Delta], 0, u/(1 + u)}];
This will produce a smooth curve and will actually take less time, not
only because numerical differentiation won't be used but also because
Mathematica won't have to subdivide the problematic intervals many times
(controlled by the PlotDivision setting) trying to smooth the
irregularities.
Maxim Rytin
m.r at inbox.ru