I want to do some integration and simulation on an ode similar to
x_ddot + d*x_dot + k*x = a*sin(omega*t)
with one modification: I want to replace the external force a * sin(omega * t) by some measurement data which shall be repeated periodicly. This data can for example look like
f = [1,1,1,1,1,1,1,1,1,1,1,1,1,1,1.1,1.3,1.7,1.9,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2.1,2.3,2.9,4.1,4.3,4.4,4.5,4.4,4,3.6,3.2,2.8,2.4,2.0,1.6,1.3,1.2,1.1,1,1,1,1,1]
My code for integrating with pydstool is working perfectly, when I replace the external force by a constant value. In principle it is built up analog to the Calcium channel model (http://www.ni.gsu.edu/~rclewley/PyDSTool/Tutorial/Tutorial_Calcium.html). Using the data as a parameter f and using something like
DSargs.fnspecs = {'force' : (['t'], 'f[int(t)]')
produces the TypeError:
only length-1 arrays can be converted to Python scalars
Is it possible to use discrete data as input-parameter?
I would be very happy if anyone could help me!
Thanks a lot!
Emil
Last edit: Emil Knut 2016-06-29
If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
First of all, thank you very much, my problem is solved now!
I went through example interp_vode_test. I saw how interpolation works, but I didn't find a way to apply it periodicly (except for creating a long vector, containing the same data again and again first of all). But I found another, in my case even better solution: Combining if-statements with the mod() function (which gives the rest of a division), for example:
You found a bug in the parsing of embedded if statements for C-based integrators. Could you update your PyDSTool to the latest version on github, or at least grab PyDSTool/core/codegenerators/c.py? You should be able to use Radau with your setup now, which is the best available for stiff systems.
Although, I don't really understand the motivation for your choice of auxiliary function. It's not periodic, if that's what you were going for. You should plot it (e.g. from using VODE) and see for yourself.
If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
The updated tool works fine! Dopri and Radau both are performing well. I have one more question, i'm plotting results using traj.sample() and I've been playing around with the options dt, doEvents and precise. But the resulting sample can naver be more precise than the underlaying time mesh, which depends on the integrator, where vode produces a very coarse mesh and radau a very fine mesh. How can I take influence on this time mesh? I've tried Vode and Dopri in combination with
DSargs.algparams = {'max_step': 0.0001}
which didn't make anython better.
Best regards,
Emil
Last edit: Emil Knut 2016-07-05
If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
The updated tool works fine! Dopri and Radau both are performing well. I
have one more question, i'm plotting results using traj.sample() and I've
been playing around with the options dt, doEvents and precise. But the
resulting sample can naver be more precise than the underlaying time mesh,
which depends on the integrator, where vode produces a very coarse mesh and
radau a very fine mesh. How can I take influence on this time mesh?
I think I'm very close to my aim now. Integrating my system works perfectly, it gives the same results in comparison to my old tool using odeint. Now I want to start doing some bifurcation analysis of my system using LC-C. As I know the exact period-duration for the initial free-parameter value, it should be easy to generate the initcycle: Integrate for example 30 periods and use the last period (=tlo...thi) as "cycle"
In some cases I get 10-20 bifurcation points (containing RG, BP, LP) but then the output turns to NaN NaN NaN....'Warning: NaNs in solution from AUTO.'
In some other cases bifurcation stops at an MX, sometimes after finding some bif-points, sometimes right at the beginning of bifurcation.
Do you have an idea, whats wrong here? Which problems are causing NaN and MX? Unfortunately I can't give you my system, but I guess I'm tradiing with a stiff system here, depending on the parameters.
I would really appreciate another answer, sorry for asking that much. :-/
Thanks a lot!
Last edit: Emil Knut 2016-07-11
If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
To add to Drew's points: If you plan to do continuation and
bifurcation analysis, you are unlikely to succeed if you have a system
with a non-smooth external input as a time-varying parameter. Those
algorithms expect smooth systems only. At least, you should consider
fitting a combo of trig functions to your input signal and using that
instead.
I think I'm very close to my aim now. Integrating my system works perfectly,
it gives the same results in comparison to my old tool using odeint. Now I
want to start doing some bifurcation analysis of my system using LC-C. As I
know the exact period-duration for the initial free-parameter value, it
should be easy to generate the initcycle: Integrate for example 30 periods
and use the last period (=tlo...thi) as "cycle"
In some cases I get 10-20 bifurcation points (containing RG, BP, LP) but
then the output turns to NaN NaN NaN....'Warning: NaNs in solution from
AUTO.'
In some other cases bifurcation stops at an MX, sometimes after finding some
bif-points, sometimes right at the beginning of bifurcation.
Do you have an idea, whats wrong here? Which problems are causing NaN and
MX? Unfortunately I can't give you my system, but I guess I'm tradiing with
a stiff system here, depending on the parameters.
I would really apreciate another answer, sorry for asking that much. :-/
If you are able to continue and find bifurcation points, then the number of points for your initial cycle is fine. What values of NumIntervals have you tried? For stiff systems, I've had to go upwards of 1000 on this one.
If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
Both hints were helpful, increasing NumIntervals helped somewhat, I went from 40 to higher Numbers, testing up to 2000. But for some parameters I reached MX still to fast. Then I tried to replace my aux-function by a dummy funtion = f(sin(t)), so that there is no if-statement in my functions any more. Now continuation seems to work! - Thank you guys!
As usual, there is an other issue now: Limiting the free parameter seems not to work and the free parameter doesn't want to move to values which are high enough. The free parameter is only turning lower, not higher in comparison to initial value. I've tried it like described in examples:
Hi, Emil. Not sure I understand what is going on now. What do you mean by the free parameter doesn't want to move high enough? Do you mean that you can only continue in one direction? That is strange. Can you switch the order and do forward before backward and report back? Thanks.
If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
If your step size is tiny and your num_points is small, you may have
to repeat calls to .forward() many times to make O(1) size changes
in the parameter.
Hi, Emil. Not sure I understand what is going on now. What do you mean by
the free parameter doesn't want to move high enough? Do you mean that you
can only continue in one direction? That is strange. Can you switch the
order and do forward before backward and report back? Thanks.
I increased step-size and used .forward() several times. But the parameter stays always close to the initial value (I've tried different initial parameter-values). I'm getting a lot of BPs and LPs, but no PDs or LPCs. How is this possible, as I'm doing LC-C?
Thanks,
Emil
Last edit: Emil Knut 2016-07-20
If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
after trying a lot, I think I'm able to describe my problem better now. If you have a look on example PyCont_Lorenz, the line in the bifurcation diagram is limited by an LPC on the right-hand side. But what happens for higher values of the parameter? The same problem was discussed here:
In his final script, R.A. takes the LPC, calculates 0.5 * period-duration and starts continuation again. The LPC is detected as PD now an continuation goes on. But doing this for each LPC might be very complicated, as my problem shows many LPCs. Isn't there a better way to to this kind of bifurcation?
Thank you very much,
Emil
If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
Last try:
I want to examine some periodic doubling phenomena. Forward and backward continuation is working, it gives me a line in bif.diagram, where stability is lost in a periodic doubling point, comparable to the example PyCont_Lorenz:
Now I want to continue in PD1 with double period, which should give me a stable solution again, compareable to this example from literature:
which is not working, PCargs.period =... seems not to work, PyCont['LC2'].forward() starts with old period in PD1 and therefore LC2 gets identical with LC1.
Hope anyone can help me.
Thank you!
Emil
Last edit: Emil Knut 2016-07-27
If you would like to refer to this comment somewhere else in this project, copy and paste the following link:
Hi to all!
I want to do some integration and simulation on an ode similar to
x_ddot + d*x_dot + k*x = a*sin(omega*t)with one modification: I want to replace the external force a * sin(omega * t) by some measurement data which shall be repeated periodicly. This data can for example look like
f = [1,1,1,1,1,1,1,1,1,1,1,1,1,1,1.1,1.3,1.7,1.9,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2.1,2.3,2.9,4.1,4.3,4.4,4.5,4.4,4,3.6,3.2,2.8,2.4,2.0,1.6,1.3,1.2,1.1,1,1,1,1,1]My code for integrating with pydstool is working perfectly, when I replace the external force by a constant value. In principle it is built up analog to the Calcium channel model (http://www.ni.gsu.edu/~rclewley/PyDSTool/Tutorial/Tutorial_Calcium.html). Using the data as a parameter f and using something like
DSargs.fnspecs = {'force' : (['t'], 'f[int(t)]')produces the TypeError:
only length-1 arrays can be converted to Python scalarsIs it possible to use discrete data as input-parameter?
I would be very happy if anyone could help me!
Thanks a lot!
Emil
Last edit: Emil Knut 2016-06-29
You could allow the data to be linearly interpolated, if you follow a scheme like that shown here: https://sourceforge.net/p/pydstool/discussion/472291/thread/ceddb4c0/. The included
examples/interp_vode_test.pyalso shows a basic working example. Does that solve your problem?Hello Rob!
First of all, thank you very much, my problem is solved now!
I went through example interp_vode_test. I saw how interpolation works, but I didn't find a way to apply it periodicly (except for creating a long vector, containing the same data again and again first of all). But I found another, in my case even better solution: Combining if-statements with the mod() function (which gives the rest of a division), for example:
wfn_str = 'auxval1(t)'
fnspecs = {'auxval1': (['t'], 'if(mod(t,3)<1,-t,if(mod(t,3)<2,0,3))')
produces a periodic RHS. As I'm not using measurement-data directly, this works perfect for me! Integration works for this example with
testODE = Vode_ODEsystem(DSargs)
but
testODE = Dopri_ODEsystem(DSargs)
gives the error
SyntaxError: Mismatch of braces in expression rhs_if(fmod(t,3)<1,-t,if(fmod(t,3)<2,0,3))rhs_if(fmod(t,3)<2,0,3))
...why is this? Which integrator do you recommend for stiff ODE-systems? Is there a list, showing the available integrators?
Thank you very much for your fast help and you great tool!
Best regards
Emil
You found a bug in the parsing of embedded
ifstatements for C-based integrators. Could you update your PyDSTool to the latest version on github, or at least grabPyDSTool/core/codegenerators/c.py? You should be able to use Radau with your setup now, which is the best available for stiff systems.Although, I don't really understand the motivation for your choice of auxiliary function. It's not periodic, if that's what you were going for. You should plot it (e.g. from using VODE) and see for yourself.
Thank you for your help rob!
The updated tool works fine! Dopri and Radau both are performing well. I have one more question, i'm plotting results using traj.sample() and I've been playing around with the options dt, doEvents and precise. But the resulting sample can naver be more precise than the underlaying time mesh, which depends on the integrator, where vode produces a very coarse mesh and radau a very fine mesh. How can I take influence on this time mesh? I've tried Vode and Dopri in combination with
DSargs.algparams = {'max_step': 0.0001}
which didn't make anython better.
Best regards,
Emil
Last edit: Emil Knut 2016-07-05
It's all here, in the controls for time step:
http://www.ni.gsu.edu/~rclewley/PyDSTool/Generators.html
On Tue, Jul 5, 2016 at 3:11 PM, Emil Knut pydsuser@users.sf.net wrote:
Again, thank you for your help!
I think I'm very close to my aim now. Integrating my system works perfectly, it gives the same results in comparison to my old tool using odeint. Now I want to start doing some bifurcation analysis of my system using LC-C. As I know the exact period-duration for the initial free-parameter value, it should be easy to generate the initcycle: Integrate for example 30 periods and use the last period (=tlo...thi) as "cycle"
DSargs.algparams = { 'max_pts': 1000000} DS = Generator.Radau_ODEsystem(DSargs) traj = DS.compute('demo') cycle = traj.sample(dt=0.00001, tlo, thi, precise=True)Based on this, I'm using
for continuation. Unfortunately the continuation doesn't work properly. I've spent a lot of time playing with the following options:
In some cases I get 10-20 bifurcation points (containing RG, BP, LP) but then the output turns to NaN NaN NaN....'Warning: NaNs in solution from AUTO.'
In some other cases bifurcation stops at an MX, sometimes after finding some bif-points, sometimes right at the beginning of bifurcation.
Do you have an idea, whats wrong here? Which problems are causing NaN and MX? Unfortunately I can't give you my system, but I guess I'm tradiing with a stiff system here, depending on the parameters.
I would really appreciate another answer, sorry for asking that much. :-/
Thanks a lot!
Last edit: Emil Knut 2016-07-11
To add to Drew's points: If you plan to do continuation and
bifurcation analysis, you are unlikely to succeed if you have a system
with a non-smooth external input as a time-varying parameter. Those
algorithms expect smooth systems only. At least, you should consider
fitting a combo of trig functions to your input signal and using that
instead.
On Mon, Jul 11, 2016 at 9:08 AM, Emil Knut pydsuser@users.sf.net wrote:
If you are able to continue and find bifurcation points, then the number of points for your initial cycle is fine. What values of NumIntervals have you tried? For stiff systems, I've had to go upwards of 1000 on this one.
Thank you very much Drew and Rob!
Both hints were helpful, increasing NumIntervals helped somewhat, I went from 40 to higher Numbers, testing up to 2000. But for some parameters I reached MX still to fast. Then I tried to replace my aux-function by a dummy funtion = f(sin(t)), so that there is no if-statement in my functions any more. Now continuation seems to work! - Thank you guys!
As usual, there is an other issue now: Limiting the free parameter seems not to work and the free parameter doesn't want to move to values which are high enough. The free parameter is only turning lower, not higher in comparison to initial value. I've tried it like described in examples:
pardict['alpha'] = 4. ... DSargs.pars = pardict ... DSargs.pdomain = {'alpha':[1., 8.]} ... PCargs.freepars = ['alpha'] ... PCargs.StopAtPoints = 'B' ... PyCont['LC1'].backward() PyCont['LC1'].forward()I've tried Dopri and Radau.
Any idea what I'm doing wrong here?
Thank you!
Emil
Last edit: Emil Knut 2016-07-12
Hi, Emil. Not sure I understand what is going on now. What do you mean by the free parameter doesn't want to move high enough? Do you mean that you can only continue in one direction? That is strange. Can you switch the order and do forward before backward and report back? Thanks.
If your step size is tiny and your num_points is small, you may have
to repeat calls to
.forward()many times to make O(1) size changesin the parameter.
On Tue, Jul 12, 2016 at 2:09 PM, Drew LaMar mdlamar@users.sf.net wrote:
I increased step-size and used .forward() several times. But the parameter stays always close to the initial value (I've tried different initial parameter-values). I'm getting a lot of BPs and LPs, but no PDs or LPCs. How is this possible, as I'm doing LC-C?
Thanks,
Emil
Last edit: Emil Knut 2016-07-20
Hi again!
after trying a lot, I think I'm able to describe my problem better now. If you have a look on example PyCont_Lorenz, the line in the bifurcation diagram is limited by an LPC on the right-hand side. But what happens for higher values of the parameter? The same problem was discussed here:
https://sourceforge.net/p/pydstool/discussion/472291/thread/aa15f0ee/
In his final script, R.A. takes the LPC, calculates 0.5 * period-duration and starts continuation again. The LPC is detected as PD now an continuation goes on. But doing this for each LPC might be very complicated, as my problem shows many LPCs. Isn't there a better way to to this kind of bifurcation?
Thank you very much,
Emil
This was nonsense again - sorry for that.
Last try:
I want to examine some periodic doubling phenomena. Forward and backward continuation is working, it gives me a line in bif.diagram, where stability is lost in a periodic doubling point, comparable to the example PyCont_Lorenz:
Now I want to continue in PD1 with double period, which should give me a stable solution again, compareable to this example from literature:
How can this be done? I've tried
PCargs.name = 'LC2' PCargs.initpoint = 'LC1:PD1' PCargs.period = PyCont['LC1'].getSpecialPoint('PD1')['_T']*2 PyCont.newCurve(PCargs) PyCont['LC2'].forward()which is not working, PCargs.period =... seems not to work, PyCont['LC2'].forward() starts with old period in PD1 and therefore LC2 gets identical with LC1.
Hope anyone can help me.
Thank you!
Emil
Last edit: Emil Knut 2016-07-27