[SciPy-User] Troubles with odeint or ode

ANDREA ARMAROLI rmrndr@unife...
Thu Nov 5 04:17:15 CST 2009


Hi guys,

I think it is worth placing in the cookbook. If you like, I can provide a few details of the framework it comes from.

Thank you VERY much for the quick response.

Andrea.
---------- Original Message -----------
From: Jim Vickroy <Jim.Vickroy@noaa.gov> 
To: SciPy Users List <scipy-user@scipy.org> 
Sent: Wed, 04 Nov 2009 08:03:39 -0700 
Subject: Re: [SciPy-User] Troubles with odeint or ode

> Michael Aye wrote:I think this example is worth to be kept in the cookbook, what do

you
guys

think?

+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...@gmail.com>

wrote:

2009/11/3 ANDREA ARMAROLI <rmr...@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...@scipy.org
http://mail.scipy.org/mailman/listinfo/scipy-user

 foo.py
1KViewDownload

_______________________________________________
SciPy-User mailing

list
SciPy-U...@scipy.orghttp://mail.scipy.org/mailman/listinfo/scipy-user

_______________________________________________
SciPy-User mailing

list
SciPy-User@scipy.org
http://mail.scipy.org/mailman/listinfo/scipy-user

------- End of Original Message -------
 
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mail.scipy.org/pipermail/scipy-user/attachments/20091105/cb75d565/attachment.html 


More information about the SciPy-User mailing list