# How to perturb a Dynamic System?

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.

-
Could you post some code, please? –  belisarius Oct 15 '12 at 18:06
Edited, sorry about that. –  Ashok M Rao Oct 15 '12 at 18:11
Oh, no need to apologize! Posting code serves three purposes: 1) Focus on your problem 2) relieve answerers from typing it and 3) understand your Mma abilities so the answer could be "downearthed" in case you were a beginner –  belisarius Oct 15 '12 at 18:14
@AshokMRao Have you tried Input? Is that what you were looking for to make it interactive and require user input?? –  drN Oct 15 '12 at 20:31
That's one way of doing it, but what I really want to know is how to perturb the system from a particular value in the middle of a function. Like a break. After I know how to do this, I can program it in with either Manipulate or Input, I think? –  Ashok M Rao Oct 15 '12 at 21:24

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]

-
Hey, thanks for the help! I'll try this out with a normal differential equation and try to work it with a delayed one. I don't think it should be all too different. –  Ashok M Rao Oct 16 '12 at 14:10
Tell me how it goes! –  Gabriel Oct 16 '12 at 15:45
Though I'm wondering, is there a way to create multiple perturbations like this, for example some perturbation every N seconds? Maybe this could be kind of like a piecewise function, f or f + perturbation if a given time? –  Ashok M Rao Oct 16 '12 at 18:21
Yeah I do the continual perturbations most of the time, you just make a list of the conditions and functions and use Apply with Which I could clean up a notebook and make a link if you want. –  Gabriel Oct 16 '12 at 20:07
I have added a more involved example that does multiple perturbations –  Gabriel Oct 18 '12 at 21:14

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!

-

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.

-