1
$\begingroup$

I am trying to solve a system of ODE. I am new to Mathematica enter image description here enter image description here enter image description here enter image description here 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] == 0, y2[0] == 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] == 0, y2[0] == 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: NDSolve::nrnum1

enter image description here

Does anyone have any idea how to resolve these issues? Any help would be appreciated. :)

$\endgroup$
8
  • $\begingroup$ I guess "/.s" should read "/.sol" And what plot are you expecting? $\endgroup$ Feb 22 at 11:04
  • $\begingroup$ @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 $\endgroup$
    – Chris
    Feb 22 at 11:50
  • $\begingroup$ 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. $\endgroup$
    – bbgodfrey
    Feb 22 at 13:11
  • $\begingroup$ @bbgodfrey Thank you for the tip :) $\endgroup$
    – Chris
    Feb 23 at 2:18
  • $\begingroup$ @Chris Your model is not same as in a paper linked. Why do you expect to get same picture as in Figure 1? $\endgroup$ Feb 23 at 3:47

1 Answer 1

0
$\begingroup$

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[0][0] == 10^-10}, Table[y[i][0] == 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[[1]]], {t, 0, 10}, 
 PlotLabel -> "Layer Coverage Dynamics", AxesLabel -> {t, Coverage}, 
 PlotRange -> All, 
 PlotLegends -> Table[Subscript[\[Theta], i], {i, 0, 10}]]

Figure 1

$\endgroup$
2
  • $\begingroup$ 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 $\endgroup$
    – Chris
    Feb 23 at 15:52
  • $\begingroup$ @Chris Sorry, code has been updated. Try again. $\endgroup$ Feb 23 at 16:35

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.