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]
Input
? Is that what you were looking for to make it interactive and require user input?? $\endgroup$