Re: Errors in Mathematica

*To*: mathgroup at smc.vnet.net*Subject*: [mg69258] Re: Errors in Mathematica*From*: Peter Pein <petsie at dordos.net>*Date*: Sun, 3 Sep 2006 23:47:13 -0400 (EDT)*References*: <ed3qkr$rnh$1@smc.vnet.net>*Sender*: owner-wri-mathgroup at wolfram.com

Chris schrieb: > NDSolve[{P'[x] == Q[x], Q'[x] == 10^(-10) , > P[10^(-6)]==0, Q[10^(-6)] == 0, > n[x] == Exp[-A[x]] , > A'[x] == -A[x]/236, > A[10^(-6)]== -.0001}, > {P,Q,A,n}, {x, 10^(-6),1000}, > MaxSteps rightarrow Infinity] > Plot[Evaluate[P[x]/. %], {x,10^(-6), 10}] > > Here 2 equations link the functions P and Q. > Two other equations link the functions n and A. > There is no coupling between the two sets, and yet: > > Note the denominator 236 that appears in the second set. > As this number is changed the solution for P[x], and in particular the > value P[6] changes. > Of course this should not happen. > There is a discontinuity between 236 and 237 and another between 274 > and 275. > After that, further increase has no effect on P[6]. > > The problem goes away if the small numbers are increased. > It also goes away if the exponential function is replaced by unity or > some elementary function. > > If the small numbers are made much smaller it causes Mathematica to > shut down. > > Can anyone tell me how to handle this? > Hi Chris, I guess the change from 236 to 237 changes the stepsize of the algorithm at some point. You solve for (approx.) the range 0..1000 but use the value at 6. Using a different StartingStepSize eliminates this problem: P6 gives P[6] as a function of the 'magic' denominator and the StartingStepSize: P6[denom_, sss_:Automatic] := P[6] /. First[NDSolve[{P'[x] == Q[x], Q'[x] == 10^(-10), P[10^(-6)] == 0, Q[10^(-6)] == 0, n[x] == Exp[-A[x]], A'[x] == -A[x]/denom, A[10^(-6)] == -10^(-4)}, {P, Q, A, n}, {x, 10^(-6), 1000}, MaxSteps -> Infinity, StartingStepSize -> sss]] P6[237] - P6[236] --> 3.896888607874218*^-11 this is the effect, you have observed. Let's try a smaller sss: P6[237, 10^(-6)] - P6[236, 10^(-6)] --> 0. To see how small the stepsize has to be: Plot[Log[10, $MachineEpsilon + Abs[P6[237, Re[10^(-x)]] - P6[236, Re[10^(-x)]]]], {x, 0, 1}, AxesOrigin -> {0, Log[10, $MachineEpsilon]}]; The plot shows that for sss ~10^0.7 (ca. 5) and smaller values the effect vanishes. HTH, Peter