neurophysik / jitcdde

Just-in-time compilation for delay differential equations
Other
56 stars 14 forks source link

Updating Solution of DDE at regular Time #23

Closed drsaroj closed 5 years ago

drsaroj commented 5 years ago

Suppose my model equations are

f = [
y(0)*(r0-r2*y(3)-r1*y(0))-alpha*y(0)*y(1)/(1+e1*y(1)),
y(1)*(s0-s2*y(2)-s1*y(1))+(g*alpha*y(0)*y(1))/(1+e1*y(1)),
-delta0*y(2)+sigma*y(1,t-tau)*y(3,t-tau)-beta*y(1)*y(2),
varphi-delta1*y(3)-sigma*y(1)*y(3)+nu*beta*y(1)*y(2),
]

We want to solve the equation with a history say (10,10,10,10) in [0,T]. When time say t1 in [0,T] is of the form n*k where k is fixed and n in non-negative integer, I want to update the solutions obtained by solver at the time just before n*k by the following:

  y(0)=y(0)+(1-p)y(0)
  y(1)=y(1)
  y(2)=y(2)
  y(3)=y(3)

i.e. the change is only in y(0) but all other remains as it is obtained from the solver. The solution to the DDE again will start with the newly changed variable at t1=n*k to next time t2=(n+1)*k where there again changes in the solution by the same above equations.

This solution process will go till final time T for each n=0,1,2,3, ......

The detailed problems are available in the attached file. impulses.pdf

Wrzlprmft commented 5 years ago

The main challenge is this: All DDE solvers require the past to be continuously differentiable and end in the current state and derivative – which again depends on the past, yielding a chicken–egg problem. JiTCDDE further only allows to store the past in a way that complies with this. This is the reason why you need to address initial discontinuities. (I strongly suggest that you read and understand these two sections from the documentation at this point, if you haven’t done so yet.)

What you want to do, introduces a discontinuity per construction. Thus you can only approximate this. The best way is probably the following:

  1. Integrate with a sufficiently small fixed step size (integrate_blindly) exactly up to t₁.
  2. Get the complete state with get_state. This gives you the anchors of the cubic Hermite interpolation which JiTCDDE uses to describe the state.
  3. Add a new anchor at t₁+ε, applying your impulse. Compute the derivative for the new anchor taking into account the impulse.
  4. Set a new past using purge_past and add_past_points.
  5. Repeat.

I will look into the options for providing a service function that does this (and also works with adaptive integration).

drsaroj commented 5 years ago

I know the difficulties in these because the interpolation not only of function but also for its derivatives is used. The method of steps that you described is possible if the solver asks for not only the past history function but also its derivative.

How can I instruct the solver to use my function derivative instead of it's own in every step. Please correct me if I am wrong.

Wrzlprmft commented 5 years ago

First, I forgot a step in the above which I added now; this may clarify things.

You can define the initial past (with derivative) in arbitrary details using add_past_point. This is the lower-function used by constant_past and past_from_function.

Wrzlprmft commented 5 years ago

So, I added a function jump which should address your needs. Please read its entire documentation and test it carefully, as you will be the first one using this. There also is a small example.

To use this, you need to install the latest version from GitHub, e.g., via:

pip3 install git+git://github.com/neurophysik/jitcdde --user