0
$\begingroup$

I am trying to plot the solution given in the code with respect to "delc", not w.r.t time "t". I don't know what to do for getting the plot between solution and value "delc" which vary 0 to 2. The initial condition are also included in the text. The only thing that i want is, plot between solution and "delc" for any arbitrary value of "t". If anyone can solve this is most welcome.

w1 = 1;
gma1 = 0.005;
n1 = 1;
  gma2 = 0.005;
  G1 = 0.005;
  k1 = .1;
k2 = 0.1;
a1 = 0.07;
a2 = 0.58;
k0 = 0.1;
Q1 = 1.268;
del0 = 1;
 N1 = 1;
    ome = 1;
     M1 = del0*(1 - Cos[ome*t]);
s = ParametricNDSolveValue[{V11'[t] - V21[t]*w1 - V12[t]*w1 == 0, 
V12'[t] - V22[t]*w1 + w1*V11[t] + gma1*V12[t] - 
  Sqrt[2]*G1*a1*V13[t] - Sqrt[2]*G1*a2*V14[t] == 0, 
V13'[t] - V23[t]*w1 + k1*V13[t] + 
  Sqrt[2]*G1*a2*V11[t] - (-G1*Q1 + delc)*V14[t] == 0, 
V14'[t] - V24[t]*w1 + k1*V14[t] - 
  Sqrt[2]*G1*a1*V11[t] - (G1*Q1 - delc)*V13[t] == 0, 
V21'[t] + V11[t]*w1 + gma1*V21[t] - Sqrt[2]*G1*a1*V31[t] - 
  Sqrt[2]*G1*a2*V41[t] - w1*V22[t] == 0, 
V22'[t] + V12[t]*w1 + gma1*V22[t] - Sqrt[2]*G1*a1*V32[t] - 
  Sqrt[2]*G1*a2*V42[t] + w1*V21[t] + gma1*V22[t] - 
  Sqrt[2]*G1*a1*V23[t] - Sqrt[2]*G1*a2*V24[t] - gma2*(2*n1 + 1) ==
  0, V23'[t] + w1*V13[t] + gma1*V23[t] - Sqrt[2]*G1*a1*V33[t] - 
  Sqrt[2]*G1*a2*V43[t] + k1*V23[t] + 
  Sqrt[2]*G1*a2*V21[t] - (-G1*Q1 + delc)*V24[t] == 0,
V24'[t] + V14[t]*w1 + gma1*V24[t] - Sqrt[2]*G1*a1*V34[t] - 
  Sqrt[2]*G1*a2*V44[t] + k1*V24[t] - 
  Sqrt[2]*G1*a1*V21[t] - (G1*Q1 - delc)*V23[t] == 0, 
V31'[t] + k1*V31[t] + 
  Sqrt[2]*G1*a2*V11[t] - (-G1*Q1 + delc)*V41[t] - w1*V32[t] == 0, 
V32'[t] + k1*V32[t] + 
  Sqrt[2]*G1*a2*V12[t] - (-G1*Q1 + delc)*V42[t] + w1*V31[t] - 
  Sqrt[2]*G1*a2*V34[t] - Sqrt[2]*G1*a1*V33[t] + gma1*V32[t] == 0, 
V33'[t] + k1*V33[t] + 
  Sqrt[2]*G1*a2*V13[t] - (-G1*Q1 + delc)*V43[t] + k1*V33[t] + 
  Sqrt[2]*G1*a2*V31[t] - (-G1*Q1 + delc)*V34[t] - k0 == 0, 
V34'[t] + k1*V34[t] + 
  Sqrt[2]*G1*a2*V14[t] - (-G1*Q1 + delc)*V44[t] + k1*V34[t] - 
  Sqrt[2]*G1*a1*V31[t] - (G1*Q1 - delc)*V33[t] == 0, 
V41'[t] + k1*V41[t] - 
  Sqrt[2]*G1*a1*V11[t] - (G1*Q1 - delc)*V31[t] - w1*V42[t] == 0, 
V42'[t] + k1*V42[t] + 
  Sqrt[2]*G1*a1*V12[t] - (G1*Q1 - delc)*V32[t] + w1*V41[t] - 
  Sqrt[2]*G1*a2*V44[t] - Sqrt[2]*G1*a1*V43[t] + gma1*V42[t] == 0, 
V43'[t] + k1*V43[t] - 
  Sqrt[2]*G1*a1*V13[t] - (G1*Q1 - delc)*V33[t] + k1*V43[t] + 
  Sqrt[2]*G1*a2*V41[t] - (-G1*Q1 + delc)*V44[t] == 0, 
V44'[t] + k1*V44[t] - 
  Sqrt[2]*G1*a1*V14[t] - (G1*Q1 - delc)*V34[t] + k1*V44[t] - 
  Sqrt[2]*G1*a1*V41[t] - (G1*Q1 - delc)*V43[t] - k0 == 0, 
V11[0] == 1, V12[0] == 1, V13[0] == 0, V14[0] == 0, V21[0] == 0, 
V22[0] == 1, V23[0] == 0, V24[0] == 0, V31[0] == 0, V32[0] == 0, 
V33[0] == 0, V34[0] == 0, V41[0] == 0, V42[0] == 0, V43[0] == 0, 
V44[0] == 0}, {V11, V12, V13, V14, V21, V22, V23, V24, V31, V32, 
V33, V34, V41, V42, V43, V44},{t, 0, 100},g0];
P2 = Plot[{Evaluate[1/2*(V11[t] + V22[t] - 2*V12[t])^(-1) /. s]}, {t, 
   0, 60}, PlotRange -> {0, 1}, Frame -> True,    
FrameLabel -> {Style["Time", Bold, 20], 
Style[" \!\(\*SubscriptBox[\(S\), \(q\)]\)", Bold, 20]}, 
 FrameTicksStyle -> Directive[FontSize -> 20], 
 PlotStyle -> {Thickness[0.0005], Thickness[0.008]}]
$\endgroup$
4
  • $\begingroup$ Can't find $g_0$ in your code for the DEs. $\endgroup$ – LouisB Jan 21 '20 at 11:34
  • $\begingroup$ g0 is the parameter inside ParametricNDSolveValue . $\endgroup$ – Ulrich Neumann Jan 21 '20 at 11:35
  • $\begingroup$ It's unclear what you want to plot "about g0". P2=... is a time plot (which doesn't work) $\endgroup$ – Ulrich Neumann Jan 21 '20 at 11:37
  • $\begingroup$ @LouisB , sorry its not g0 , it is delc. Now please again see the text. $\endgroup$ – vini Jan 21 '20 at 11:46
3
$\begingroup$

To plot the solution for a range of delc, make these changes to your code:

(1) Remove delc = 1; Leave delc undefined. It will be your parameter.

(2) Keep the ParametricDSolveValue command and its first argument, but change the other arguments to get

s = ParametricNDSolveValue[  ...   ,
    {1/2*(V11[t] + V22[t] - 2*V12[t])^(-1)}, {t, 0, 100}, delc];

Notice that the second argument is just the expression that we want to plot. This is not the only way to do it, but it has its advantages. Evaluate the code and plot the solution over the range of delc at t = 60 like this

Plot[s[delc] /. t -> 60, {delc, 0, 2},
 Frame -> True, PlotRange -> {All, {0.545, .5475}},
 FrameLabel -> {Style["\[CapitalDelta]c", Bold, 20], 
   Style["\!\(\*SubscriptBox[\(S\), \(q\)]\)", Bold, 20],
   Style["t = 60", Bold, 20]},
 FrameTicksStyle -> Directive[FontSize -> 20],
 GridLines->Automatic]

Figure 1

We can also plot the solution versus time for delc = 1.25 like this

Plot[s[1.25], {t, 0, 100},
 PlotRange -> {All, {0, 15}}, Frame -> True,
 FrameLabel -> {Style["t", Bold, 20], 
   Style["\!\(\*SubscriptBox[\(S\), \(q\)]\)", Bold, 20],
   Style["\[CapitalDelta]c = 1.25", Bold, 20]},
 FrameTicksStyle -> Directive[FontSize -> 20],
 GridLines -> Automatic]

Figure 2

Edit to Plot an Average

One way to average $S_q(\Delta c)$ for $\Delta c$ range from 0 to 2 is to simply evaluate s[delc] at various values of $\Delta c$ and find their mean. Here is the (complete) code that accomplishes that:

ClearAll["Global`*"]
w1 = 1;
gma1 = 0.005;
n1 = 1;
gma2 = 0.005;
G1 = 0.005;

k1 = .1;
k2 = 0.1;
a1 = 0.07;
a2 = 0.58;
k0 = 0.1;
Q1 = 1.268;
del0 = 1;
N1 = 1;
ome = 1;
M1 = del0*(1 - Cos[ome*t]);

s = ParametricNDSolveValue[{
    V11'[t] - V21[t]*w1 - V12[t]*w1 == 0,
    V12'[t] - V22[t]*w1 + w1*V11[t] + gma1*V12[t] -
      Sqrt[2]*G1*a1*V13[t] - Sqrt[2]*G1*a2*V14[t] == 0,
    V13'[t] - V23[t]*w1 + k1*V13[t] +
      Sqrt[2]*G1*a2*V11[t] - (-G1*Q1 + delc)*V14[t] == 0,
    V14'[t] - V24[t]*w1 + k1*V14[t] -
      Sqrt[2]*G1*a1*V11[t] - (G1*Q1 - delc)*V13[t] == 0,
    V21'[t] + V11[t]*w1 + gma1*V21[t] -
      Sqrt[2]*G1*a1*V31[t] - Sqrt[2]*G1*a2*V41[t] - w1*V22[t] == 0,
    V22'[t] + V12[t]*w1 + gma1*V22[t] -
      Sqrt[2]*G1*a1*V32[t] - Sqrt[2]*G1*a2*V42[t] +
      w1*V21[t] + gma1*V22[t] - Sqrt[2]*G1*a1*V23[t] -
      Sqrt[2]*G1*a2*V24[t] - gma2*(2*n1 + 1) == 0,
    V23'[t] + w1*V13[t] + gma1*V23[t] -
      Sqrt[2]*G1*a1*V33[t] - Sqrt[2]*G1*a2*V43[t] + k1*V23[t] +
      Sqrt[2]*G1*a2*V21[t] - (-G1*Q1 + delc)*V24[t] == 0, 
    V24'[t] + V14[t]*w1 + gma1*V24[t] - Sqrt[2]*G1*a1*V34[t] -
      Sqrt[2]*G1*a2*V44[t] + k1*V24[t] -
      Sqrt[2]*G1*a1*V21[t] - (G1*Q1 - delc)*V23[t] == 0,
    V31'[t] + k1*V31[t] + Sqrt[2]*G1*a2*V11[t] -
      (-G1*Q1 + delc)*V41[t] - w1*V32[t] == 0,
    V32'[t] + k1*V32[t] + Sqrt[2]*G1*a2*V12[t] -
      (-G1*Q1 + delc)*V42[t] + w1*V31[t] -
      Sqrt[2]*G1*a2*V34[t] - Sqrt[2]*G1*a1*V33[t] + gma1*V32[t] == 0,
    V33'[t] + k1*V33[t] + Sqrt[2]*G1*a2*V13[t] -
      (-G1*Q1 + delc)*V43[t] + k1*V33[t] +
      Sqrt[2]*G1*a2*V31[t] - (-G1*Q1 + delc)*V34[t] - k0 == 0,
    V34'[t] + k1*V34[t] + Sqrt[2]*G1*a2*V14[t] -
      (-G1*Q1 + delc)*V44[t] + k1*V34[t] -
      Sqrt[2]*G1*a1*V31[t] - (G1*Q1 - delc)*V33[t] == 0,
    V41'[t] + k1*V41[t] - Sqrt[2]*G1*a1*V11[t] -
      (G1*Q1 - delc)*V31[t] - w1*V42[t] == 0,
    V42'[t] + k1*V42[t] + Sqrt[2]*G1*a1*V12[t] -
      (G1*Q1 - delc)*V32[t] + w1*V41[t] -
      Sqrt[2]*G1*a2*V44[t] - Sqrt[2]*G1*a1*V43[t] + gma1*V42[t] == 0,
    V43'[t] + k1*V43[t] - Sqrt[2]*G1*a1*V13[t] -
      (G1*Q1 - delc)*V33[t] + k1*V43[t] +
      Sqrt[2]*G1*a2*V41[t] - (-G1*Q1 + delc)*V44[t] == 0,
    V44'[t] + k1*V44[t] - Sqrt[2]*G1*a1*V14[t] -
      (G1*Q1 - delc)*V34[t] + k1*V44[t] -
      Sqrt[2]*G1*a1*V41[t] - (G1*Q1 - delc)*V43[t] - k0 == 0,
    V11[0] == 1, V12[0] == 1, V13[0] == 0, V14[0] == 0,
    V21[0] == 0, V22[0] == 1, V23[0] == 0, V24[0] == 0,
    V31[0] == 0, V32[0] == 0, V33[0] == 0, V34[0] == 0,
    V41[0] == 0, V42[0] == 0, V43[0] == 0, V44[0] == 0},
   1/2*(V11[t] + V22[t] - 2*V12[t])^(-1), {t, 0, 100}, delc];

With[{avg = Mean@Table[s[dc], {dc, Subdivide[0, 2, 10]}]},
 Plot[avg, {t, 0, 100}, PlotRange -> {0, 15}]]

enter image description here

Note that Subdivide[0, 2, 10] gives 11 evenly spaced sample points. You can generate that same list using Range[0, 2, 2/10] or any other method you choose. Eleven is enough points for the plot. It turns out s[delc] gives practically the same function of $t$ for any value of delc in the interval from 0 to 2.

$\endgroup$
17
  • $\begingroup$ Please clarify 2nd point you mentioned above because its still not giving any plot. $\endgroup$ – vini Jan 22 '20 at 6:44
  • $\begingroup$ Here's what I did: Start with a fresh kernel. Copy / paste the code from the question into a new notebook. Remove the line containing delc = 1;. Remove all the code beginning at {V11, V12, V13 .. to the end. Copy / paste the second line of code from item (2) in answer. Copy / paste the plot command from the answer. Runs fine for me. If you can't get it to work, maybe you should edit your question to show exactly what code you are using and what results or error message you are getting. $\endgroup$ – LouisB Jan 22 '20 at 6:59
  • $\begingroup$ Thanks very much. its working now. $\endgroup$ – vini Jan 22 '20 at 7:07
  • $\begingroup$ This plot that you have plotted, is for a particular value of time.i.e t=60, can you tell me how can i plot for t=30,10,20,70,.... and after that average of all plots? $\endgroup$ – vini Jan 30 '20 at 14:34
  • $\begingroup$ Plot[s[delc] /. t -> {10, 20, 30, 70}, {delc, 0, 2}] gives a quick plot at the 4 times. The "curves" look like straight lines because of the scale on the vertical axis. I don't think you really want to average these curves. Please consider asking a new question and tell us what you are trying to achieve, show or study, so we can provide answers that better hit the target. You could average those curves, however, with Plot[Mean[s[delc] /. t -> {10, 20, 30, 70}], {delc, 0, 2}] $\endgroup$ – LouisB Jan 30 '20 at 22:21

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.