NIntegrate::inumr: "The integrand E^(-0.9τ1) DiracDelta[t-τ1] v[t-τ1] x[t-τ1] has evaluated to non-numerical values "

I am facing an interesting situation over here. The aim is to solve a system of IVP having an integral over a delay range. Here is my try

beta = 0.0012;
lambda = 2;
d = .10;
alpha = 0.002;
a = 0.5;
p = 5.6 ;
k = 70;
b = 2;
c = 40;
m1 = 0.9;
m2 = 0.9;
q = 5.6;
sigma = .0005;

First[NDSolve[{
x'[t] == lambda - d x[t] - beta x[t] v[t], x[t /; t <= 0] == 3,
y'[t] ==
beta *NIntegrate[
Exp[-m1*τ1]*DiracDelta[t - τ1]*x[t - τ1]*
v[t - τ1], {τ1, 0, Infinity}] - a y[t] -
alpha w[t] y[t], y[t /; t <= 0] == 6,
z'[t] == a w[t] y[t] - b z[t], z[t /; t <= 0] == 3,
v'[t] ==
k NIntegrate[
Exp[-m2*τ1]*DiracDelta[τ1]*y[t - τ1], {τ1,
0, Infinity}] - p v[t], v[t /; t <= 0] == 149,
w'[t] == c  z[t] - q w[t], y5[t /; t <= 0] == 1
},
{x, y, z, v, w},
{t, 0, 20}]];


which gives this error,

NIntegrate::inumr: The integrand 0.9 E^(-0.9 [Tau]1) v[t-[Tau]1] x[t-[Tau]1] has evaluated to non-numerical values for all sampling points in the region

I have no clue what is going on. Please assist.

• "non-numerical values" comes from NIntegrate , change for Integrate if symbolic and not numerical. Somebody else should comment about the validity of z[t /; t <= 0] inside DSolve. Commented Oct 7, 2015 at 14:13
• NIntegrate can't handle DiracDelta (at least currently), this is mentioned in Possible Issues of the document of DiracDelta. Commented Oct 8, 2015 at 3:22
• @xzczd If we replace 'DiracDelta' by a 'Sin' or 'Cos', I am getting the same error message.
– zhk
Commented Oct 8, 2015 at 5:46
• Because NDSolve also has trouble in handling nonlinear delay integro-differential equations 囧. If the equation is linear, one can make use of LaplaceTransform(this is a example). But as far as I can tell, no one has ever managed to resolve a nonlinear case in this site. This answer is slightly related. Commented Oct 8, 2015 at 9:11

The problem my solution is that I do not know whether it is correct.Maybe someone will correct my solution.

Clear["Global*"]
beta = 3/2500;
lambda = 2;
d = 1/10;
alpha = 1/500;
a = 1/2;
p = 28/5;
k = 70;
b = 2;
c = 40;
m1 = 9/10;
m2 = 9/10;
q = 28/5;
sigma = 1/2000;


I calculate first integral separate:

Integrate[Exp[-m1*τ1]*DiracDelta[t - τ1]*x[t - τ1]*
v[t - τ1], {τ1, 0, Infinity},GenerateConditions -> False]


E^(-9 t/10) HeavisideTheta[t] v[0] x[0]

The [t /; t <= 0] isn't necessary because it's no longer a delay DE.

Initial conditions are: v[0]=149, x[0]=3 and assume that HeavisideTheta[t] is Piecewise[{{0, t < 0}, {1, t > 0},{Infinity, t == 0}}](Yes I now is "undefined" at point t=0' but it does not change anything) then:

 V[0] = 149;
X[0] = 3;
intY = E^(-9 t/10)*Piecewise[{{0, t < 0}, {1, t > 0},{Infinity, t == 0}}]*V[0]*X[0]


I calculate second integral:

  intV = Integrate[Exp[-m2*τ1]*DiracDelta[τ1]*y[t - τ1], {τ1, 0, Infinity}, GenerateConditions -> False]


We have:

-(-1 + HeavisideTheta[0]) y[t]

Assume that Limit[HeavisideTheta[x], x -> 0]=1 then intV=0.

Or HeavisideTheta[0] is often taken to be 1/2 then intV=y[t]/2

  intV = 0;
sol =
With[{ϵ = 10^-10},
First[NDSolve[{x'[t] == lambda - d x[t] - beta x[t] v[t],
x[ϵ] == 3,
y'[t] == beta*intY - a y[t] - alpha w[t] y[t],
y[ϵ] == 6, z'[t] == a w[t] y[t] - b z[t],
z[ϵ] == 3, v'[t] == k*intV - p v[t],
v[ϵ] == 149, w'[t] == c z[t] - q w[t],
w[ϵ] == 1}, {x, y, z, v, w}, {t, ϵ,
20}]]];


Plots with assuming that: intV=0

  Plot[Evaluate[{x[t], y[t], z[t], v[t], w[t]} /. sol], {t, 0, 20},
PlotRange -> All, PlotLegends -> {"x[t]", "y[t]", "z[t]", "w[t]"}]


  Plot[Evaluate[{x[t], y[t], z[t], v[t], w[t]} /. sol], {t, 0,
20}, PlotRange -> {{0, 7}, {0, 500}},
PlotLegends -> {"x[t]", "y[t]", "z[t]", "w[t]"}]


Space curve for: {x[t],y[t],z[t]}

  ParametricPlot3D[{x[t], y[t], z[t]} /. sol, {t, 0, 20},
PlotStyle -> {Orange, Thickness[0.015]}, BoxRatios -> {1, 1, 1},
PlotLegends -> {"space curve"}, AxesLabel -> {x, y, z}]


• Oh, I didn't expect the integration can be calculated symbolicly! … But why you assume v[0]=1; x[0]=2? Commented Oct 11, 2015 at 12:53
• HeavisideTheta[0] is often taken to be 1/2, although I believe there is not universal agreement about that (hence it being undefined in M). That would make intV to be y[t]/2. Commented Oct 11, 2015 at 13:03
• @MichaelE2 ,Thanks .Improving the post. Commented Oct 11, 2015 at 13:06
• @xzczd.Without this assumption NDSolve can't solve.Other numbers can be set up v[0]=?,x[0]=?. ? = Any number. Commented Oct 11, 2015 at 13:24
• So there's no special meaning to assume the values? Then you don't need to assume, they're already given as the initial conditions: V[0] = 149; X[0] = 3;. Also, the WorkingPrecision option can be taken away, and the [t /; t <= ϵ] isn't necessary because it's no longer a delay DE :) Commented Oct 11, 2015 at 13:33