[SciPy-User] Troubles with odeint or ode

[Jim Vickroy]

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... 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
>   

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20091104/65674cc3/attachment.html>