Time dependent constant in function

Hello, if shortly the problem is the following - I am trying to implement a mathematical model of induction motor [https://www.mathworks.com/matlabcentral/fileexchange/63759-dynamic-model-of-3-ph-im-in-synchronously-rotating-reference-frame?s_tid=prof_contriblnk]. Based on previous discussions in the topic "problems about time-dependent delay and jitcdde parameter f_sym", I build a piecewise implementation (shown below), but unfortunately it was not functional. You can explain what the problem is, I think it would be useful to many.

from numpy import sin,cos,sqrt,dot,array,angle,linspace
from jitcdde import y,t,jitcdde
import matplotlib.pyplot as plt
import math
plt.rcParams['figure.figsize'] = (10.0, 4.0)
#===============================================================
#Voltage Source
#===============================================================
def volt_source(f,v_rms,t_i):
    omega=2*math.pi*f
    v_r=v_rms*sin(omega*t_i)
    v_y=v_rms*sin(omega*t_i-2*math.pi/3)
    v_b=v_rms*sin(omega*t_i+2*math.pi/3)
    A=array([[1,-0.5,-0.5],[0,math.sqrt(3)/2,-math.sqrt(3)/2]])
    [v_a,v_b]=dot(A,array([v_r,v_y,v_b]))*sqrt(2/3)
    theta=angle(complex(v_a,v_b))
    v_d=v_a*cos(theta)+v_b*sin(theta)
    v_q=-v_a*sin(theta)+v_b*cos(theta)
    return [v_d,v_q,theta]
#===============================================================
#Parameters
#===============================================================
f=50
v_rms=230
R_s=0.435 #Stator resistance
R_r=0.816 #Rotor resistance
L_ls=0.002 #Stator inductance
L_lr=0.002 #Rotor inductance
L_m=0.0693
M=L_m # Mutual inductance
P=2 #Number of Poles
J=0.089 #Inertia
B=0.01 #friction coefficient
T_l=0.0 #Load torque
W_e=2*math.pi*f
L_ss=L_ls+L_m
L_rr=L_lr+L_m
v_d=[]
v_q=[]
#===============================================================
#Initial conditions
#===============================================================
delta_t=1e-4
initial_state=array([.0,.0,.0,.0,.0])
n=1
num_sampl=10000
times=linspace(0,n,n*num_sampl)
for time in times:
    v_d.append(volt_source(f,v_rms,time)[0])
    v_q.append(volt_source(f,v_rms,time)[1])
#===============================================================
# Solution (piecewise)
#===============================================================
result=[]
for i,valu_time in enumerate(times):
   #==============================================================
    #Equations dynamic model induction motor with synchronously rotating reference frame
  #===============================================================
    v_q_t=v_q[i]
    v_d_t=v_d[i]
    equat_dyna_model=[
            (-((y(3)-y(3,t-delta_t))/delta_t)*L_m+v_q_t-W_e*(y(1)*
               (L_ls+L_m)+y(2)*L_m)-R_s*y(0))/(L_ls+L_m),
            (-((y(2)-y(2,t-delta_t))/delta_t)*L_m+v_d_t+W_e*(y(0)*
               (L_ls+L_m)+y(3)*L_m)-R_s*y(1))/(L_ls+L_m),
            (-((y(1)-y(1,t-delta_t))/delta_t)*L_m+(W_e-y(4)*2/P)*
             (y(3)*(L_lr+L_m)+y(0)*L_m)-R_r*y(2))/(L_lr+L_m),
            (-((y(0)-y(0,t-delta_t))/delta_t)*L_m-(W_e-y(4)*2/P)*
             (y(2)*(L_lr+L_m)+y(1)*L_m)-R_r*y(3))/(L_lr+L_m),
            ((((y(0)*(y(1)+y(2))*L_m+y(2)*L_lr-y(1)*(y(0)+y(3))*
                L_m+y(3)*L_lr)*(L_m/(L_lr+L_m))*(2/3*P/2))-T_l)*(1/J)-y(4)*B/J)
            ]
    solut=jitcdde(equat_dyna_model)
    solut.constant_past(initial_state,time=valu_time)
    result.append(solut.integrate(valu_time))
    initial_state=result[-1]
    time=valu_time

The above code runs for me, but it probably does not do what you want. You have the initialisation of the integration within the time loop, which means that you intialise the integrator again and again. For most applications, you want a rough layout like this:


result = []
equat_dyna_model = […]
solut = jitcdde(equat_dyna_model)
solut.constant_past(initial_state)
solut.step_on_discontinuities()
for i,valu_time in enumerate(times):
    result.append(solut.integrate(valu_time))

Also, importing sin and similar from numpy is usually a bad sign. Please read this.

Finally, your link tells me little about what your model is.

Before I delve deeper into this, can you please check the above? Also, please state what you expect and how the result differs? “unfortunately it was not functional” can mean a lot of things.

Thanks for the feedback) And so let's start: 1) I got acquainted with the materials by reference (and got acquainted with the article on this library ) - as a result, I rewrote the implementation (induction motor model - IMM) of the function def volt_source(f,v_rms,t_i) on sympy and symengine (the code is below)

##=============================================================================
##SAGAR KUMAR DASH M.TECH in Power Electronics;
##INDIAN INSTITUTE OF TECHNOLOGY, TRIVANDRUM
##Dynamic Model of I.M. in Synchronously rotating Reference Frame
##=============================================================================
import symengine
import sympy
from numpy import linspace,array
from jitcdde import y,t,jitcdde
import matplotlib.pyplot as plt
import math
plt.rcParams['figure.figsize'] = (10.0, 4.0)
#==============================================================================
#Voltage Source
#==============================================================================
def volt_source(f,v_rms,t_i):
    omega=2*symengine.pi*f
    v_r=v_rms*symengine.sin(omega*t_i)
    v_y=v_rms*symengine.sin(omega*t_i-2*symengine.pi/3)
    v_b=v_rms*symengine.sin(omega*t_i+2*symengine.pi/3)
    A=symengine.DenseMatrix([[1,-0.5,-0.5],[0,symengine.sqrt(3)/2,
                             -symengine.sqrt(3)/2]])
    [v_a,v_b]=(symengine.Matrix.dot(A,symengine.DenseMatrix([v_r,v_y,v_b]))*
    symengine.sqrt(2/3))
#    theta=sympy.arg((v_a+v_b*sympy.I))
    theta=sympy.log(v_a+v_b*sympy.I).as_real_imag()[1]
    v_d=v_a*symengine.cos(theta)+v_b*symengine.sin(theta)
    v_q=-v_a*symengine.sin(theta)+v_b*symengine.cos(theta)
    return [v_d,v_q,theta]
#==============================================================================
#Parameters
#==============================================================================
f=50
v_rms=230
R_s=0.435 #Stator resistance
R_r=0.816 #Rotor resistance
L_ls=0.002 #Stator inductance
L_lr=0.002 #Rotor inductance
L_m=0.0693
M=L_m # Mutual inductance
P=2 #Number of Poles
J=0.089 #Inertia
B=0.01 #friction coefficient
T_l=0.0 #Load torque
W_e=2*math.pi*f
L_ss=L_ls+L_m
L_rr=L_lr+L_m
v_d=[]
v_q=[]
#==============================================================================
#Initial conditions
#==============================================================================
delta_t=1e-3
initial_state=array([.0,.0,.0,.0,.0])
n=1
num_sampl=1000
times=linspace(0,n,n*num_sampl)
#==============================================================================
# Solution
#==============================================================================
result=[]
equat_dyna_model=[
        (-((y(3)-y(3,t-delta_t))/delta_t)*L_m+volt_source(f,v_rms,t)[1]
        -W_e*(y(1)*(L_ls+L_m)+y(2)*L_m)-R_s*y(0))/(L_ls+L_m),
        (-((y(2)-y(2,t-delta_t))/delta_t)*L_m+volt_source(f,v_rms,t)[0]
        +W_e*(y(0)*(L_ls+L_m)+y(3)*L_m)-R_s*y(1))/(L_ls+L_m),
        (-((y(1)-y(1,t-delta_t))/delta_t)*L_m+(W_e-y(4)*2/P)*
         (y(3)*(L_lr+L_m)+y(0)*L_m)-R_r*y(2))/(L_lr+L_m),
        (-((y(0)-y(0,t-delta_t))/delta_t)*L_m-(W_e-y(4)*2/P)*
         (y(2)*(L_lr+L_m)+y(1)*L_m)-R_r*y(3))/(L_lr+L_m),
        ((((y(0)*(y(1)+y(2))*L_m+y(2)*L_lr-y(1)*(y(0)+y(3))*
            L_m+y(3)*L_lr)*(L_m/(L_lr+L_m))*(2/3*P/2))-T_l)*(1/J)-y(4)*B/J)
        ]
solut = jitcdde(equat_dyna_model)
solut.constant_past(initial_state)
#solut.step_on_discontinuities()
for i,valu_time in enumerate(times):
    result.append(solut.integrate(valu_time))

2) Below I will give a description of the work that I do - yes I think your library is suitable for solving this problem default 3) new implementation is compiled, but in the end it "falls". At the moment, the main task for me is to complete this solution and debug it, in the future I plan to use this implementation with a genetic algorithm to restore parameters IMM based on experimental data

============================================================
READ ME FIRST
============================================================
Generating compiled integrator failed; resorting to lambdified functions. If you can live with using the Python backend, you can call generate_lambdas to explicitly do this and bypass the compile attempt and error messages. Otherwise, you want to take care of fixing the above errors.
============================================================

============================================================

  warn(line + "READ ME FIRST" + line + "Generating compiled integrator failed; resorting to lambdified functions. If you can live with using the Python backend, you can call generate_lambdas to explicitly do this and bypass the compile attempt and error messages. Otherwise, you want to take care of fixing the above errors." + 2*line)
Using default integration parameters.
Traceback (most recent call last):

  File "<ipython-input-4-da33e78f7cd7>", line 1, in <module>
    runfile('C:/Users/user/Desktop/6_INDUC_MOTOR/1_IMPLEMENTATION/1_PYTHON/dynamic_dq_axis_3phase_matlab_frame.py', wdir='C:/Users/user/Desktop/6_INDUC_MOTOR/1_IMPLEMENTATION/1_PYTHON')

  File "C:\Users\user\Anaconda3\lib\site-packages\spyder_kernels\customize\spydercustomize.py", line 704, in runfile
    execfile(filename, namespace)

  File "C:\Users\user\Anaconda3\lib\site-packages\spyder_kernels\customize\spydercustomize.py", line 108, in execfile
    exec(compile(f.read(), filename, 'exec'), namespace)

  File "C:/Users/user/Desktop/6_INDUC_MOTOR/1_IMPLEMENTATION/1_PYTHON/dynamic_dq_axis_3phase_matlab_frame.py", line 78, in <module>
    result.append(solut.integrate(valu_time))

  File "C:\Users\user\Anaconda3\lib\site-packages\jitcdde\_jitcdde.py", line 776, in integrate
    self.DDE.get_next_step(self.dt)

  File "C:\Users\user\Anaconda3\lib\site-packages\jitcdde\_python_core.py", line 250, in get_next_step
    k_2 = self.eval_f(self.t + 0.5 *delta_t, self.y + 0.5 *delta_t*k_1)

  File "C:\Users\user\Anaconda3\lib\site-packages\jitcdde\_python_core.py", line 240, in eval_f
    return self.f(t, y, *self.parameters)

  File "C:\Users\user\Anaconda3\lib\site-packages\jitcdde\_python_core.py", line 189, in <lambda>
    self.f = lambda *args: np.array(F(*args)).flatten()

  File "<lambdifygenerated-3>", line 2, in _lambdifygenerated
    return ([-6.10098176718091*Y[0] - 314.159265358978*Y[1] - 305.346943750031*Y[2] - 971.949509116407*Y[3] - 1316.92996923826*sqrt(3)*(sin(pi*(100*t + 1/3)) + cos(pi*(100*t + 1/6)))*cos(arg(-563.382640840131*1j*(-cos(50*pi*t) + 1)**2 + 281.691320420065*sin(100*pi*t) - 1126.76528168026*1j*cos(50*pi*t) + 845.073961260196*1j)) + 971.949509116407*past_y(t - 0.001, 3, anchors(t - 0.001)) - 3950.78990771479*sin(100*pi*t)*sin(arg(-563.382640840131*1j*(-cos(50*pi*t) + 1)**2 + 281.691320420065*sin(100*pi*t) - 1126.76528168026*1j*cos(50*pi*t) + 845.073961260196*1j)), 314.159265358978*Y[0] - 6.10098176718091*Y[1] - 971.949509116407*Y[2] + 305.346943750031*Y[3] - 1316.92996923826*sqrt(3)*(sin(pi*(100*t + 1/3)) + cos(pi*(100*t + 1/6)))*sin(arg(-563.382640840131*1j*(-cos(50*pi*t) + 1)**2 + 281.691320420065*sin(100*pi*t) - 1126.76528168026*1j*cos(50*pi*t) + 845.073961260196*1j)) + 971.949509116407*past_y(t - 0.001, 2, anchors(t - 0.001)) + 3950.78990771479*sin(100*pi*t)*cos(arg(-563.382640840131*1j*(-cos(50*pi*t) + 1)**2 + 281.691320420065*sin(100*pi*t) - 1126.76528168026*1j*cos(50*pi*t) + 845.073961260196*1j)), -971.949509116407*Y[1] - 11.4446002805049*Y[2] - 14.0252454417952*(0.0693*Y[0] + 0.0713*Y[3])*(Y[4] - 314.159265358979) + 971.949509116407*past_y(t - 0.001, 1, anchors(t - 0.001)), -971.949509116407*Y[0] - 11.4446002805049*Y[3] + 14.0252454417952*(0.0693*Y[1] + 0.0713*Y[2])*(Y[4] - 314.159265358979) + 971.949509116407*past_y(t - 0.001, 0, anchors(t - 0.001)), 0.504540082260428*Y[0]*Y[2] - 0.504540082260428*Y[1]*Y[3] + 0.0145610413350773*Y[2] + 0.0145610413350773*Y[3] - 0.112359550561798*Y[4]])

NameError: name 'arg' is not defined

In the previous version of the code, I did not want to rewrite the function def volt_source(f,v_rms,t_i) because of this, I made piecewise solutions in a loop with varying parameters for the equations, and each initialization began with the last decision default

I only have time for a quick guess now, but here we go:

JiTCDDE cannot really deal with complex numbers natively. This is because it requires some effort and never came up so far. The easiest solution is to boil down everything to pairs of real numbers. In your case, it probably suffices to calculate theta with triangular functions instead of complex numbers.

Sidenotes:

If runtime turns out to be an issue, you may want to try to streamline your calculations of the derivative using helpers (see the documentation).
in the future I plan to use this implementation with a genetic algorithm to restore parameters IMM based on experimental data

I am not exactly sure what this means, but note that changing control parameters in a running simulation may result in some difficulties: You would have to address initial discontinuities whenever you do this.
If you have any explicit time dependencies, the best way is to implement them literally. You do not need anything piecewise for this; the algorithm can deal with this. This is no problem for the algorithm and what happens right now in your code (ignoring the issue with complex numbers).

PS: Looking at your equations, you seem to only use the delay to approximate the temporal derivative of your components. I think you can avoid this by inserting your respective equations and solving them for the time derivative. (If not, this still doesn’t seem like a good use of DDEs and may cause a lot of trouble.)

For example, your equations for y₀ and y₃ are of the following forms:

ẏ₀ = −c₀ẏ₃+f₀(y), ẏ₃ = −c₃ẏ₀+f₃(y),

where c₀ and c₃ are some constants. If you insert the second equation into the first one, you get an equation that you can solve for ẏ₀. Do something similar with the other equations and you end up with a system of ordinary differential equations, which should make your life considerably easier.

And so, the model was rewritten into sympy, the original system of equations was not changed, because if you make substitutions, then the number of variables will be more than the number of equations. The resulting code works steadily, there is also the question of how to measure the time to execute the code correctly (see are set correctly times.times?).

##=============================================================================
##SAGAR KUMAR DASH M.TECH in Power Electronics
##INDIAN INSTITUTE OF TECHNOLOGY, TRIVANDRUM
##Dynamic Model of I.M. in Synchronously rotating Reference Frame
##=============================================================================
import sympy
import time
from numpy import linspace,array
from jitcdde import y,t,jitcdde
import matplotlib.pyplot as plt
plt.rcParams['figure.figsize'] = (10.0, 4.0)
#==============================================================================
#Voltage Source
#==============================================================================
def volt_source(f,v_rms,t_i):
    omega=2*sympy.pi*f
    v_r=v_rms*sympy.sin(omega*t_i)
    v_y=v_rms*sympy.sin(omega*t_i-2*sympy.pi/3)
    v_b=v_rms*sympy.sin(omega*t_i+2*sympy.pi/3)
    A=sympy.Matrix([[1,-0.5,-0.5],[0,sympy.sqrt(3)/2,
                             -sympy.sqrt(3)/2]])
    B=sympy.Matrix([v_r,v_y,v_b])*sympy.sqrt(2/3)
    [v_a,v_b]=A*B
    p=sympy.Piecewise((0,v_a>=0),(sympy.pi,v_a<0))
    theta=sympy.atan(v_b/v_a)+p
    v_d=v_a*sympy.cos(theta)+v_b*sympy.sin(theta)
    v_q=-v_a*sympy.sin(theta)+v_b*sympy.cos(theta)
    return [v_d,v_q,theta]
#==============================================================================
#Parameters
#==============================================================================
f=50
v_rms=230
R_s=0.435 #Stator resistance
R_r=0.816 #Rotor resistance
L_ls=0.002 #Stator inductance
L_lr=0.002 #Rotor inductance
L_m=0.0693
M=L_m # Mutual inductance
P=2 #Number of Poles
J=0.089 #Inertia
B=0.01 #friction coefficient
T_l=0.0 #Load torque
W_e=2*sympy.pi*f
L_ss=L_ls+L_m
L_rr=L_lr+L_m
v_d=[]
v_q=[]
#==============================================================================
#Initial conditions
#==============================================================================
delta_t=1e-4
initial_state=array([.0,.0,.0,.0,.0])
n=1
times=linspace(0,n,int(n/delta_t))
#==============================================================================
# Solution
#==============================================================================
result=[]
equat_dyna_model=[
        (-((y(3)-y(3,t-delta_t))/delta_t)*L_m+volt_source(f,v_rms,t)[1]
        -W_e*(y(1)*(L_ls+L_m)+y(2)*L_m)-R_s*y(0))/(L_ls+L_m),
        (-((y(2)-y(2,t-delta_t))/delta_t)*L_m+volt_source(f,v_rms,t)[0]
        +W_e*(y(0)*(L_ls+L_m)+y(3)*L_m)-R_s*y(1))/(L_ls+L_m),
        (-((y(1)-y(1,t-delta_t))/delta_t)*L_m+(W_e-y(4)*2/P)*
         (y(3)*(L_lr+L_m)+y(0)*L_m)-R_r*y(2))/(L_lr+L_m),
        (-((y(0)-y(0,t-delta_t))/delta_t)*L_m-(W_e-y(4)*2/P)*
         (y(2)*(L_lr+L_m)+y(1)*L_m)-R_r*y(3))/(L_lr+L_m),
        ((((y(0)*(y(1)+y(2))*L_m+y(2)*L_lr-y(1)*(y(0)+y(3))*
            L_m+y(3)*L_lr)*(L_m/(L_lr+L_m))*(2/3*P/2))-T_l)*(1/J)-y(4)*B/J)
        ]
solut=jitcdde(equat_dyna_model)
solut.constant_past(initial_state)
#solut.set_integration_parameters(1e-4,1e-2)
solut.step_on_discontinuities()
solut.integrate_blindly(delta_t)
t_solut=time.time()
for valu_time in times:
    result.append(solut.integrate(valu_time))
elapsed=time.time()-t_solut

Compare Matlab and Python implementation: 1) Omega (rotation speed) default 2)Iqs,Ids - currents on stator in dq axis - But at the initial moment of time there are some differences in the amplitude of the current, but closer to the middle the differences disappear, this requires checking default

because if you make substitutions, then the number of variables will be more than the number of equations.

I fail to see this. You get one equation for each dynamical variable.

As for timing, it really depends on your application. Sometimes you have to take compilation times into account (if you really need to compile multiple times); sometimes you do not (because you only only compile once at the end of the day).

Do you have any questions issues demonstrated in your plots? If yes, I need more information, e.g., what is plotted, what do you expect to happen, etc.

I took a closer look at your code and found several issues. I strongly recommend that you try to fully understand everything I will be telling you (ask me if you have questions) and double-check all the mathematics of what you are doing (since I found several severe flukes and I only scratched the surface of your work).

Here are the main issues (in no particular order):

For any complex number z, you have:

Re(z) cos(arg(z)) + Im(z) sin(arg(z)) = |z| −Re(z) sin(arg(z)) + Im(z) cos(arg(z)) = 0

This simplifies and obviates several of your computations.
If you need to numerically or symbolically calculate the angle of a complex number, use atan2.
In symbolic computing, use actual fractions whenever possible, e.g., use Rational(2,3) instead of 2/3. The former will be handled exactly and possibly simplified; the latter will be only be approximated.
Most importantly: You can reduce your system to an ODE and you should do so. You do not get any additional dynamic variables or similar problems. Below I let SymPy solve the system for me, but you can also do this on paper (which may provide valuable insights). I described how to do this above.

Here is why using DDEs in your case is a computational nightmare: You initiate the DDE with a constant past and thus your initial estimate of all dynamical variables (y₀, y₁, …) via terms like (y(3)-y(3,t-delta_t))/delta_t is 0. Now, if you compute the derivative at this point using equat_dyna_model, this will usually give you a different value. Note that I used the constant past and the first time step only as an example: This is a problem that you have at every single time step.

If you are lucky, it may be that the two derivatives equalise over time, but I see no straightforward way to estimate how long this takes and whether it happens at all. Until the, the estimated derivative caught up with the real one, it will always lag behind. It will also always be a bit too low, being roughly delta_t/2 behind in time. This makes your results incorrect (and incomparable between different implementations – not that this really matters).

While this problem is similar in nature to that of initial discontinuities, it is not the same; it is worse. The problem of initial continuities is completely gone after a finite and usually short time. Your problem stays longer (and possibly forever) and has worse consequences. This is also the cause of severe numerical instabilities, e.g., some values going to −10⁸ during the first time step. I wouldn’t be surprised if this causes many other problems that I haven’t even thought of. It’s a miracle that you got this to produce some remotely reasonable solution in the first place (not that this helps you in any way).

Below is a solution of your problem that uses ODEs. Note that it is much faster (in fact, most of the time is spent on symbolic stuff before the actual integration right now) and there are no computational hiccups. Again, I strongly suggest that you do not continue working with this before understanding everything that happens here and all the changes I made:

from sympy import (
        Symbol, Rational,
        sin, cos, atan2, sqrt, pi,
        Eq, Matrix, Piecewise,
        solve
    )
import numpy as np
from jitcode import y,t,jitcode

# Voltage Source
def volt_source(f,v_rms,t_i):
    omega = 2*pi*f
    v_r = v_rms*sin(omega*t_i)
    v_y = v_rms*sin(omega*t_i-2*pi/3)
    v_b = v_rms*sin(omega*t_i+2*pi/3)
    A = Matrix([[1,-Rational(1,2),-Rational(1,2)],[0,sqrt(3)/2,-sqrt(3)/2]])
    B = Matrix([v_r,v_y,v_b])*sqrt(Rational(2,3))
    [v_a,v_b] = A*B
    theta = atan2(v_b,v_a)
    length = sqrt(v_a**2+v_b**2)
    return [length.simplify(),0,theta.simplify()]

# Parameters
f = 50
v_rms = 230
R_s = 0.435 # Stator resistance
R_r = 0.816 # Rotor resistance
L_ls = 0.002 # Stator inductance
L_lr = 0.002 # Rotor inductance
L_m = 0.0693
M = L_m # Mutual inductance
P = 2 # Number of Poles
J = 0.089 # Inertia
B = 0.01 # friction coefficient
T_l = 0.0 # Load torque
W_e = 2*pi*f
L_ss = L_ls+L_m
L_rr = L_lr+L_m
v_d = []
v_q = []

# Initial conditions

ydot = [ Symbol(f"ydot{i}") for i in range(5) ]
g = [ Symbol(f"g_{i}") for i in range(5) ]

vs = volt_source(f,v_rms,t)

helpers = [
    ( g[0], vs[1]-W_e*(y(1)*(L_ls+L_m)+y(2)*L_m)-R_s*y(0) ),
    ( g[1], vs[0]+W_e*(y(0)*(L_ls+L_m)+y(3)*L_m)-R_s*y(1) ),
    ( g[2], (W_e-y(4)*2/P)*(y(3)*(L_lr+L_m)+y(0)*L_m)-R_r*y(2) ),
    ( g[3],-(W_e-y(4)*2/P)*(y(2)*(L_lr+L_m)+y(1)*L_m)-R_r*y(3) ),
    ( g[4], ((((y(0)*(y(1)+y(2))*L_m+y(2)*L_lr-y(1)*(y(0)+y(3))*L_m+y(3)*L_lr)*(L_m/(L_lr+L_m))*(Rational(2,3)*P/2))-T_l)*(1/J)-y(4)*B/J) ),
]

# Solution
result = []
equations = [
        Eq( ydot[0], (-ydot[3]*L_m+g[0])/(L_ls+L_m) ),
        Eq( ydot[1], (-ydot[2]*L_m+g[1])/(L_ls+L_m) ),
        Eq( ydot[2], (-ydot[1]*L_m+g[2])/(L_lr+L_m) ),
        Eq( ydot[3], (-ydot[0]*L_m+g[3])/(L_lr+L_m) ),
        Eq( ydot[4], g[4] ),
    ]

eqs_solved = solve(equations,ydot)
rhs = [ eqs_solved[ydoti] for ydoti in ydot ]

delta_t = 1e-4
initial_state = np.zeros(5)
times = np.arange(0,1,delta_t)

solut = jitcode(rhs,helpers,verbose=False)
solut.set_integrator("dopri5")
solut.set_initial_value(initial_state)
for time in times:
    print( time, *solut.integrate(time) )

Thanks for the help and your clarification. But there are a number of points (for understanding): 1) "Now, if you compute the derivative at this point using equat_dyna_model, this will usually give you a different value. Note that I used the constant past and the first time step only as an example: This is a problem that you have at every single time step." - I agree with this, this is my mistake - the minimum solution is solut.constant_past (initial_state, time = 0.0) 2) But to be honest, I don’t understand what the problem is when writing not through a differential, but as a difference of a function with a delay, while the delay is equal to the step - solut.integrate_blindly (delta_t)

"Now, if you compute the derivative at this point using equat_dyna_model, this will usually give you a different value. Note that I used the constant past and the first time step only as an example: This is a problem that you have at every single time step." - I agree with this, this is my mistake - the minimum solution is solut.constant_past (initial_state, time = 0.0)

I fail to understand what your remark has to do with the quote or anything. What do you mean by minimum solution?

But to be honest, I don’t understand what the problem is when writing not through a differential, but as a difference of a function with a delay, while the delay is equal to the step - solut.integrate_blindly (delta_t)

First of all, the step you are talking about is the sampling step. The integration step is chosen adaptively and can be much lower. Even if the step were equal to the integration step, you would still suffer from the problem, unless you make your integration step very small.

The following analogy may help to understand your problem: Suppose you want to solve the equation x = 0.9·x + a. Instead of solving this analytically, you do this numerically and iteratively via the sequence x_n+1 = 0.9·x_n + a. Yes, eventually, x will approach the solution, but for most initial values x₀, it will take quite some time.

You are doing something like this: The x on the left-hand side is your ẏ; the x on the right-hand side is your approximation of ẏ via the finite difference. So, at best, it will take some time, till ẏ is right. To make things worse, a is changing all the time and it may change so quickly that your approach will never converge to the right solution, but only “chase” after it all the time.

Also please consider that if you apply some new numerical approach to find a solution to some problem, the default assumption is that it doesn’t work and it is up to you to convince your boss, peer-reviewers, etc. that your approach is actually working – it’s not up to them to elaborate a flaw to the last detail. Even then, the fact that your approach yields values of −10⁸ in the first integration step is a pretty strong argument against it.

I did not put it correctly - in the last iteration of the code I made a mistake (solut.constant_past (initial_state)) that caused initialization [.0, .0, .0, .0, .0] for each iteration, I decided to make a"minimal change" to solve this - solut.constant_past (initial_state, time = 0.0) (for time<=0.0)

I did not put it correctly - in the last iteration of the code I made a mistake (solut.constant_past (initial_state)) that caused initialization [.0, .0, .0, .0, .0] for each iteration, I decided to make a"minimal change" to solve this - solut.constant_past (initial_state, time = 0.0) (for time<=0.0)

There is no difference between the two. If the time argument for constant_past is not given, it defaults to 0.

default strange, in this case it is not equivalent for I.set_initial_value(initial,time=0.0) (jitcode)?

strange, in this case it is not equivalent for I.set_initial_value(initial,time=0.0) (jitcode)?

Except for the fact, that you have to provide not only an initial value, but an entire initial history when working with DDEs (which you choose to be constant with constant_past) the two are the same.

However, I have the strong feeling that I fail to understand your problem. Can you please elaborate where you see a contradiction?

neurophysik / jitcdde

Time dependent constant in function #16