Re: numerical derivatives at very high x-values: precision problem
- To: mathgroup at smc.vnet.net
- Subject: [mg64043] Re: [mg64021] numerical derivatives at very high x-values: precision problem
- From: "Carl K. Woll" <carlw at wolfram.com>
- Date: Sat, 28 Jan 2006 02:24:18 -0500 (EST)
- References: <200601271838.k0RIcAaA014635@mx1.wolfram.com>
- Sender: owner-wri-mathgroup at wolfram.com
Nadir Kaplan wrote:
[snip]
>>>The values are f.e. (from a set of recursion relations) :
>>>
>>>------------------------------------------------------------------------
>>
>>----
>>
>>>J[40]= 1.12970494019605678505771742460*10^12
>>>
>>>mu[40]= 3.38911605114868836337349826143*10^12
>>>
>>>V[40]= -8.47278910237611083984375000000*10^11
>>>
>>>T1[40]= 5.64852470100879638671875000000*10^11
>>>
>>>
>>>
>>>The function to be partial-differentiated:
>>>
>>>------------------------------------------------------------
>>>
>>>lngamma1[J_, mu_, V_, t1_]:=
>>>
>>> Module[{mu2, U1, U2, Umax}, mu2=SetPrecision[mu/(2*d), 30];
>>>
>>> U1=SetPrecision[0, 30];
>>>
>>> U2=SetPrecision[1.5*mu2+func1[0.5*mu2, t1], 30];
>>>
>>> Umax= Max[U1, U2];
>>>
>>>
>>>SetPrecision[Umax+Log[Exp[U1-Umax]+2*Exp[U2-Umax]*func2[0.5*mu2, t1]],
>>
>>30]
>> ^^^^^^^^^^^^^^^^
>> ||||||||||||||||
>>
>>why not just U1-Umax?
>
> lngamma1 is one of the five recursion relations, and I established a general
> style to easily control them. I am able to see the exponents immediately by
> writing so.
>
Well, Umax + Log[Exp[U1-Umax] equals U1, so including this complication
in your function is odd. It's also odd that your function has arguments
V and J which don't appear in the function definition, and the function
definition has the symbol d which is undefined (perhaps it is 40?).
>>> ]
[snip]
>>
>>First, eliminate all of the SetPrecision statements, they aren't
>>helpful. Next, you don't tell us what func1 and func2 are. The rest of
>>your lngamma1 function can be symbolically differentiated by Mathematica
>>without problems, so unless func1 and func2 are undifferentiable
>>functions, there is no need to use numerical derivatives.
>>
>>If func1 or func2 is symbolically undifferentiable, then you should
>>consider using the ND standard package to do differentiation. It uses
>>Richardson extrapolation to find derivatives, which is much more
>>reliable than the derivative formula you give above.
>>
>
> Sorry, I forgot to paste them:
>
> func1[A_, t1_] := Sqrt[2*t1^2 + A^2];
>
> func2[A_, t1_] := A/Sqrt[2*t1^2 + A^2]*0.5*(1 -
> Exp[-2*Sqrt[2*t1^2 + A^2]]) + 0.5*(1 + Exp[-2*Sqrt[2*t1^2 + A^2]]);
>
These functions are differentiable, so there is no reason to worry about
numerical derivatives.
> (SetPrecision statements eliminated as you suggested above.)
>
Perhaps the presence of the machine numbers 0.5 and 1.5 led you to mess
with SetPrecision. Instead, just replace them with 1/2 and 3/2.
> The problem is, when symbollically differentiating the recursion relations
> without finding and factorizing the largest exponent Umax, it definitely
> leads to an overflow error ( I tried it in another system, not in Mathematica but
> the problem is already at the exponents, their components grow too much too
> early. This is why I should keep the log of the recursion relations, on the
> other hand, d/dx (log f(x))=(df/dx)/f -->overflow because of the large
> exponents )
>
Assuming d is 40, I get underflow errors because of things like
In[20]:=
Log[Exp[-10^12]+2]//N
During evaluation of In[20]:=
General::"unfl": ""Underflow occurred in computation."
Out[20]=
0.693147
I don't know how to avoid these errors, but they seem benign.
When I tried your example, Mathematica had no problems taking
derivatives of lngamma1 and evaluating it at numerical values, other
than issuing Underflow warnings.
Carl Woll
Wolfram Research
>
>>Finally, instead of writing lngamma1 as
>>
>>lngamma1[J_, mu_, V_, t1_]:=...
>>
>>you should consider writing it as
>>
>>lngamma1[J_, mu_, V_][t1_]:=...
>
> I'm new at Mathematica; what is the difference?
>
> Nadir Kaplan
>