[SciPy-User] Troubles with odeint or ode

[Anne Archibald]

2009/11/4 Jim Vickroy <Jim.Vickroy at noaa.gov>:
> Michael Aye wrote:
>
> I think this example is worth to be kept in the cookbook, what do you
> guys think?

In the short term it's probably worth putting something like this in
the cookbook (though don't use the OP's code without permission!) but
really the solution is to fix odeint (which shouldn't be too
difficult). I put a ticket in Trac, but haven't had time to actually
fix the problem yet.

Anne

> +1
>
> I thought the same but did not speak up.
> -- jv
>
> It is showing how to do the oscillator modelling and can show the
> 'danger' of lack of complex support of odeint.
>
> Just my feeling it could be quite helpful.
> I myself certainly copied this into my treasure of worth-to-remember
> evernotes.. ;)
>
> Regards,
> Michael
>
> On Nov 3, 8:28 pm, Anne Archibald <peridot.face... at gmail.com> wrote:
>
>
> 2009/11/3 ANDREA ARMAROLI <rmr... at unife.it>:
>
>
>
> Dear users,
>
>
> I'm trying to solve an ODE system that models a parametric oscillator with
> complex amplitudes at two freqencies.
>
>
> I'm new to python. This problem is simply solved in matlab using ode45.
>
>
> If I try to use odeint I cannot set parameters like tolerances or max step.
>
>
> You can in fact control these using the (optional) parameters hmax,
> rtol, and atol.
>
>
>
> The result is constant-valued solutions.
>
>
> After a bit of experimentation (in particular, a print statement
> inside your derivative function) it turns out that the problem is that
> odeint does not support complex values (it silently discards imaginary
> parts). This is not a major obstacle, since you can just pack and
> unpack the values yourself. When I do that, the plot I get is two
> oscillatory results (the absolute values), and one nice flat line (for
> the total, which I'm guessing you need to be conserved).
>
> I have to say, it would be very good if odeint either reported an
> error or worked with complex values, so you would have found this
> easier to track down. But it looks like it does a reasonable job of
> solving your problem once you work around its lack of complex support.
>
> Anne
>
>
>
>
>
> Trying with ode class and ZVODE integrator, I have that total intensity is
> not
> conserved. I have weird quasi-periodic oscillations.
>
>
> I do know that my equations have excess degrees of freedom and I know from
> symmetries what the integrals of motion are, but in Matlab this works pretty
> well.
>
>
> Here is the code with odeint
>
>
> import numpy as N
> import scipy as S
> import scipy.integrate
> import pylab as P
>
>
> def deriv(Y,t):
>
>
>    A1 = Y[0]
>    A2 = Y[1]
>    A1d = 1j*A2*N.conj(A1)
>    A2d = 1j*(A1**2/2.0 + dk*A2)
>    return [A1d, A2d]
>
>
> nplot = 10000
> zmax = 10.0
> zstep = zmax/nplot
>
>
> dk = 0.5 # normalised detuning/dispersion
> phi0 = 0.0*N.pi # initial dephasing
> eta0 = 0.3 # initial pump intensity
>
>
> # initial values
> u20 = N.sqrt(eta0)*N.exp(1j*phi0)
> u10 = N.sqrt(2*(1-eta0))
>
>
> H0=dk*eta0+2*N.sqrt(eta0)*(1-eta0)*N.cos(phi0);
>
>
> Y,info =
> scipy.integrate.odeint(deriv,[u10,u20],N.arange(0,zmax+zstep,zstep),full_ou
> tput=True,printmessg
> = True)
> # are conserved quantities constant?
> Ptot = N.abs(Y[:,0])**2/2+ N.abs(Y[:,1])**2
> P.plot(N.arange(0,zmax+zstep-0.0001,zstep),Ptot)
> # then Hamiltonian...
> #...
>
>
> p1  = P.plot(N.arange(0,zmax+zstep-0.0001,zstep),N.abs(Y[:,0]))
> p2  = P.plot(N.arange(0,zmax+zstep-0.0001,zstep),N.abs(Y[:,1]))
>
>
> Thank you very much for your help.
>
>
> Andrea Armaroli
>
>
> _______________________________________________
> SciPy-User mailing list
> SciPy-U... at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>
>
>  foo.py
> 1KViewDownload
>
> _______________________________________________
> SciPy-User mailing list
> SciPy-U... at scipy.orghttp://mail.scipy.org/mailman/listinfo/scipy-user
>
>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>
>
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>
>