Re: NDSolve::ndsz question
- To: mathgroup at smc.vnet.net
- Subject: [mg63840] Re: NDSolve::ndsz question
- From: Maxim <m.r at inbox.ru>
- Date: Wed, 18 Jan 2006 03:03:03 -0500 (EST)
- References: <dpnqoa$6rh$1@smc.vnet.net> <dpvmjv$ltf$1@smc.vnet.net> <dq53pi$aak$1@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
On Thu, 12 Jan 2006 08:25:22 +0000 (UTC), Peter Pein <petsie at dordos.net>
wrote:
>
> Regarding the use of finite/infinite precision using NDSolve,NSum etc.,
> I do
> not fully understand, why the effect is so immense. I encountered the
> following some time ago:
>
> s[k_, x_] := ((2*k - x)/(2*k + x))^k - Exp[-x]
> Series[s[k, x], {k, Infinity, 3}]//Normal
> --> -(x^3/(E^x*(12*k^2)))
>
> Therefore Sum[s[k,x],{k,0,Infinity}] should converge at least for x>0
> (due to
> s[k,x] there are singularities at nonpositive even integervalues of x).
>
> SetOptions[NIntegrate, SingularityDepth -> 50, MaxRecursion -> 50];
>
> f[x_?NumericQ] := NSum[
> s[k, x], {k, 0, }, WorkingPrecision -> $MachinePrecision]
>
> without setting WorkingPrecision, the evaluation of f[] fails (*why?*).
>
> f[0.25]
> --> ComplexInfinity
>
> f[1/4]
> --> 0.2195215
>
> For most Reals x>0, f returns ComplexInfinity, but f[Rationalize[x,0]]
> returns
> a reasonable value.
>
> IMO this behaviour is neither predictable nor desirable.
>
> Peter
>
WorkingPrecision -> $MachinePrecision still makes Mathematica use
significance arithmetic, because $MachinePrecision is just the number
15.95459, which is 'arbitrary' for Mathematica. To use hardware
floating-point arithmetic we need to specify WorkingPrecision ->
MachinePrecision. MachinePrecision is something similar to Pi or E: it is
numeric but not a number.
In version 5.2 we can pass options to NIntegrate like this:
NSum[s[k, 1/4], {k, 0, Infinity},
Method -> {NIntegrate, MaxRecursion -> 15}]
Method -> NIntegrate is slightly different from Method -> Integrate in
that the latter tries to carry out the integration symbolically before
calling NIntegrate.
In this case NIntegrate will complain that the integrand is not numerical
at k = 3.4`12*^1296. The reason is that for large values of k (2*k -
x)/(2*k + x) will be close to 1 and precision control will fail:
In[1]:= 1`20^1`20*^100
Out[1]= Overflow[]
Two possible ways to avoid this is either to set the precision high enough
(so that the base will be distinguishable from 1) or to set $MinPrecision
= $MaxPrecision, in which case precision control seems to work in a
slightly different way:
In[2]:= s[k_, x_] = ((2*k - x)/(2*k + x))^k - Exp[-x];
In[3]:= Module[{prec = 20},
Block[{$MinPrecision = prec, $MaxPrecision = prec},
NSum[s[k, 1/4], {k, 0, Infinity}, WorkingPrecision -> prec]]]
Out[3]= 0.21958908271401804675393479219441155635`20.
This will give a more precise answer:
In[4]:= Block[{$MinPrecision = 70, $MaxPrecision = 70},
NSum[s[k, 1/4], {k, 0, Infinity},
NSumTerms -> 10^4, WorkingPrecision -> 50,
Method -> {NIntegrate, MinRecursion -> 5, MaxRecursion -> 15}]]
Out[4]= 0.2195214776056095400104236039965324040333381067530140390406\
6026688871851600645857299981`70.
With this choice of settings for NIntegrate we don't need to modify
$MinPrecision and $MaxPrecision, but it allows us to get more correct
digits. But also this method has to be used with caution:
In[5]:= Block[{$MinPrecision = 20, $MaxPrecision = 20}, 2^0``20]
Out[5]= 0
For small values of x there is a more efficient way to evaluate the sum:
expand s[k, x] in powers of k and sum term by term.
In[6]:= Module[{x = .25`70},
s[0, x] + Normal[s[k, x] + O[k, Infinity]^70] /.
k^p_ -> Zeta[-p]]
Out[6]= 0.2195214776056095400104236039965324040333381067530140430389\
0299164553992110211438464347`69.57381232331642
So In[4] gave 53 correct digits.
Maxim Rytin
m.r at inbox.ru