crunching a ODE with Euler's Method doesn't work

Question

So I have to crunch a ODE with different methods and Δt, but suddenly my function just, stops working with Δt<1. The ODE is : T(t)'=R*T(t) - K*T(t-δ), with R and K constants, decimal and small constants at that.

The first part is crunching with δ=0, so it went smoothly and no issues, but now I'm crunching it with δ=5, and suddenly the ODE just doesn't work. With Δt=1 the program crunches the ODE just fine with Euler's method, but with Δt=.1, it suddenly stops working.

It doesn't crash or anything, but instead of calculating the numbers correctly, it doesn't seem to recognize the δ, or it assumes it as zero, and goes on like that, to the point that the graphic is equal to those generated with δ=0.

I'm working with Numpy and Matplotlib, and here is the code:

r = 0.1
R = 0.1
CeCw = 0.27
γ1 = 1
γ2 = 0.164
α = 0.612
δ = 5
Te0 = 1
tf = 240
Δt2 = np.arange(0, tf+.1, .1)

def d(Te, n):
    return R*Te[n] - α*γ1*CeCw*Te[n-δ]

Te2 = np.array([Te0])
for n in range(len(Δt2)-1):
    if n*.1<=δ:
        Te2 = np.append(Te2, 1)
    else:
        Te2 = np.append(Te2, Te2[n] + .1*d(Te2, n))

If I change that .1 to simply 1, it suddenly starts working again (albeit, it generates an incorrect graph since it isn't the proper method). I tried wrapping the numbers with Decimal() and using different functions, but somehow that .1 comes back to delete the δ.

EDIT: Moved forward and tried implementing Heun's Method with δ=5, and again, with Δt=.1 or lower, it behaves as if δ=0

Answer 1

This kind of equation is usually called Delay-differential equation and is not an ODE. Specialized solvers usually use an interpolation of sufficiently high error order to get the values at past times.

More specifically to your code, you are not implementing the delay correctly. The delay δ is a time difference, not an index difference. For the index difference m_δ you want for your kind of implementation that m_δ*h = δ. Then use m_δ to determine the length of the initial segment and the offset in the derivative function.

This could look like

δ = 5
mδ = 50
h = δ/mδ
Te0 = 1
tf = 240
Δt2 = np.arange(0, tf+h, h)

def d(Te, n):
    return R*Te[n] - α*γ1*CeCw*Te[n-mδ]

Te2 = np.array(len(Δt2)*[Te0], dtype=float)

for n in range(mδ,len(Δt2)-1):
    Te2[n+1] = Te2[n] + h*d(Te2, n)

This then indeed produces an oscillation

The same visual image is obtained using the Heun method with the time loop

for n in range(mδ,len(Δt2)-1):
    k1 = d(Te2, n)
    Te2[n+1] = Te2[n] + h*k1
    k2 = d(Te2, n+1)
    Te2[n+1] = Te2[n] + 0.5*h*(k1+k2)

The numerical values differ from the second decimal place on, as is to be expected with first and second order methods and h=0.1.