CODEWITHSUNDEEP
pythonmathdifferential-equationsnumber-theory
Simd
Simd

Reputation: 21233

How to solve a delay differential equation numerically

I would like to compute the Buchstab function numerically. It is defined by the delay differential equation:

enter image description here

How can I compute this numerically efficiently?

Upvotes: 1

Views: 770

Answers (2)

AmirHosein Sadeghimanesh
AmirHosein Sadeghimanesh

Reputation: 478

Well, first you can simplify your equation. Note that

enter image description here

Now if I denote "dw(u)/du" with a w that has a dot above it, your equation simplifies to the following one.

enter image description here

Now it should be fairly easy to use a simple numerical method for ODEs such as Euler method (https://en.wikipedia.org/wiki/Euler_method) to write a simple for-loop in Python to solve your DDE numerically. But here I provide a Matlab code to solve this and generate the same plot that Lutz Lehmann and Wikipedia page of Buchstab function have shown. There are several commands for DDEs in Matlab, but I like to use ddesd even though your delay terms are not state dependent. For a short video about how to use this command in Matlab see this YouTube video https://youtu.be/s3G1aCCE-4Q.

sol = ddesd( @ddefun, @delay, @history, [ 2, 4 ] );
u = sol.x;
w = sol.y;
uhistory = linspace( 1, 2 );
whistory = 1 ./ uhistory;
plot( uhistory, whistory, u, w, 'b' );
xlabel( 'u' );
ylabel( 'w(u)' );

% local functions
function dxdt = ddefun( u, w, Z )
    dxdt = ( Z( 1, 1 ) - w ) / u;
end
function d = delay( u, w )
    d = u - 1;
end
function v = history( u )
    v = 1 / u;
end

And here is the output figure.

enter image description here

Upvotes: 3

Lutz Lehmann
Lutz Lehmann

Reputation: 25992

To get a general feeling of how DDE integration works, I'll give some code, based on the low-order Heun method (to avoid uninteresting details while still being marginally useful).

In the numerical integration the previous values are treated as a function of time like any other time-depending term. As there is not really a functional expression for it, the solution so-far will be used as a function table for interpolation. The interpolation error order should be as high as the error order of the ODE integrator, which is easy to arrange for low-order methods, but will require extra effort for higher order methods. The solve_ivp stepper classes provide such a "dense output" interpolation per step that can be assembled into a function for the currently existing integration interval.

So after the theory the praxis. Select step size h=0.05, convert the given history function into the start of the solution function table

u=1
u_arr = []
w_arr = []
while u<2+0.5*h:
    u_arr.append(u)
    w_arr.append(1/u)
    u += h

Then solve the equation, for the delayed value use interpolation in the function table, here using numpy.interp. There are other functions with more options in `scipy.interpolate.

Note that h needs to be smaller than the smallest delay, so that the delayed values are from a previous step. Which is the case here.

u = u_arr[-1]
w = w_arr[-1]
while u < 4:
    k1 = (-w + np.interp(u-1,u_arr,w_arr))/u
    us, ws = u+h, w+h*k1
    k2 = (-ws + np.interp(us-1,u_arr,w_arr))/us
    u,w = us, w+0.5*h*(k1+k2)
    u_arr.append(us)
    w_arr.append(ws)

Now the numerical approximation can be further processed, for instance plotted.

plt.plot(u_arr,w_arr); plt.grid(); plt.show()

enter image description here

Upvotes: 3

Related Questions