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)=\mathrm{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.

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.*

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`

$\text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}}\to \text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathit{y}\left(\frac{\mathit{t}}{2}\right)$`ypdel`

$\text{\hspace{0.17em}}\to \text{\hspace{0.17em}}{\mathit{y}}^{\prime}\left(\mathit{t}-\pi \right)\text{\hspace{0.17em}}$

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

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

Finally, define the interval of integration $\left[{\mathit{t}}_{0\text{\hspace{0.17em}}}\text{\hspace{0.17em}\hspace{0.17em}}{\mathit{t}}_{\mathit{f}}\right]$ and solve the DDE using the `ddensd`

solver.

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

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])

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 %-------------------------------------------

[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.

`dde23`

| `ddensd`

| `ddesd`

| `deval`