Time dependent constant in the equation

[rongboyang]

Is there a way to add a time-dependent constant in the model? Let's say for some reason I want to add a function of time t as a constant in the mackey_glass_jump example, notice that I added sin(t) at the end of the function, but that would not work. What is the correct syntax for it?

def time_dependent_constant(x):
    return 1.2 * x

τ = 15
n = 10
β = 0.25
γ = 0.1

f = [ β * y(0,t-τ) / (1 + y(0,t-τ)**n) - γ*y(0) + time_dependent_constant(t) ]
DDE = jitcdde(f,verbose=False)

DDE.constant_past([1.0])

DDE.step_on_discontinuities()

for time in numpy.arange(DDE.t, DDE.t+1000, 1):
    print( time, *DDE.integrate(time) )

DDE.jump(1,DDE.t)

for time in numpy.arange(DDE.t, DDE.t+1000, 1):
    print( time, *DDE.integrate(time) )

[Wrzlprmft]

After imports, your example works, so I can only guess what the problem is:

As the input for JiTCDDE is symbolic, you need to use symbolic functions, i.e., symengine.sin(t). Using numerical implementations of the sine won’t work, since you cannot compute the numerical sine of a symbol (e.g., numpy.sin(t)). Please also see this.

[rongboyang]

I see. What if I have a time-dependent constant which is a big array, something like [200, 320, 130, ...., 230], it's the result of some complex computation. What do you suggest in order to make it symbolic, is there a way? I am trying to utilize this DDE solver for an engineering problem, so my question may seem a little bit different.

result = [200, 320,  130, ....,  230]
def time_dependent_constant(t):
    return result[t]
...

data = []
for time in numpy.arange(DDE.t, DDE.t+1000, 1):
    data.append([time, DDE.integrate(time)])

[Wrzlprmft]

I see. To do something like this, virtually move your input data to the past of a dummy component. Then create an initial past that features your result. Please see this issue for JiTCODE for details (second bullet point, and maybe some of the follow-up questions).

[rongboyang]

Hi Wrzlprmft,

Thanks for your previous reply!

Right now what I need to solve is a non-autonomous DDE, in other words, the right-hand side of the equation contains a function that explicitly depends on the time t. The function is numerically defined by an array. I couldn't find an example in the sample code. Do you mind providing me a way to solve it? Thank you very much!

[Wrzlprmft]

Do you mind providing me a way to solve it?

Okay, here is an example code:

import numpy as np
from jitcdde import y,t,jitcdde
from chspy import CubicHermiteSpline
from matplotlib import pyplot as plt

# Replace with something reasonable:
input_times = np.linspace(0,70,20)
input_data = np.random.normal(size=20)
# The input_times must be sorted:
assert np.all(np.diff(input_times)>0)

n = 2 # number of dynamic variables of your DDE

# Function to use the input stored in an extra dummy dynamic variable (the n-th one):
max_time = input_times[-1]
my_input = lambda t: y(n,t-max_time)

τ = 10
f = [
        -y(0) + y(1,t-τ) + my_input(t),
        -y(1) + y(0,t-τ),
    ]
assert len(f)==n

# Add dynamics for dummy dynamic variable (do nothing).
f.append(0)

past = CubicHermiteSpline(n=n+1)
for time,datum in zip(input_times,input_data):
    # Here I assume that your initial past for the DDE is all zero:
    regular_past = np.zeros(n)

    past.append((
            time-max_time,
            np.hstack((regular_past,datum)),
            np.zeros(n+1), # zero derivative at the anchors
        ))

times = np.linspace(-max_time,0,1000)
plt.plot(times+max_time,past.get_state(times)[:,n])
plt.title("This is your input as the integrator will see it:")
plt.show()

DDE = jitcdde(f,verbose=False)
DDE.add_past_points(past)
DDE.adjust_diff()

times = np.linspace(DDE.t+1,max_time,100)
result = np.vstack([DDE.integrate(time)[:n] for time in times])
plt.plot(times,result)
plt.show()

However, I strongly urge you to understand every step of this instead of blindly adapting this. You need to apply your knowledge of your system and your requirements.

In particular, note all the saddles in the implementation of the past. The reason for these is that I used a zero derivative at all the anchors – which is where the saddles are. This can usually be done better if you have a somewhat continuous input (instead of random numbers as in the example) and can apply your knowledge about it.

[rongboyang]

This just worked! Words can't express my gratefulness!!!

[rongboyang]

Hi Wrzlprmft,

Sorry for bringing this thread up again. I was wondering what should I do if I want to add multiple constants in the equation. I tried something based on your suggestion but it's not working properly. Do you mind taking a loot at it?

import numpy as np
from jitcdde import y, t, jitcdde
from chspy import CubicHermiteSpline
from matplotlib import pyplot as plt

# Replace with something reasonable:
input_times = np.linspace(0, 70, 20)

input_data = [0.42935087, 1.92088671,  0.41712295, -0.25007707, -0.18954231, -1.24214212, -1.07214407, -0.91995987,
              0.51654649, -1.25133288,  0.29660553, -0.30552719, 0.29201498, 1.10732322, -0.18839248, -0.99940753,
              -1.95084375, 0.91911038, -1.14697187, 0.86489447]

input_data2 = [0.42935087, 1.92088671,  0.41712295, -0.25007707, -0.18954231, -1.24214212, -1.07214407, -0.91995987,
              0.51654649, -1.25133288,  0.29660553, -0.30552719, 0.29201498, 1.10732322, -0.18839248, -0.99940753,
              -1.95084375, 0.91911038, -1.14697187, 0.86489447]

regular_past_values = [[5, 2], [4, 2], [5, 2], [5, 3], [1, 2], [5, 2], [4, 2], [5, 2], [5, 3], [1, 2],
                [5, 2], [4, 2], [5, 2], [5, 3], [1, 2],[5, 2], [4, 2], [5, 2], [5, 3], [10, 3]]

# The input_times must be sorted:
assert np.all(np.diff(input_times) > 0)

n = 2  # number of dynamic variables of your DDE

# Function to use the input stored in an extra dummy dynamic variable (the n-th one):
max_time = input_times[-1]
start_time = 0
my_input = lambda t: y(n, (t + start_time) - max_time)
my_input_2 = lambda t: y(n, (t + start_time) - max_time)
τ = 10
f = [
    -y(0) + y(1, t - τ) + my_input(t),
    -y(1) + y(0, t - τ) + my_input_2(t)
]

assert len(f) == n
f.append(0)
f.append(0)
past = CubicHermiteSpline(n=n+2)

for time, datum, datum_2, regular_past in zip(input_times, input_data, input_data2, regular_past_values):

    past.append((
        (time + start_time) - max_time,
        np.hstack((regular_past, datum, datum_2)),
        np.zeros(n + 2),  # zero derivative at the anchors
    ))

DDE = jitcdde(f, verbose=False)
DDE.add_past_points(past)
DDE.adjust_diff()

times = np.linspace(DDE.t, max_time, 100)
result = np.vstack([DDE.integrate(time)[:n] for time in times])
# plt.plot(input_times, input_data)
plt.plot(times, result)
plt.show()

[Wrzlprmft]

I could spot an error in this line:

 my_input_2 = lambda t: y(n, (t + start_time) - max_time)

The number n is the number of the (dummy) dynamic variable where your input is stored. n is the first beyond the regular dynamic variables (due to zero-based counting). What you want here is the second, which is n+1:

 my_input_2 = lambda t: y(n+1, (t + start_time) - max_time)

In general, this can look like this (skipping some lines):

m = 2 # number of inputs

my_input = [
        lambda t: y(n+i, (t + start_time) - max_time)
        for i in range(m)
    ]

f = [
    -y(0) + y(1, t - τ) + my_input[0](t),
    -y(1) + y(0, t - τ) + my_input[1](t)
]

There may be further issues I missed, but without any specification what kind of output you expect or what goes wrong, it’s difficult to spot them.

[Wrzlprmft]

I now added an interface that comprises the above and allows you to add input to JiTCDDE without fiddling with indices. In addition, I added an class method to CHSPy to create better splines from data. Finally, CHSPy now has a plotting method.

It would be great if you could test it.

Here is an example very similar to the above using all these new features:

import numpy as np
from jitcdde import y, t, jitcdde_input, input
from chspy import CubicHermiteSpline
from matplotlib import pyplot as plt

# Defining Input
# --------------

input_times = np.linspace(0, 70, 20)
max_time = input_times[-1]

input_data = np.random.normal(size=(20,2))
input_spline = CubicHermiteSpline.from_data(input_times,input_data)

fig,axes = plt.subplots()
input_spline.plot(axes)
axes.set_title("input")
plt.show()

# Defining Dynamics
# -----------------

τ = 10
f = [
    -y(0) + y(1,t-τ) + input(0),
    -y(1) + y(0,t-τ) + input(1),
]

# Actual Integration
# ------------------

DDE = jitcdde_input(f,input_spline)
DDE.constant_past(np.zeros(len(f)))
DDE.adjust_diff()
times = np.linspace(DDE.t, max_time, 100)
result = np.vstack([DDE.integrate(time) for time in times])

# Plotting
# --------

fig,axes = plt.subplots()
axes.plot(times, result)
axes.set_title("result")
plt.show()

[Wrzlprmft]

A few days ago, I added a feature to use callbacks, which is an alternative way to implement your input. However, you have to take care of he continuity and it is probably slower.