Re: False divergence of the NDSolve solution: how to avoid

*To*: mathgroup at smc.vnet.net*Subject*: [mg101842] Re: False divergence of the NDSolve solution: how to avoid*From*: Roland Franzius <roland.franzius at uos.de>*Date*: Sun, 19 Jul 2009 07:14:14 -0400 (EDT)*References*: <h3n5uh$2gp$1@smc.vnet.net> <h3pesm$oib$1@smc.vnet.net>

Roland Franzius schrieb: > Alexei Boulbitch schrieb: >> Dear Community, >> >> I am simulating a system of ODE using v6. Here are the equations: >> >> eq1 = x'[t] == y[t]; >> eq2 = y'[t] == 1/x[t] - 1.4 - (4.5 + y[t])*(1 + z[t]^2); >> eq3 = z'[t] == 18*z[t] - 0.75*(4.5 + y[t])^2*z[t] - z[t]^3; >> >> It is simulated at x>0. This system at x>0 seems to be globally stable. >> To understand it observe that at large x, y, and z one finds >> y' ~ - y*z^2 and z' ~ - z^3. In other words, there is a kind of a >> non-linear "returning force" for y and z, while x follows the dynamics >> of y. Since x occurs in eq2 this is not correct. Elimination of dummy y gives x'' ~ x' z^2 and z'= -z^3 but the problem is at x=0 x'' = 1/x z' = z( 18 -0.75(4.5+x') -z^2) If by error you have x<0 and x' <4.5 you have a problem. >> However, when solving it on Mathematica I sometimes find trajectories >> that counterintuitively diverge. >> Check this for example: >> >> NDSolve[{eq1, eq2, eq3, x[0] == 0.669, y[0] == 0.881, >> z[0] == 0.988}, {x, y, z}, {t, 0, 40}]; >> Plot[{Evaluate[x[t] /. s], Evaluate[y[t] /. s], >> Evaluate[z[t] /. s]}, {t, 0, 45}, PlotRange -> All, >> PlotStyle -> {Red, Green, Blue}, >> AxesLabel -> {Style["t", 16], Style["x, y, z", 16]}] >> My guess is that this is due to some peculiarity in the numeric method >> used, and the method should be probably changed, or its parameters >> specified. I am however, not experienced in numeric approaches for >> solving ODEs. >> >> Now comes the question: >> Can you give me a hint, of >> (i) what may be the reason of such a behavior? >> and >> (ii) What should I do to avoid such a false divergence? > > Its always extremely difficult to control solutions of osciallating > nonlinear systems for a long time:-( > > At least in this case use Abs[z[t]] on the right to cross the line > immediately if an error of z<0 is occurs or choose a t-step adaptive > method if z is approaching 0. > ############# After a hint of Alois Steindl I reconsidered the full second order system once again. It is possible to integrate beyond t=28 if x is restricted to x>0 somehow. The second order system eq1 = x''[t] == 1/x[t] - 1.4 - (4.5 + x'[t])*(1 + z[t]^2) eq2 = z'[t] == 18*z[t] - 0.75*(4.5 + x'[t])^2 *z[t] - z[t]^3 explodes numerically somewhere at t=28 s = NDSolve[{eq1, eq2, x[0] == 0.669, x'[0] == 0.881, z[0] == 0.988}, {x, z}, {t, 0, 28.5}]; To see the point of explosion at a logarithmic scale look at Plot[{Evaluate[ArcSinh[x[t] /. s]], Evaluate[ArcSinh[z[t] /. s]]}, {t, 27., 28.1}, PlotRange -> All, PlotStyle -> {Red, Green}, AxesLabel -> {Style["t", 16], Style["x, y, z", 16]}] The reason is x becoming negative and a moment later z follows. This is analytically impossible because of the hard reflection at x=0 by 1/x. Numerically of course a big z-term with x'<-4.5 may override it. A try with absolute values eq1 = x''[t] == 1/x[t] - Abs[1.4 - (4.5 + x'[t])*(1 + z[t]^2)] eq2 = z'[t] == 18*z[t] - 0.75*(4.5 + x'[t])^2 *z[t] - z[t]^3 s = NDSolve[{eq1, eq2, x[0] == 0.669, x'[0] == 0.881, z[0] == 0.988}, {x, z}, {t, 0, 40}]; Plot[{Evaluate[x[t] /. s], Evaluate[z[t] /. s]}, {t, 0, 3}, PlotRange -> All, PlotStyle -> {Red, Green, Blue}, AxesLabel -> {Style["t", 16], Style["x, y, z", 16]}] The apparent singularities at x=0 is caused by the hard reflection of x at x=0 by the non-Lipshitz x'' = 1/x force term. The Abs-approach seamingly eliminates the problem for the given starting values. If it is working for all starting values, who knows. -- Roland Franzius

**Re: Player Pro Limitations**

**Re: eigenvector centrality**

**Re: False divergence of the NDSolve solution: how to avoid**

**eigenvector centrality**