New feature to access past values of function derivative

neurophysik / jitcdde

Just-in-time compilation for delay differential equations

Other

57 stars 15 forks source link

New feature to access past values of function derivative #24

Closed nchasovnikov closed 5 years ago

nchasovnikov commented 5 years ago

I would like to use JiTCDDE to solve neutral DDE system of equations. This system represents a physical model of sound production in a recorder (musical instrument). I would like to reproduce results presented in an article by Auvray et al. https://asa.scitation.org/doi/10.1121/1.3672815 The system of equations looks like this:

$y(t) = \sum{y_n(t)}\\ \ddot{y}_n + \varepsilon_n\nu_n \dot{y}_n + \nu^2_n y_n =\\ =\mu_n \frac{d^2}{dt^2}tanh(y(t-\tau)-y_0) - \zeta_n \frac{d}{dt} |y|^2$

And is a set of damped harmonic oscillator coupled by non-linear functions.

It would be great to be able to access past values of function derivative although it would lead to necessity of designing initial past values.

Wrzlprmft commented 5 years ago

I will try to tend to this in the next week.

Until then, can you give me (as far as possible):

Some values for the control parameters.
The explicit first-order dynamical equations.
Some initial conditions.
Some idea what the solution looks like.

… so I don’t have to extract this from the paper.

nchasovnikov commented 5 years ago

Solution should normally be wave like with frequencies quite close to respective nu_n. It was mentioned in the article they used arbitrary values for history.

Control parameters: $\varepsilon_n = [0.03966, 0.03184, 0.02847]\\ \nu_n = [1.000, 2.033, 3.066]\\ \mu_n =[0.16115667456085775, 0.14093420256851111, 0.11465065353644151]\\ \bar{y}_0 = 0.0\\ \tau = 1.7735\\ \zeta_n = [0.017940997406325931, 0.015689701773967984, 0.012763648066925721]$

System of equations: $\begin{array}{lcl} \ddot{y_0} &=& \frac{d}{dt}\dot{y_0}\\ \ddot{y_1} &=& \frac{d}{dt}\dot{y_1}\\ \ddot{y_2} &=& \frac{d}{dt}\dot{y_2}\\ \ddot{y_0} &=& 2 \mu_0 sech(y(t-\tau)-\bar{y}_0)^2 \cdot\\ && (\ddot{y}(t-\tau)-2 \dot{y}(t-\tau)^2tanh(y(t-\tau)-\bar{y}_0)) \\ && -2 \zeta_0 |y(t)| \dot{y}(t) -\varepsilon_0\nu_0 \dot{y}_0 - \nu^2_0 y_0 \\ \ddot{y_1} &=& 2 \mu_1 sech(y(t-\tau)-\bar{y}_0)^2 \cdot\\ && (\ddot{y}(t-\tau)-2 \dot{y}(t-\tau)^2tanh(y(t-\tau)-\bar{y}_0)) \\ && -2 \zeta_1 |y(t)| \dot{y}(t) -\varepsilon_1\nu_1 \dot{y}_1 - \nu^2_1 y_1 \\ \ddot{y_2} &=& 2 \mu_2 sech(y(t-\tau)-\bar{y}_0)^2 \cdot\\ && (\ddot{y}(t-\tau)-2 \dot{y}(t-\tau)^2tanh(y(t-\tau)-\bar{y}_0)) \\ && -2 \zeta_2 |y(t)| \dot{y}(t) -\varepsilon_2\nu_2 \dot{y}_2 - \nu^2_2 y_2 \\ \end{array}$ Here: $y(t - \tau) = y_0(t-\tau) + y_1(t-\tau) + y_2(t-\tau)\\ \dot{y}(t - \tau) = \dot{y}_0(t-\tau) + \dot{y}_1(t-\tau) + \dot{y}_2(t-\tau)\\ \ddot{y}(t - \tau) = \ddot{y}_0(t-\tau) + \ddot{y}_1(t-\tau) + \ddot{y}_2(t-\tau)\\ y(t) = y_0(t) + y_1(t) + y_2(t)\\ \dot{y}(t) = \dot{y}_0(t) + \dot{y}_1(t) + \dot{y}_2(t)$

Using six functions instead of three allows me to access past and present values of first derivative. But in the first term of equations 4-6 of a system I need past value of second derivative (in another words derivative of functions 1-3 of a system).

Wrzlprmft commented 5 years ago

I finally worked on this. You can find everything in this new branch (and this is what you have to check out to run the code). In particular, here is an implementation of your code.

The problem is that I get escalating oscillations with this, no matter what I do. I therefore suspect that the fact that the DDE is neutral is not the problem here. Can you please double-check the above equations as well as my implementation?

nchasovnikov commented 5 years ago

Thank you. I will look them through. I've managed to solve this problem (system of equations) using PyDDE package, but performance of a program is not what I'm aiming for. I found that escalation begins somewhere before first period of lag time elapses. Will try to look deeper into "Dissipative Approximations" method by Shampine.

PS: I've managed to find some typos in system of equations I provided, but even with them fixed problem is not solved. Providing harmonic past function (using frequencies nu_n) doesn't solve it too.

Wrzlprmft commented 5 years ago

I've managed to find some typos in system of equations I provided, but even with them fixed problem is not solved.

Can you please share them, so I can perhaps get an idea where the problem is?

nchasovnikov commented 5 years ago

Functions should look like these: f = {y(i): ydot[i] for i in range(3)} f.update({ydot[i]: μ[i] * sech(ybar_0 - y_tot(t-τ))**2 *\ (ydotdot_tot(t-τ) + 2 * ydot_tot(t-τ)**2*tanh(ybar_0 - y_tot(t-τ)))\ - 2*ζ[i] * abs(y_tot(t)) * ydot_tot(t)\ - ε[i] * ν[i] * ydot[i]\ - ν[i]**2 * y(i) for i in range(3)})

Is it a good idea to use a built-in Abs function from Symengine?

Wrzlprmft commented 5 years ago

Functions should look like these: f = {y(i): ydot[i] for i in range(3)} f.update({ydot[i]: μ[i] * sech(ybar_0 - y_tot(t-τ))**2 *\ (ydotdot_tot(t-τ) + 2 * ydot_tot(t-τ)**2*tanh(ybar_0 - y_tot(t-τ)))\ - 2*ζ[i] * abs(y_tot(t)) * ydot_tot(t)\ - ε[i] * ν[i] * ydot[i]\ - ν[i]**2 * y(i) for i in range(3)})

I still get escalations with these. However, everything is perfectly smooth as intended, so I am still suspecting an error in the equations rather than with the fact that the DDE is neutral. Can you show me the working PyDDE code, so I can compare myself?

Is it a good idea to use a built-in Abs function from Symengine?

Non-smooth functions are generally a problem when integrating; hence I routinely replace them with smooth approximations.

nchasovnikov commented 5 years ago

FOUND!

lambda for sech is incorrect. should be: sech = lambda x: 2 / (exp(x)+exp(-x))

nchasovnikov commented 5 years ago

Compared results by JITCDDE and pyDDE - nearly equivalent.

Thank you, Garrit!

Wrzlprmft commented 5 years ago

lambda for sech is incorrect. should be: sech = lambda x: 2 / (exp(x)+exp(-x))

Damnit. Mea culpa.

Now that it works, I optimised your example to make maximal use of JiTCDDE’s speed-optimisation. Your case provides plenty of opportunity for this; I could halve the integration time. I found and fixed a bug in the process (so you need to update before you can use this).

Please report any problems you encounter while working with this. As of now certain features like Lyapunov exponents may not work with neutral DDEs. I will try to address these soon and craft everything into a release.

Compared results by JITCDDE and pyDDE - nearly equivalent.

Great to know. That saves me a considerable amount of testing work.