Comments (9)
First, some clarification: f_sym
is not called; it is translated to C code, which is then called. It is a list of (or generator function yielding) SymEngine expressions. You cannot replace these expressions with functions or something similar as SymEngine won’t know what to do with them.
-
I cannot run your example code because several things are undefined. If you want to define something piecewise, you have to build it with a sharp sigmoidal for two reasons:
- SymEngine doesn’t support an analogue to SymPy’s Piecewise.
- Integrators do not like discontinuities.
-
Creating multiple
jitcdde
objects hardly ever makes sense. To access the current time, just uset
(imported withfrom jitcdde import t
) like you do in a normal delay term. -
Integrating over the past is not natively supported (see Issue #2). What you can try is to implement a reasonable discretisation of your integral into
f_sym
. -
There is an example of a state-dependent delay here. Time-dependent delays are not that different. That being said, if you can provide a good example (in form of equations), ideally with some expectations about the outcome, I might implement it as an example.
from jitcdde.
Thank you for your reply!
I know f_sym can't be called and it needs a SymEngine expressions. Maybe I should use an example to explain.
If I want to implement following piecewise function (just an example, not my real model),
As you see, I want the expression to be changed with current time t, how can I implement it? I tried following code, but it doesn't work.
import numpy as np
from jitcdde import jitcdde, y, t
t1=5.0
t2=10.0
def f(): # I just use this function to define the delay.
if t<t1: # Use t or DDE.t?
return y(t)
elif t>t1 and t<t2:
return y(t-t1)
else:
return (t-t2)*y(t-t1)+y(t-t2)
fdelay = [f()]
DDE = jitcdde(fdelay)
DDE.constant_past([1.0])
DDE.step_on_discontinuities()
data = []
tlist = []
for time in np.arange(DDE.t, DDE.t+1000, 10):
data.append( DDE.integrate(time)[0] )
tlist.append(time)
As you said in your last reply, I use t to get current time. However, I don't understand how can translated C codes know the changed expression of delay term if I don't update list fdelay.
1.I just put a part of my previous code so it indeed lost some definition. Just want to illustrate my question.
2.As you see above, I doubted that whether jitcdde could get the changed expression from a single list fdelay, so I try to create multiple jitcdde objects in my previous codes.
3.This Issue #2 is similar to my real model (), and I have done discretisation on my integral using scipy.integrate.simps. Although it seems awkward, it does give right symbolic expression. Thank you for your explanation.
4.I know this state-dependent example, but it just add a y as a augment of delay function. Sometimes we need to change the whole function form of delay, as I show in above example. Besides, I think an integral delay like my real model and Issue #2 could be common, maybe they are good examples.
5.How can I add some extra parameters? I may have to sweep those parameters to get multiple solutions.
Thank you again!
Wance
from jitcdde.
-
The problem with your code is that the conditional is evaluated in the line
fdelay = [f()]
, i.e., before you even calljitcdde
. The best way to implement your condition are probably sharp sigmoidals. Also note that you have to give an argument to indicate which component of y you want to access – even if it is one-dimensional.import numpy as np from jitcdde import jitcdde, y, t from symengine import tanh, Rational eps = Rational(1,100) sigmoid = lambda x: (tanh(x/eps)+1)/2 t1 = 5.0 t2 = 10.0 f = [ sigmoid(-t+t1)*y(0,t) + sigmoid( t-t1)*sigmoid(-t+t2)*y(0,t-t1) + sigmoid( t-t2)*( (t-t2)*y(0,t-t1)+y(0,t-t2) ) ] DDE = jitcdde(f,verbose=False) DDE.constant_past([1.0]) DDE.step_on_discontinuities() times = np.arange(DDE.t,DDE.t+20,1) data = [ DDE.integrate(time)[0] for time in times ]
-
An alternative to the above is to integrate up to
t1
, useget_state
to obtain the past, initialise a new integrator with it, integrate up tot2
and so on. However, after every re-intialisation, you would have to address initial discontinuities again (unless the cases are joint in a sufficiently smooth fashion). This approach may be computationally more feasible if you dwell in one case long enough, as you can avoid computations of complicated and irrelevant zeros (or almost zeros) likesigmoid(-100)
. -
simps
seems to handle symbolic input well, if you operate it correctly (e.g.,simps( [ y(0,t-τ) for τ in range(-10,0) ])
). The only drawback is that it uses floats for factors, which prevents certain optimisations (which should not matter much though). Anyway, right now, the best I can do is to supply a utility function that is similar tosimps
. The only better alternative is to implement a semi-symbolic integration on a low level, which is a lot of work. -
First note that neither of your examples so far requires time-dependent delays. Still, if you want such a thing, just write a respective term like you would on paper, e.g.,
y(0,t-2-sin(t))
would give you a time-dependent delay. -
Use the argument
control_pars
ofjitcdde
(look at its documentation for details).
from jitcdde.
I forgot to reply since I am busy recently. Thank you for your solution!
Your "sharp sigmoidals" solution is quite clever! But it's impractical for my case: the delay in my model is (I didn't show all details in the last comment.)
When t>delta, the delay is easy to handle. The question is when 0<t<delta, it requires a dynamic discretization of the integral, which can't be described by sharp sigmoidals. And I also want to make the discretization more precise by adding more points at first, because the integral region is small when t is small.
I did think up a solution: I added a list of anchors whose value is 0 when t<0, and one anchor valued 1 at t=0.
anchors = []
for tim in np.linspace(-delta+tstart,tstart,n+1):
anchors.append((tim,[0.0,0.0],[0.0,0.0]))
anchors[-1] = (tstart,[re0,im0],[func.we(0.0)*im0,-func.we(0.0)*re0])
Then I only use this delay form,
But in this way I introduce a discontinuity at t=0, the solution will lose a part after I step_on_discontinuty. What's worse, the output result doesn't show what the approximated theory expected. I doubt if the discontinuity cause some mistakes.
So I still hope to find a dynamic discretation which can be accepted by jitcdde. I will try get_state() when I have time. Maybe it will work. Or I will try to use other packages to solve 0<t<delta first, then input it into jitcdde.
Thank you!
from jitcdde.
Before I answer this:
This equation contains τ outside the integral as well as a the integration variable. This is at the very least bad notation, but it also suggests that your actual conditions are something like 0<t<δ and t>δ. This would make the entire problem much easier to handle.
from jitcdde.
Oh, you're right. The "tau" outside the integral should be "delta". I made mistakes when I was writing the comment.
from jitcdde.
Okay, here are some ways to do it:
Preamble
from symengine import cos, tanh, Symbol, sympify, Rational
from jitcdde import t, y, jitcdde, past_y, anchors
import numpy as np
n = 1
def trapezoid(integrand,variable,lower,upper,steps=21):
step_size = (sympify(upper-lower)/sympify(steps-1)).simplify()
def evaluations():
yield integrand.subs(variable,lower)/2
for i in range(1,steps-1):
position = lower + i*step_size
yield integrand.subs(variable,position)
yield integrand.subs(variable,upper)/2
result = sum(x for x in evaluations())*step_size
result = result.subs({
past_y(t,i,anchors(t)): y(i)
for i in range(n)
})
return result
eps = Rational(1,100)
sigmoid = lambda x: (tanh(x/eps)+1)/2
τ = Symbol("τ")
T = cos(τ)
δ = 2
The main part of this is a symbolic implementation of the trapezoidal method. The complex substitution at the end is to avoid unnecessary overlap. (I will implement a similar function soon.)
The rest are imports, definitions, and an example problem.
Variant 1: One Integral
This is based on the insight that the first case of your integral definition is not necessary if T(t,τ) = 0 for τ<0. Thus we only need apply the integral once and can avoid time-dependent delays altogether. We need a sigmoid to ensure T(t,τ) = 0 for τ<0 though:
f = [ trapezoid( sigmoid(τ)*T*y(0,τ), τ, t-δ, t ) ]
DDE = jitcdde(f)
DDE.compile_C(simplify=False)
DDE.constant_past([1.0])
for time in np.arange(DDE.t,20,0.05):
print(time,DDE.integrate(time)[0])
Note that we do not need to address initial discontinuities as there are none per construction: f(t=0)=0 since the integral is over a vanishing interval. The line DDE.compile_C(simplify=False)
exists to avoid a lengthy and pointless simplification attempt of the integral.
This variant means that the integral is less precise at the beginning. This in turn only is a problem if the integral is large with respect to rest of the derivative and the initial value. The former is the case here; the latter isn’t.
Variant 2: Swapping Integrals with a Sigmoid
Here, we use an adaptive integral for the beginning and swap it with a fixed one (like in the previous variant) via sigmoidals later.
f = [
sigmoid(-t+δ)*trapezoid(T*y(0,τ),τ, 0 ,t)
+ sigmoid( t-δ)*trapezoid(T*y(0,τ),τ,t-δ,t)
]
DDE = jitcdde(f,max_delay=1e8)
DDE.compile_C(simplify=False)
DDE.constant_past([1.0])
for time in np.arange(DDE.t,20,0.05):
print(time,DDE.integrate(time)[0])
The problem with this is that we have to calculate the first integral all the time, even if it contributes (almost) nothing. Since the first integral always starts at 0, this also means that we have to drag a lot of the past data around.
Variant 3: Reinitialising the Integrator
This is like the above, just that we do not use sigmoids to switch, but instead reinitialise the integrator with a different f:
f_1 = [trapezoid(T*y(0,τ),τ, 0 ,t)]
f_2 = [trapezoid(T*y(0,τ),τ,t-δ,t)]
DDE_1 = jitcdde(f_1,max_delay=δ,verbose=False)
DDE_1.compile_C(simplify=False)
DDE_1.constant_past([1.0])
for time in np.linspace(0.05,δ,21):
print(time,DDE_1.integrate_blindly(time)[0])
DDE_2 = jitcdde(f_2,verbose=False)
DDE_2.compile_C(simplify=False)
DDE_2.add_past_points(DDE_1.get_state())
for time in np.arange(DDE_2.t,20,0.05):
print(time,DDE_2.integrate(time)[0])
This is arguably the cleanest variant and comes with the fewest problems. Note that we need to integrate blindly in the first part to make sure we finish the last integration step precisely at δ.
Note that all three variants produce comparable results for this example.
from jitcdde.
I now implemented a utility function for integration called quadrature
– which is similar to trapezoid
from the above examples. You can now do something like:
from symengine import cos, Symbol
from jitcdde import t, y, jitcdde, quadrature
import numpy as np
δ = 2
τ = Symbol("τ")
f_1 = [quadrature(cos(τ)*y(0,τ),τ, 0 ,t)]
f_2 = [quadrature(cos(τ)*y(0,τ),τ,t-δ,t)]
DDE_1 = jitcdde(f_1,max_delay=δ,verbose=False)
DDE_1.compile_C(simplify=False)
DDE_1.constant_past([1.0])
for time in np.linspace(0.05,δ,21):
print(time,DDE_1.integrate_blindly(time)[0])
DDE_2 = jitcdde(f_2,verbose=False)
DDE_2.compile_C(simplify=False)
DDE_2.add_past_points(DDE_1.get_state())
for time in np.arange(DDE_2.t,20,0.05):
print(time,DDE_2.integrate(time)[0])
I would be interested in any feedback you have on this.
from jitcdde.
Hi,
I am sorry to not reply you for a long time, and really thanks for your solutions! It must cost you a lot of time.
In fact, I have had tried your Variant 1 last month, and meanwhile, I also implemented my own solution (First compute the region (0, δ) using scipy.integrate.solve_ivp, then import the data as anchors into jitcdde. In this way I don't have to change the expression of integrand.) Except for some regions that jitcdde step on discontinuity, the results are perfectly matched (although my method is quite slow.)
Your Variant 1:
My solution:
I haven't tried your quadrature, now I am still using a modified variant 1. If I have time (probably after this semester), I will try it and give you feedback. Thank you again!
Best regards!
from jitcdde.
Related Issues (20)
- I need your help to solve my delay differential equations HOT 35
- How do I use non-symbolic functions? HOT 1
- Issue with C++ distribution HOT 2
- How to pass variable parameter to function HOT 1
- Weird MacOS vs. Linux / Ubuntu behaviour HOT 8
- Invalid pointer when importing matplotlib HOT 3
- How to speed up integration? HOT 2
- euler maruyama scheme for stochastic delay differential equation HOT 3
- How to speed up setting parameters? HOT 3
- Cryptic error and traceback when compiling model with control_pars HOT 7
- [Docu] Kuramoto example should include `delays` and `max_delay` HOT 1
- Setting delays results in weird crash HOT 2
- Temp files sometimes aren't deleted HOT 11
- The Lyapunov exponent spectrum is not smooth HOT 14
- How to simulate the famous Lang-Kobayashi semiconductor laser model with delayed optical feedback HOT 1
- How to get the first 30 Lyapunov exponents of this modle? HOT 6
- Intended way to make time-dependent parameters HOT 5
- Solving delay differential equation with variable delay using Jitcdde HOT 1
- clang: error: the clang compiler does not support '-march=native' HOT 1
- warning: pointer is missing a nullability type specifier HOT 6
Recommend Projects
-
React
A declarative, efficient, and flexible JavaScript library for building user interfaces.
-
Vue.js
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
-
Typescript
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
-
TensorFlow
An Open Source Machine Learning Framework for Everyone
-
Django
The Web framework for perfectionists with deadlines.
-
Laravel
A PHP framework for web artisans
-
D3
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
-
Recommend Topics
-
javascript
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
-
web
Some thing interesting about web. New door for the world.
-
server
A server is a program made to process requests and deliver data to clients.
-
Machine learning
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
-
Visualization
Some thing interesting about visualization, use data art
-
Game
Some thing interesting about game, make everyone happy.
Recommend Org
-
Facebook
We are working to build community through open source technology. NB: members must have two-factor auth.
-
Microsoft
Open source projects and samples from Microsoft.
-
Google
Google ❤️ Open Source for everyone.
-
Alibaba
Alibaba Open Source for everyone
-
D3
Data-Driven Documents codes.
-
Tencent
China tencent open source team.
from jitcdde.