[SciPy-User] Troubles with odeint or ode

[Michael Aye]

I think this example is worth to be kept in the cookbook, what do you
guys think?
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