0

I'm trying to solve a 2D delay differential equation with some parameters. The problem is that I can’t get the right solution (which I know) and I suspect that it comes from the integration step, but I'm not sure and I don't really understand how the JiTCDDE works.

This is the DDE:

enter image description here

This is my model:

def model(p, q, r, alpha, T, tau, tmax, ci):
    f = [1/T * (p*y(0)+alpha*y(1, t-tau)), 1/T * (r*y(0)+q*y(1))]
            
    DDE = jitcdde(f)
    
    
    DDE.constant_past(ci)
    
    DDE.step_on_discontinuities()
    
         
    data = []
    for time in np.arange(DDE.t, DDE.t+tmax, 0.09):
        data.append( DDE.integrate(time)[1])
    return data

Where I'm only interested in the y(1) solution

And the parameters:

T=32        #escala temporal
p=-2.4/T
q=-1.12/T
r=1.5/T
alpha=.6/T
tau=T*2.4     #delay
tmax=400
ci = np.array([4080, 0])

This is the plot I have with that model and parameters:

enter image description here

And this is (the blue line) the correct solution (someone give me the plot not the data)

enter image description here

Wrzlprmft
  • 4,234
  • 1
  • 28
  • 54
DiegoT
  • 29
  • 6
  • Can please [edit] your question to give us some more information as to what problem this is and what the correct solution is (so we can compare)? Right now, all we know is that it’s not doing what you want, which is not much to work with. – Wrzlprmft Dec 12 '20 at 21:53
  • That being said, an apparent problem is that `model` only returns values but not the time points at which these values occur. – Wrzlprmft Dec 12 '20 at 21:54
  • I havent the correct solution in array, but it plotted is like the las image – DiegoT Dec 12 '20 at 22:04
  • I want an array that i can plot and give me the curve of the last image I'm getting something like a parabola The problem is – DiegoT Dec 12 '20 at 22:07
  • from what I see, the DDE.integrate only gives the value of the functions (y(0) and y(1)) – DiegoT Dec 12 '20 at 22:08
  • *from what I see, the DDE.integrate only gives the value of the functions (y(0) and y(1))* – Well, it doesn’t return the time because you already give this as an input. It’s simply `time`. – Wrzlprmft Dec 12 '20 at 22:32
  • Note that you get a solution that is more similar to your control when you integrate for a longer time (`tmax=10000`). Since your control uses “arbitrary units”, this may be the cause. However, you don’t get a second local maximum. I strongly doubt that your control solution is the correct solution of the DDE in your code. However, without further information (e.g., what equation are you trying to implement, why you think the control is correct, …), it’s impossible to say where the problem is. You might as well have a typo in your definition of `f`. – Wrzlprmft Dec 12 '20 at 22:37
  • There was a wrong 1/T in the equation parameters. The control solution is correct because it was integrated in matlab by someone with more expertice in solving delay differential equations and he give me those parameters to me so I can get the same solution to te DDE by myself as a practice of data modeling. The model is a 2D DDE which represent the evolution of some volume in time. Now I'm getting a UnsuccessfulIntegration error and I think is because the DDE is ill-posed, How do I handle thing like this ? Is the step or the integration parameters some important to this case ? – DiegoT Dec 13 '20 at 23:46
  • If DDE23 successfully integrates your DDE, but your JiTCDDE implementation doesn’t, the most likely case is that there is a further mistake in your implementation of the derivative. But without knowing what the derivative is supposed to be or your code, it’s impossible to tell. Why don’t you just show us these things? – Wrzlprmft Dec 14 '20 at 08:14
  • What you mean with "what the derivative is supposed to be or your code" ? I almost copy the Matlab DDE but with Python sintaxis: `ddex1de= @(t,y,ylag) [1/T*(p*y(1)+alpha*ylag(2)) ; 1/T*(r*y(1)+q*y(2))]` – DiegoT Dec 14 '20 at 13:16
  • *What you mean with "what the derivative is supposed to be or your code" ?* – Right now, we do not know the intended derivative (equations, not code). We also do not know what your current code (producing the `UnsuccessfulIntegration`) looks like. – Wrzlprmft Dec 14 '20 at 13:18
  • I've edited the question with the equations – DiegoT Dec 14 '20 at 14:34

1 Answers1

0

The following code works for me and produces a result similar to your control:

import numpy as np
from jitcdde import y,t,jitcdde

T = 1
p = -2.4/T
q = -1.12/T
r = 1.5/T
alpha = .6/T
tau = T*2.4
tmax = 10
ci = np.array([4080, 0])

f = [
        1/T * (p*y(0)+alpha*y(1, t-tau)),
        1/T * (r*y(0)+q*y(1))
    ]

DDE = jitcdde(f)
DDE.constant_past(ci)
DDE.adjust_diff()

times = np.linspace( DDE.t, DDE.t+tmax, 1000 )
values = [DDE.integrate(time)[1] for time in times]

from matplotlib.pyplot import subplots
fig,axes = subplots()
axes.plot(times,values)
fig.show()

Note the following:

  • I set T=1 (and adjusted tmax accordingly). I presume that there still is a mistake here.

  • I used adjust_diff instead of step_on_discontinuities. The problem with your model is that it has a heavy discontinuity in the derivative at t=0. (A discontinuity is normal, but none of this kind). This causes problems with the adaptive step size control at the very beginning of the integration. Such a discontinuity suggests that there is something wrong with either your model or your initial past. The latter doesn’t matter if you only care about long-term behaviour, but this doesn’t seem to be the case here. I added a passage to the documentation about this kind of issue.

Wrzlprmft
  • 4,234
  • 1
  • 28
  • 54
  • Well, it doesnt work to me, it looks like the correct shape but not the values, since the one I've posted was from t=0 to 270.Even so, keep trow me the following warning: The target time is smaller than the current time. No integration step will happen. The returned state will be extrapolated from the interpolating Hermite polynomial for the last integration step. You may see this because you try to integrate backwards in time, in which case you did something wrong. You may see this just because your sampling step is small, in which case there is no need to worry. – DiegoT Dec 14 '20 at 17:59
  • @DiegoT: As I said, there is probably something wrong with your handling of `T` (which appears in your code but not the equations). As long as you do not provide full information, it’s impossible to tell you where the problem arises. — The warning you cite is (hopefully) very explicit: There is nothing to worry about if your sampling step is small, which is the case here (to obtain a detailed plot). – Wrzlprmft Dec 14 '20 at 18:52
  • I integrated your problem with [`solve_dde`](https://github.com/jrmejansen/scipy/tree/ddeSolver) and the solution matches that of JiTCDDE. So the problem lies almost certainly with your equations or parameters. – Wrzlprmft Dec 14 '20 at 22:08
  • Thank you for your help – DiegoT Dec 16 '20 at 13:25
  • I've solved the problem, the 1/T goes on the equations or on the parameters, nos both.I ask a more question: – DiegoT Apr 27 '21 at 18:10