DDE of Neutral Type

This example shows how to use ddensd to solve a neutral DDE (delay differential equation), where delays appear in derivative terms. The problem was originally presented by Paul [1].

The equation is

${\mathit{y}}^{\prime }\left(\mathit{t}\right)=1+\mathit{y}\left(\mathit{t}\right)-2{\mathit{y}\left(\frac{\mathit{t}}{2}\right)}^2-{\mathit{y}}^{\prime }\left(\mathit{t}-\pi \right)$.

The history function is$\mathit{y}\left(\mathit{t}\right)=\cos \left(\mathit{t}\right)$ for $\mathit{t}\le 0$.

Since the equation has time delays in a ${\mathit{y}}^{\prime }$ term, the equation is called a neutral DDE. If the time delays are only present in $\mathit{y}$ terms, then the equation would be a constant or state-dependent DDE, depending on what form the time delays have.

To solve this equation in MATLAB®, you need to code the equation, delays, and history before calling the delay differential equation solver ddensd. You either can include these as local functions at the end of a file (as done here), or save them as separate files in a directory on the MATLAB path.

Code Delays

First, write functions to define the delays in the equation. The first term in the equation with a delay is $\mathit{y}\left(\frac{\mathit{t}}{2}\right)$.

function dy = dely(t,y) 
    dy = t/2;
end

The other term in the equation with a delay is ${\mathit{y}}^{\prime }\left(\mathit{t}-\pi \right)$.

function dyp = delyp(t,y) 
    dyp = t-pi;
end

In this example, only one delay for $\mathit{y}$ and one delay for ${\mathit{y}}^{\prime }$ are present. If there were more delays, then you can add them in these same function files, so that the functions return vectors instead of scalars.

Note: All functions are included as local functions at the end of the example.

Code Equation

Now, create a function to code the equation. This function should have the signature yp = ddefun(t,y,ydel,ypdel), where:

• t is time (independent variable).

• y is the solution (dependent variable).

• ydel contains the delays for $\mathit{y}$.

• ypdel contains the delays for ${\mathit{y}}^{\prime }=\frac{\mathrm{dy}}{\mathrm{dt}}$.

These inputs are automatically passed to the function by the solver, but the variable names determine how you code the equation. In this case:

• ydel$\to \mathit{y}\left(\frac{\mathit{t}}{2}\right)$

• ypdel$\to {\mathit{y}}^{\prime }\left(\mathit{t}-\pi \right)$

function yp = ddefun(t,y,ydel,ypdel) 
    yp = 1 + y - 2*ydel^2 - ypdel;
end

Code Solution History

Next, create a function to define the solution history. The solution history is the solution for times $\mathit{t}\le {\mathit{t}}_0$.

function y = history(t)
    y = cos(t);
end

Solve Equation

Finally, define the interval of integration $\left[{\mathit{t}}_0{\mathit{t}}_{\mathit{f}}\right]$ and solve the DDE using the ddensd solver.

tspan = [0 pi];
sol = ddensd(@ddefun, @dely, @delyp, @history, [0,pi]);

Plot Solution

The solution structure sol has the fields sol.x and sol.y that contain the internal time steps taken by the solver and corresponding solutions at those times. However, you can use deval to evaluate the solution at the specific points.

Evaluate the solution at 20 equally spaced points between 0 and pi.

tn = linspace(0,pi,20);
yn = deval(sol,tn);

Plot the calculated solution and history against the analytical solution.

th = linspace(-pi,0);
yh = history(th);
ta = linspace(0,pi);
ya = cos(ta);

plot(th,yh,tn,yn,'o',ta,ya)
legend('History','Numerical','Analytical','Location','NorthWest')
xlabel('Time t')
ylabel('Solution y')
title('Example of Paul with 1 Equation and 2 Delay Functions')
axis([-3.5 3.5 -1.5 1.5])

Local Functions

Listed here are the local helper functions that the DDE solver ddensd calls to calculate the solution. Alternatively, you can save these functions as their own files in a directory on the MATLAB path.

function yp = ddefun(t,y,ydel,ypdel) % equation being solved
    yp = 1 + y - 2*ydel^2 - ypdel;
end
%-------------------------------------------
function dy = dely(t,y) % delay for y
    dy = t/2;
end
%-------------------------------------------
function dyp = delyp(t,y) % delay for y'
    dyp = t-pi;
end
%-------------------------------------------
function y = history(t) % history function for t < 0
    y = cos(t);
end
%-------------------------------------------

References

[1] Paul, C.A.H. “A Test Set of Functional Differential Equations.” Numerical Analysis Reports. No. 243. Manchester, UK: Math Department, University of Manchester, 1994.