[SciPy-User] Troubles with odeint or ode

[Anne Archibald]

2009/11/3 ANDREA ARMAROLI <rmrndr 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_output=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-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: foo.py
Type: text/x-python
Size: 1235 bytes
Desc: not available
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20091103/d0a63da4/attachment.py>