1
$\begingroup$

I'm trying to solve a coupled DDE. Here is my code

beta = 0.0005;
lambda = 1;
d = 0.1;
alpha = 0.000001;
a = 0.3;
p = 1 ;
k = 200;
u = 8;
c = 0.2;
b = 0.15;

Module[{sold, τ1 = 10, τ2 = 10}, 
  sold = 
    First[NDSolve[{
      Subscript[y, 1]'[t] == lambda - d Subscript[y, 1][t] - 
        beta Subscript[y, 1][t] Subscript[y, 3][t]/(1 + alpha Subscript[y, 3][t]), 
      Subscript[y, 1][t /; t <= 0] == 1,
      Subscript[y, 2]'[t] == 
        beta Subscript[y, 1][t - τ1] Subscript[y, 3][t - τ1]/
          (1 + alpha Subscript[y, 3][t - τ1]) - 
        a Subscript[y, 2][t] - p Subscript[y, 2][t] Subscript[y, 4][t],
      Subscript[y, 2][t /; t <= 0] == 1,
      Subscript[y, 3]'[t] == k Subscript[y, 2][t - τ2] - u Subscript[y, 3][t], 
      Subscript[y, 3][t /; t <= 0] == 5,
      Subscript[y, 4]'[t] == 
        c Subscript[y, 2][t] Subscript[y, 4][t] - b Subscript[y, 4][t], 
     Subscript[y, 4][t /; t <= 0] == 1
     }, 
     {Subscript[y, 1], Subscript[y, 2], Subscript[y, 3], Subscript[y, 4]}, 
     {t, 0, 200}]];

Plot[
  Evaluate[{
    Subscript[y, 1][t], Subscript[y, 2][t], 
    Subscript[y, 3][t], Subscript[y, 4][t]} /. sold], 
  {t, 0, 200}, 
  Filling -> Axis, AspectRatio -> 1]]

I want to get four separate, framed plots. Also, can anybody help me to rewrite this code in more easily understandable way?

Improve this question
$\endgroup$

1 Answer 1

Reset to default
3
$\begingroup$

I changed the subscripts to numerated variables - much easier to read IMO.

Module[{sold, τ1 = 10, τ2 = 10, plot}, 
 sold = First[
   NDSolve[{y1'[t] == 
      lambda - d y1[t] - beta y1[t] y3[t]/(1 + alpha y3[t]), 
     y1[t /; t <= 0] == 1, 
     y2'[t] == 
      beta y1[t - τ1] y3[
          t - τ1]/(1 + alpha y3[t - τ1]) - a y2[t] - 
       p y2[t] y4[t], y2[t /; t <= 0] == 1, 
     y3'[t] == k y2[t - τ2] - u y3[t], y3[t /; t <= 0] == 5, 
     y4'[t] == c y2[t] y4[t] - b y4[t], y4[t /; t <= 0] == 1}, {y1, 
     y2, y3, y4}, {t, 0, 200}]];

 plot = Plot[Evaluate[# /. sold], {t, 0, 200}, Filling -> Axis, 
     AspectRatio -> 1] & /@ {y1[t], y2[t], y3[t], y4[t]};

 GraphicsGrid[Partition[plot, 2], Frame -> All]]

enter image description here

Changes to get the boxed separate plots are to map the plot over the desired components, collecting them in the list plot, then using GraphicsGrid to display.

You can read up on the options for display using it, and check the "see also" there for other options to format sets of plots.

If you want to add labels, you can do something like this for the plot part (which can be extended to arbitrary items distinct to each result):

plot = MapThread[
   Plot[Evaluate[# /. sold], {t, 0, 200}, Filling -> Axis, AspectRatio -> 1, PlotLabel -> #2] &, 
        {{y1[t], y2[t], y3[t], y4[t]}, {"Y1", "Y2", "Y3", "Y4"}}];

Giving:

enter image description here

Improve this answer
$\endgroup$
4
  • $\begingroup$ Thanks dear for your time and useful comments. I'm still wondering about the separate plot. Looking for a more attractive graphs. Each one should be boxed separately, with line thickness higher than the default one, with axes clearly visible and can adjust the range more easily. $\endgroup$
    – zhk
    Feb 20, 2015 at 9:46
  • $\begingroup$ It hard to see the full picture because the plots are not showing the full output. $\endgroup$
    – zhk
    Feb 20, 2015 at 9:51
  • $\begingroup$ @mmm: Surely you're not suggesting readers iteratively fix the code at your direction? Consult the documentation re: options for Plot and GraphicsGrid (and the 'see also' within for other ways to format) and give it a whirl. If you get stuck, post back, we're here to help. $\endgroup$
    – ciao
    Feb 20, 2015 at 20:33
  • $\begingroup$ @rasher My issue has been resolved, thanks dear for your helpful responses. Now, to extract data from the above multiple line plot, I used 'Cases' but couldn't succeed. Any comment? $\endgroup$
    – zhk
    Feb 22, 2015 at 17:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.