Mathematica 9 is now available
Services & Resources / Wolfram Forums
-----
 /
MathGroup Archive
2006
*January
*February
*March
*April
*May
*June
*July
*August
*September
*October
*November
*December
*Archive Index
*Ask about this page
*Print this page
*Give us feedback
*Sign up for the Wolfram Insider

MathGroup Archive 2006

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: Periodic Rebirth of Hyperbolic Functions by ODE in Mathematica

  • To: mathgroup at smc.vnet.net
  • Subject: [mg66003] Re: Periodic Rebirth of Hyperbolic Functions by ODE in Mathematica
  • From: Maxim <m.r at inbox.ru>
  • Date: Wed, 26 Apr 2006 04:38:17 -0400 (EDT)
  • References: <e2i7qg$8tg$1@smc.vnet.net>
  • Sender: owner-wri-mathgroup at wolfram.com

On Mon, 24 Apr 2006 10:04:00 +0000 (UTC), Narasimham  
<mathma18 at hotmail.com> wrote:

> Solutions of non-linear second order differential equations formed for
> r = Sech[th],Tanh[th] after two differentiations  are repeated  at
> about th=8.0 which value should actually be Infinity. After this value
> they do not tend to the expected asymptotic value but become  periodic
> functions !.
>
> The NDSolve Runge-Kutta numerical algorithm appears to involve a small
> positive error at large argument values that allows build-up of
> function back to its full value as a reflected periodic function.This
> is seen more clearly in a polar plot of r = Sech[th] and Tanh[th] where
> periodic lobes reappear, simulating  periodic functions ( 0.02 added in
> plots just to distinguish the curves).
>
> BTW, I would be much enthused if one can identify the ( Jacobi Elliptic
> like ??) functions labeled as Reborn, that result by this Mathematica
> command/program.  Or please indicate whether there is a command to
> convert asymptotic functions with introduced second order errors to
> their periodic counterparts where  the time period can be numerically
> defined /specified.  It appears as an unacceptable error for Sech or
> Tanh because the hyperbolic monotonic nature of the non-linear second
> order differential equation gets  numerically depicted  or compromised
> as an elliptic periodic one.
>
> Regards,  Narasimham
>
> "-- Sech --"
> Clear[r,th,R];
> eq = { r''[th] == r[th]*(1 - 2 r[th]^2), r[0]==1, r'[0]==0};
> NDSolve[ eq, r,{th,0,15}] ;
> R[th_]=r[th]/.First[%];
> Reborn=Plot[%,{th,0,15}];
> soln=Plot[Sech[th]+.02,{th,0,15}];
> Show[Reborn,soln];
> {x,y}={R[th]Cos[th],R[th]Sin[th]};
> ParametricPlot[{x,y},{th,0,15},AspectRatio->Automatic];
> "-- Tanh --"
> Clear[r,th,R];
> eq = { r''[th] == 2*(r[th]^3 - r[th]), r[0]==0, r'[0]==1};
> NDSolve[ eq, r,{th,0,38}] ;
> R[th_]=r[th]/.First[%];
> Reborn=Plot[%,{th,0,38}];
> soln=Plot[Tanh[th]+.02,{th,0,38}];
> Show[Reborn,soln];
> {x,y}={R[th]Cos[th],R[th]Sin[th]};
> ParametricPlot[{x,y},{th,0,15},AspectRatio->Automatic];
>

One way to obtain a more accurate solution is to figure out an invariant  
of the system:

In[1]:= inv = Integrate[(r''[th] - r[th]*(1 - 2*r[th]^2))*r'[th], th]

Out[1]= (-(1/2))*r[th]^2 + r[th]^4/2 + (1/2)*Derivative[1][r][th]^2

Hence this quantity should remain constant. Then

sol = First@ NDSolve[
   {r''[th] - r[th]*(1 - 2*r[th]^2) == 0, r[0] == 1, r'[0] == 0},
   r, {th, 0, 100},
   Method -> {Projection, Invariants -> inv},
   MaxStepFraction -> .001]

Plot[r[th] - Sech[th] /. sol, {th, 0, 100}, PlotRange -> All]

shows that the error remains bounded (in fact, the accuracy even increases  
linearly). We could also get the exact solution from the system {inv == 0,  
r[0] == 1} but in general that may be problematic, as the solution to that  
system is not unique, there is another solution, namely (r -> (1&)).

This method doesn't work very well for your second system of ODEs. In that  
case we can use (1 - r[th]^2)/r'[th] as an invariant but I don't know how  
to obtain this form without guessing.

Maxim Rytin
m.r at inbox.ru


  • Prev by Date: Compile&Which
  • Next by Date: Re: error function
  • Previous by thread: Periodic Rebirth of Hyperbolic Functions by ODE in Mathematica
  • Next by thread: Object-Oriented Paradigm in Mathematica?