2
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ "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. $\endgroup$ – rhermans Oct 7 '15 at 14:13
  • 3
    $\begingroup$ NIntegrate can't handle DiracDelta (at least currently), this is mentioned in Possible Issues of the document of DiracDelta. $\endgroup$ – xzczd Oct 8 '15 at 3:22
  • $\begingroup$ @xzczd If we replace 'DiracDelta' by a 'Sin' or 'Cos', I am getting the same error message. $\endgroup$ – zhk Oct 8 '15 at 5:46
  • 1
    $\begingroup$ 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. $\endgroup$ – xzczd Oct 8 '15 at 9:11
4
+50
$\begingroup$

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

$\endgroup$
  • 1
    $\begingroup$ Oh, I didn't expect the integration can be calculated symbolicly! … But why you assume v[0]=1; x[0]=2? $\endgroup$ – xzczd Oct 11 '15 at 12:53
  • $\begingroup$ 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. $\endgroup$ – Michael E2 Oct 11 '15 at 13:03
  • $\begingroup$ @MichaelE2 ,Thanks .Improving the post. $\endgroup$ – Mariusz Iwaniuk Oct 11 '15 at 13:06
  • $\begingroup$ @xzczd.Without this assumption NDSolve can't solve.Other numbers can be set up v[0]=?,x[0]=?. ? = Any number. $\endgroup$ – Mariusz Iwaniuk Oct 11 '15 at 13:24
  • 1
    $\begingroup$ 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 :) $\endgroup$ – xzczd Oct 11 '15 at 13:33

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.