How to perturb a Dynamic System?

Question

I'm trying to model a basic feedback system with delayed feedback. I've done the initial setup and now want to add a few more advanced features to my system.

Currently, it's just a simple delayed-differential equation with manipulatable delay to show how increasing the delay can send the system out of equilibrium.

What I'd like to add, however, is the ability for the system to be perturbed at a particular point in time. For example, in the Predator-Prey model this might mean that the system acts exactly as normal at all times against the given initial condition, but suddenly at t = T, where T is the perturbation point, either the number of predators or prey is shocked.

I want to model how this affects equilibrium for various parameters. Any good places to start?

Thanks again!

Edit:

What I have initially is -

Manipulate[
 HumanFaucetSystem = 
  NDSolve[{y'[time] + a*y[time - delay] == 0, 
    y[time /; time <= 0] == deviation}, y, {time, 0, 200*delay + 1}];
 Plot[Evaluate[y[x] /. HumanFaucetSystem], {x, 0, 100*delay + 1}, 
  PlotRange -> All], {{delay, 0, "Delay"}, 0, 
  15}, {{a, 1, "Multiplier"}, -5, 
  5}, {{deviation, 25, "Initial Deviation"}, 0, 100}]

I want to add another interactive parameter to the system that allows the user to schedule a specific time point at which a certain perturbation will occur. This can be positive or negative in the value of y.

On this logic, I want to expand it to a two-agent delayed system that I've modeled already but add this feature.

Answer 1

I haven't thought about this for delay differential equations, but for initial value problems, you can just think of the perturbation as a new initial value problem, then the only issue is stitching together the interpolating functions with Witch. Since you mention predator-prey systems lets use logistic growth as the example:

sol1 = First@With[{r = 0.5, k = 10, x0 = 5},
   NDSolve[{x'[t] == r x[t] (1.0 - x[t]/k), x[0] == x0}, x, {t, 0, 200}]
]

and the solution after the perturbation:

sol2 = First@Module[{r = 0.5, k = 10, pert = 1.1, x0},
   x0 = pert ( x[200] /. sol1);
   NDSolve[{x'[t] == r x[t] (1.0 - x[t]/k), x[200] == x0}, x, {t, 200, 400}]
]

Now we use witch to put the solutions together

sol[t_] := Which[0 <= t <= 200, x[t] /. sol1, 200 < t <= 400, x[t] /. sol2]

which we can visualize with Plot

Plot[sol[t], {t, 0, 400}, PlotRange -> All]

Hope that helps. Again sorry for not dealing with the delay differential equations ... would just need to think about what to do with the more complex initial conditions.

Extension

What follows is an example that does more than one perturbation, I move to a multispecies model, but the logic will work for single dimension models as well

So we have a basic Rosenzweig-MaCarthur Predator Prey Model

Clear[fx, fy]
fx[x_, y_] := r x (1 - x/k) - (a x y)/(1 + a h x)
fy[x_, y_] := (e a x y)/(1 + a h x) - m y

Also as the information for NDSolve can be a bit awkward to specify each time lets make a function that sets the initial conditions and returns a model list in the form that the solver expects

Clear[fweb]
fweb[{x0_, y0_}] := {
  x'[t] == fx[x[t], y[t]],
  y'[t] == fy[x[t], y[t]],
  x[0] == x0,
  y[0] == y0
}

we put the parameters in a rule list so that we can change the parameters easily

param = Dispatch[{
  r -> 0.5,
  k -> 1.0,
  a -> 1.0,
  h -> 1.0,
  e -> 1.0,
  m -> 0.2
}];

Perturbation function that we call with FoldList

Clear[perturbLastValue]
perturbLastValue[sol_, pert_, deltat_] := Module[{r, c},
  {r, c} = {x[deltat], y[deltat]} /. sol;
  {
    If[r > 0, Max[r + If[Chop[pert] > 0, RandomVariate[NormalDistribution[0, pert]], 0], 0], 0],
    If[c > 0, Max[c + If[Chop[pert] > 0, RandomVariate[NormalDistribution[0, pert]], 0], 0], 0]
  }
]

Now we make the solutions for these time intervals

Clear[perturbSol]
perturbSol[x0_, pert_, tInts_] :=
  FoldList[
    First@NDSolve[fweb[perturbLastValue[#1, pert, First@#2]] /. param, {x, y}, {t, First@#2, Last@#2}, AccuracyGoal -> 20, MaxSteps -> 100000] &,
    First@NDSolve[fweb[x0] /. param, {x, y}, {t, First@tInts[[1]], Last@tInts[[1]]}, AccuracyGoal -> 20, MaxSteps -> 100000],
    Rest[tInts]
];

Now we make the time intervals between perturbation events and stick up the solution into a function we can use to get the values from

Clear[x, y, sol, cond, tInts, sols]
With[{deltat = 1},
  tInts = Table[{tstart, tstart + deltat}, {tstart, 0, 80*deltat, deltat}];
  sols = perturbSol[{0.25, 0.25}, 0.05, tInts];
  cond = Which @@ Flatten[Table[{First@tInts[[i]] <= t < Last@tInts[[i]], {x[t], y[t]} /. sols[[i]]}, {i, Length[tInts]}], 1];
  sol[tin_] := cond /. t -> tin
];

we can use the solution as follows

Plot[sol[t], {t, 0, 800}, PlotRange -> All]

Answer 2

Now that Mathematica has added WhenEvent we have the super sweet solution that requires non of this ugly boiler plate.

For the single perturbation case we have the following:

Module[{r = 0.5, k = 10, x0 = 5, perturb, sol},
  perturb = WhenEvent[Mod[t, 200], x[t] -> 1.1 x[t]];
  sol = NDSolveValue[{{x'[t] == r x[t] (1.0 - x[t]/k), x[0] == x0}, perturb}, x, {t, 0, 400}];
  Plot[sol[t], {t, 0, 400}, PlotRange -> All]
]

And if we want to make frequent perturbations we can do:

Module[{r = 0.5, k = 10, x0 = 5, perturb, sol},
  perturb = WhenEvent[Mod[t, 1], x[t] -> x[t] + RandomReal[]];
  sol = NDSolveValue[{{x'[t] == r x[t] (1.0 - x[t]/k), x[0] == x0}, perturb}, x, {t, 0, 400}];
  Plot[sol[t], {t, 0, 400}, PlotRange -> All]
]

Such an amazing feature as this tremendous simplification shows!

Answer 3

Gabriel, your answer really helped. I think I've implemented the functionality I was going for:

 Manipulate[
 System = First@
   With[{d = delay, a = multiplier, y0 = deviation}, 
    NDSolve[{y'[time] + a*y[time - d] == 0, 
      y[time /; time <= 0] == y0}, y, {time, 0, timePerturbed}]];
 PerturbedSystem = 
  First@Module[{d = delay, a = multiplier, pert = perturbation, y0}, 
    y0 = pert (y[timePerturbed] /. System);
    NDSolve[{y'[time] + a*y[time - d] == 0, 
      y[time /; time <= timePerturbed] == y0}, 
     y, {time, timePerturbed, 200*delay + 1}]];
 sol[time_] := 
  Which [ 0 <= time <= timePerturbed, y[time] /. System, 
   timePerturbed < time <= 200*delay + 1, 
   y[time] /. PerturbedSystem];
 Plot[sol[t], {t, 0, 200*delay + 1}, PlotRange -> All],
 {{delay, 0, "Delay"}, 0, 15}, {{multiplier, 1, "Multiplier"}, -5, 
  5}, {{deviation, 25, "Initial Deviation"}, 0, 100},
 {{perturbation, 0, "Amount Perturbed"}, 0, 100},
 {{timePerturbed, 0, "Time at Perturbation"}, 0, 100*delay + 1}]

The output isn't really what I was looking for, but I think that has to do with the plotrange.

Thanks!

Oh and for anyone else looking to apply this to a DDE, the logic is similar to what Gabriel posted, but you just have to give the "initial history" rather than the initial point.