2

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.

zhk
  • 11,939
  • 1
  • 22
  • 38
  • "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. – rhermans Oct 07 '15 at 14:13
  • 3
    NIntegrate can't handle DiracDelta (at least currently), this is mentioned in Possible Issues of the document of DiracDelta. – xzczd Oct 08 '15 at 03:22
  • @xzczd If we replace 'DiracDelta' by a 'Sin' or 'Cos', I am getting the same error message. – zhk Oct 08 '15 at 05:46
  • 1
    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. – xzczd Oct 08 '15 at 09:11

1 Answers1

4

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]"}]

enter image description here

  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]"}]

enter image description here

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

enter image description here

Mariusz Iwaniuk
  • 13,841
  • 1
  • 25
  • 41
  • 1
    Oh, I didn't expect the integration can be calculated symbolicly! … But why you assume v[0]=1; x[0]=2? – xzczd Oct 11 '15 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. – Michael E2 Oct 11 '15 at 13:03
  • @MichaelE2 ,Thanks .Improving the post. – Mariusz Iwaniuk Oct 11 '15 at 13:06
  • @xzczd.Without this assumption NDSolve can't solve.Other numbers can be set up v[0]=?,x[0]=?. ? = Any number. – Mariusz Iwaniuk Oct 11 '15 at 13:24
  • 1
    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 :) – xzczd Oct 11 '15 at 13:33
  • @xzczd.Yes You are right,my mistaken. – Mariusz Iwaniuk Oct 11 '15 at 14:59
  • @I_Mariusz Thanks for your efforts. How can I plot space curves x[t] vs y[t] vs z[t]? – zhk Oct 14 '15 at 05:20
  • @MMM. Use :ParametricPlot3D[{x[t], y[t], z[t]} /. sol, {t, 0, 20}, PlotStyle -> {Orange, Thickness[0.015]}, BoxRatios -> {1, 1, 1}] – Mariusz Iwaniuk Oct 14 '15 at 10:06