Re: numerical derivatives at very high x-values: precision problem
- To: mathgroup at smc.vnet.net
- Subject: [mg64042] 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:15 -0500 (EST)
- Sender: owner-wri-mathgroup at wolfram.com
Hi, > -----Original Message----- > From: Carl K. Woll [mailto:carlw at wolfram.com] To: mathgroup at smc.vnet.net > Subject: [mg64042] Re: [mg64021] numerical derivatives at very high x-values: > precision problem > > Nadir Kaplan wrote: > > Hi, > > > > When calculating partial numerical derivatives, namely the formula > [f(x+h, > > y)-f(x-h, y) ] / (2h), at very high x values, let's say x~10^12 , I > should > > choose h=(10^a)*|x| where x should be <-16 in order to get accurate > results. > > We know that the MachineEpsilon is approx. 2.2*10^(-16) , but I think I > > should somehow use a number less than that for h, otherwise the answers > are > > not satisfiying. How will I do this in order to get high precision (or > > accuracy) ? I added SetPrecision command everywhere, but it seems it > didn't > > work out well, besides it looks dumb (or not? I'm new at Mathematica ). > > > > > > > > Any help will be appreciated... > > > > > > > > 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. > > > > > ] > > > > > > > > The derivative formula: > > > > ---------------------------------------------------------------- > > > > derivt1[func_, J_, mu_, V_, t1_]:= > > > > Module[ {h}, h=SetPrecision[10^-15*Abs[t1],30] ; > > > > SetPrecision[(func[J, mu, V, t1+h]-func[J, mu, V, t1-h])/(2*h), 30] > > > > ] > > > > > > 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]]); (SetPrecision statements eliminated as you suggested above.) 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 ) > 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