# System of ODE with piecewise functionals

I am trying to solve a system of ODE. I am new to Mathematica    What I am trying to get is the plot labelled 1 & 2.

τ = 1;
A = 0.98;
equa = {y1'[t] == (
y0 - y1[t])/τ + α1[t]*(y1[t] - y2[t])/τ,
y2'[t] == (
y1[t] - y2[t])/τ + α2[t]*
y2[t]/τ - α1[t]*(y1[t] - y2[t])/τ,
y0 == 1, y1 == 0, y2 == 0};


The piecewise functions are alpha1 and alpha2. There are also piecewise functions within alpha1 and alpha2. I only tried simulating up to n=2.

d1[t_] :=
Piecewise[{{(y1[t])^0.5, y1[t] < 0.5}, {(1 - y1[t])^0.5,
y1[t] > 0.5}}];
d2[t_] :=
Piecewise[{{(y2[t])^0.5, y2[t] < 0.5}, {(1 - y2[t])^0.5,
y2[t] > 0.5}}];
α1[t_] :=
Piecewise[{{0, d1[t] = 0}, {(A* d1[t])/(d1[t] + d2[t]),
d1[t] != 0}}];
α2[t_] := Piecewise[{{0, d2[t] = 0}, {A, d2[t] != 0}}];


I use NDSolve to solve the system of ODE

sol = NDSolve[equa,{y1,y2},{t,0,10}];
Plot[Evaluate[{y1[t], y2[t]} /. s], {t, 0, 10},
PlotLabel -> "Layer Coverage Dynamics", AxesLabel -> {t, Coverage},
PlotRange -> All]


This returns a plot, however, it is not the plot that I am expecting. It seems to me that the Piecewise portion in the system of ODE is not taken into account during NDSolve. I am new to this so please bear with me.

Update: So I edited my code according to the comments

τ = 1;
A = 0.98;
d1[t_] :=
Piecewise[{{(y1[t])^0.5, y1[t] < 0.5}, {(1 - y1[t])^0.5,
y1[t] > 0.5}}];
d2[t_] :=
Piecewise[{{(y2[t])^0.5, y2[t] < 0.5}, {(1 - y2[t])^0.5,
y2[t] > 0.5}}];
α1[t_] :=
Piecewise[{{0, d1[t] == 0}, {(A* d1[t])/(d1[t] + d2[t]),
d1[t] != 0}}];
α2[t_] := Piecewise[{{0, d2[t] == 0}, {A, d2[t] != 0}}];
equa = {y1'[t] == (
y0 - y1[t])/τ + α1[t]*(y1[t] - y2[t])/τ,
y2'[t] == (
y1[t] - y2[t])/τ + α2[t]*
y2[t]/τ - α1[t]*(y1[t] - y2[t])/τ,
y0 == 1, y1 == 0, y2 == 0};
s = NDSolve[equa, {y1, y2}, {t, 0, 10}];
Plot[Evaluate[{y1[t], y2[t]} /. s], {t, 0, 10},
PlotLabel -> "Layer Coverage Dynamics", AxesLabel -> {t, Coverage},
PlotRange -> All]


Now, it returns 2 NDSolve errors:  Does anyone have any idea how to resolve these issues? Any help would be appreciated. :)

• I guess "/.s" should read "/.sol" And what plot are you expecting? Feb 22 at 11:04
• @Daniel Huber Thank you for your reply. what I am expecting can be found here: link However, the y2[t] graph in my code shows a very steep initial transient and what I am expecting to see would be a gentler slope Feb 22 at 11:50
• Piecewise[{{0, d2[t] = 0}, {A, d2[t] != 0}}] should be Piecewise[{{0, d2[t] == 0}, {A, d2[t] != 0}}] and similarly elsewhere. Do not confuse Set with Equal. Feb 22 at 13:11
• @bbgodfrey Thank you for the tip :) Feb 23 at 2:18
• @Chris Your model is not same as in a paper linked. Why do you expect to get same picture as in Figure 1? Feb 23 at 3:47

We use name for function $$\theta_i=y[i]$$,then code for $$i=0,1,...,15$$ can be written as follows

tau = 1;
A = .98; nn = 15; Y = Table[y[i][t], {i, -1, nn + 1}]; Y1 =
Table[y[i]'[t], {i, -1, nn + 1}];
y[-1][t_] := 1;
y[nn + 1][t_] := 0; ic =
Join[{y == 10^-10}, Table[y[i] == 10^-10, {i, nn}]]; var =
Table[y[i], {i, 0, nn}];
d[x_] := Piecewise[{{0, x <= 0}, {Sqrt[Abs[x]],
x <= 1/2}, {Sqrt[Abs[1 - x]], x > 1/2}, {0, x >= 1}}];
alpha[x_, z_] :=
If[d[x] == 0, 0, If[d[z] == 0, 0, (A*d[x])/(d[x] + d[z])]];
eq = Table[
Y1[[n]] == (Y[[n - 1]] - Y[[n]])/tau +
alpha[Y[[n]], Y[[n + 1]]]/tau (Y[[n]] - Y[[n + 1]]) -
alpha[Y[[n - 1]], Y[[n]]]/tau (Y[[n - 1]] - Y[[n]]), {n, 2,
nn + 2}];

soln = NDSolve[Join[eq, ic], var, {t, 0, 10}];


Visualization

Plot[Evaluate[Drop[Drop[Y, 1], -1] /. soln[]], {t, 0, 10},
PlotLabel -> "Layer Coverage Dynamics", AxesLabel -> {t, Coverage},
PlotRange -> All,
PlotLegends -> Table[Subscript[\[Theta], i], {i, 0, 10}]] • Thank you for the answer! I tried your code on my notebook but I was not able to get the same plot. I got an error NDSolve::ndsz instead. Do you have any packages loaded initially? I really appreciate your help Feb 23 at 15:52
• @Chris Sorry, code has been updated. Try again. Feb 23 at 16:35