Re: numerical derivatives at very high x-values: precision problem
- To: mathgroup at smc.vnet.net
- Subject: [mg64044] Re: [mg64021] numerical derivatives at very high x-values: precision problem
- From: "Nadir Kaplan" <nadir.kaplan at pclabs.gen.tr>
- Date: Sat, 28 Jan 2006 02:24:23 -0500 (EST)
- Sender: owner-wri-mathgroup at wolfram.com
Hi,
Firstly I want to thank you for your attention.
At the beginning I wrote the question by thinking that my problem is
qualitative, so I just left out some information by just focusing on the
precision problem.
In fact, this is a physics problem, and d=3 (dimension).
There are four parameters t, J, V and \mu, lngamma1 ( \equiv: log(\gamma{1})
) is the only one of the five recursion relations which doesn't depend on V
and J, but nevertheless I include them and those U1's and U2's even if
they're 0, in order to compare and contrast with the other recursion
relations easily, whether there's a typing mistake or not.
This recursion relations serve to several purposes, for months they have
worked properly for the other ones, but for the last duty I should calculate
the derivatives which are unstable for great numbers which I have included
below before.
> -----Original Message-----
> From: Carl K. Woll [mailto:carlw at wolfram.com]
To: mathgroup at smc.vnet.net
> Subject: [mg64044] Re: [mg64021] numerical derivatives at very high x-values:
> precision problem
>
> 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
> >