[SciPy-User] Troubles with odeint or ode
Wed Nov 4 09:03:39 CST 2009
Michael Aye wrote:
> I think this example is worth to be kept in the cookbook, what do you
> guys think?
I thought the same but did not speak up.
> 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.. ;)
> 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.
>>> 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
>>> A2 = Y
>>> 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))
>>> 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
>>> # 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 mailing list
> SciPy-User mailing list
-------------- next part --------------
An HTML attachment was scrubbed...
More information about the SciPy-User