0
$\begingroup$

This is an ODE, I want to plot P on x-axis and vars on y-axis. My code works fine but gives an empty plot

Subscript[V, 0] = -2.5;
Subscript[k, 0] = Pi/2;
\[Epsilon] = 0.05;
n = 10;
B = Sqrt[3];
Subscript[t, fin] = 25;

A = Table[Subscript[V, i], {i, n}];
For[j = 1, j < n + 1, j++, A[[j]] = 0];
A[[n/2]] = Subscript[V, 0]*(1 + \[Epsilon]); 
A[[n/2 + 1]] = Subscript[V, 0]*(1 - \[Epsilon]);
P = Table[Subscript[\[Alpha], i], {i, n}];
For[j = 1, j < n + 1, j++, P[[j]] = 0];
P[[n/2]] = 1; P[[n/2 + 1]] = 1;
Subscript[\[Phi], 0][t_] := Subscript[\[Phi], 1][t];
Subscript[\[Phi], n + 1][t_] := Subscript[\[Phi], n][t];

eqns = Table[{Sqrt[-1]*Subscript[\[Phi], i]'[t] == 
 A[[i]]*Subscript[\[Phi], i][t] - Subscript[\[Phi], i + 1][t] - 
  Subscript[\[Phi], i - 1][
   t] + (P[[i]]*Abs[Subscript[\[Phi], i][t]]^2*
    Subscript[\[Phi], i][t]), Subscript[\[Phi], i][0] == B}, {i, 
n}];
vars = Table[
Sum[Abs[\[Phi][i, Subscript[t, fin]]]^2, {i, n/2 + 2, n}], {i, 
25}]/Table[Sum[Abs[\[Phi][i, 0]]^2, {i, n/2 - 1}], {i, 25}];

sol = NDSolve[eqns, vars, {t, 0, 25}, 
Method -> {"ExplicitRungeKutta", "StiffnessTest" -> False}, 
MaxSteps -> \[Infinity], AccuracyGoal -> 8, PrecisionGoal -> 8];
Plot=ListPlot[Table[vars /. First[%], {t, 0, 250}], 
PlotRange -> All, ImageSize -> 400]

with vars explicitly given by vars$=\frac{\sum_{n>\frac{M}{2}+1}|\phi_n(t_{fin})|^2}{\sum_{n<\frac{M}{2}}|\phi_n(0)|^2}$

$\endgroup$
12
  • $\begingroup$ Check vars beceause is nonsense ? $\endgroup$ Jul 27, 2018 at 10:33
  • $\begingroup$ No, its the transmission coefficient. So essentially it will be a plot of transmission coefficient as a function of P. did I make a mistake in writing it appropriately you mean? $\endgroup$
    – AtoZ
    Jul 27, 2018 at 10:37
  • 1
    $\begingroup$ Replace by my code: vars = Table[Subscript[\[Phi], i][t], {i, 0, n}];Plot[Evaluate[vars /. sol // ReIm], {t, 0, 25}, PlotRange -> All, ImageSize -> 400, PlotLegends -> Table[\[Phi][i][t], {i, 0, n}]]? $\endgroup$ Jul 27, 2018 at 10:49
  • $\begingroup$ vars=$\frac{\sum_{n>\frac{M}{2}+1}|\phi_n(t_{fin})|^2}{\sum_{n<\frac{M}{2}}|\phi_n(0)|^2}$ $\endgroup$
    – AtoZ
    Jul 27, 2018 at 10:54
  • 1
    $\begingroup$ If I check your equation a see a lot of subscripted functions like Subscript[\[Phi], 6][t] were as your variables look like (Abs[\[Phi][7, 25]]^2) +.... Your main problem is due to the use of subscripted variables! $\endgroup$ Jul 27, 2018 at 13:03

1 Answer 1

1
$\begingroup$

It is possible that this code will be useful

Subscript[V, 0] = -2.5;
Subscript[k, 0] = Pi/2;
\[Epsilon] = 0.05;
n = 10;
B = Sqrt[3];
(*Subscript[t,fin]=25;*)

A = Table[Subscript[V, i], {i, n}];
For[j = 1, j < n + 1, j++, A[[j]] = 0];
A[[n/2]] = Subscript[V, 0]*(1 + \[Epsilon]);
A[[n/2 + 1]] = Subscript[V, 0]*(1 - \[Epsilon]);
P = Table[Subscript[\[Alpha], i], {i, n}];
For[j = 1, j < n + 1, j++, P[[j]] = 0];
P[[n/2]] = 1; P[[n/2 + 1]] = 1;
Subscript[\[Phi], 0][t_] := Subscript[\[Phi], 1][t];
Subscript[\[Phi], n + 1][t_] := Subscript[\[Phi], n][t];

eqns = Table[{Sqrt[-1]*Subscript[\[Phi], i]'[t] == 
     A[[i]]*Subscript[\[Phi], i][t] - Subscript[\[Phi], i + 1][t] - 
      Subscript[\[Phi], i - 1][
       t] + (P[[i]]*Abs[Subscript[\[Phi], i][t]]^2*
        Subscript[\[Phi], i][t]), Subscript[\[Phi], i][0] == B}, {i, 
    n}];


sol = NDSolve[eqns, 
   Table[Subscript[\[Phi], i][t], {i, 1, n}], {t, 0, 25}, 
   Method -> {"ExplicitRungeKutta", "StiffnessTest" -> False}, 
   MaxSteps -> \[Infinity], AccuracyGoal -> 8, PrecisionGoal -> 8];
f[i_, t_] := Abs[Subscript[\[Phi], i][t]]^2 /. sol
A = Sum[f[i, t] /. t -> 0, {i, 1, n/2 - 1}];
Plot[Sum[f[i, t], {i, n/2 + 2, n}]/A, {t, 0, 25}, 
 PlotLabel -> "vars(t)"]

Table[Plot[
  Evaluate[Abs[Subscript[\[Phi], i][t]]^2 /. sol], {t, 0, 25}, 
  PlotLabel -> i, AxesLabel -> Automatic], {i, 1, n}]

fig1

If there is a dependence on the parameter, then use the following code

Subscript[V, 0] = -2.5;
Subscript[k, 0] = Pi/2;
\[Epsilon] = 0.05;
n = 10;
B = Sqrt[3];
Subscript[t, fin] = 25;

A = Table[Subscript[V, i], {i, n}];
For[j = 1, j < n + 1, j++, A[[j]] = 0];
A[[n/2]] = Subscript[V, 0]*(1 + \[Epsilon]);
A[[n/2 + 1]] = Subscript[V, 0]*(1 - \[Epsilon]);
P = Table[Subscript[\[Alpha], i], {i, n}];
For[j = 1, j < n + 1, j++, P[[j]] = 0];
P[[n/2]] = p0; P[[n/2 + 1]] = p0;
Subscript[\[Phi], 0][t_] := Subscript[\[Phi], 1][t];
Subscript[\[Phi], n + 1][t_] := Subscript[\[Phi], n][t];

eqns = Table[{Sqrt[-1]*Subscript[\[Phi], i]'[t] == 
     A[[i]]*Subscript[\[Phi], i][t] - Subscript[\[Phi], i + 1][t] - 
      Subscript[\[Phi], i - 1][
       t] + (P[[i]]*Abs[Subscript[\[Phi], i][t]]^2*
        Subscript[\[Phi], i][t]), Subscript[\[Phi], i][0] == B}, {i, 
    n}];


F = Table[
   ParametricNDSolveValue[eqns, 
    Subscript[\[Phi], i], {t, 0, 25}, {p0}], {i, 1, n}];



{Plot3D[Evaluate[Abs[F[[n/2]][p0][t]]^2], {p0, 0, 10}, {t, 0, 25}, 
  Mesh -> None, ColorFunction -> Hue, 
  PlotLabel -> "\!\(\*SubscriptBox[\(\[Phi]\), \(5\)]\)(p0,t)", 
  AxesLabel -> {"p0", "t", ""}, PlotRange -> All], 
 Plot3D[Evaluate[
   Sum[Abs[F[[i]][p0][t]]^2, {i, n/2 + 2, n}]/
    Sum[Abs[F[[i]][p0][0]]^2, {i, 1, n/2 - 1}]], {p0, 0, 10}, {t, 0, 
   25}, Mesh -> None, ColorFunction -> Hue, PlotLabel -> "vars(p0,t)",
   AxesLabel -> {"p0", "t", ""}, PlotRange -> All]}

fig2

$\endgroup$
11
  • $\begingroup$ Thanks. Yes this is very helpful. Is it possible to plot vars as function of P, instead of vars(t)? $\endgroup$
    – AtoZ
    Jul 28, 2018 at 12:08
  • 1
    $\begingroup$ This is possible, but by definition P = 0. There must be parameters or at least one parameter, eg For[j = 1, j < n + 1, j++, P[[j]] = p0]; $\endgroup$ Jul 28, 2018 at 12:34
  • 1
    $\begingroup$ Ok, P= {0, 0, 0, 0, 1, 1, 0, 0, 0, 0}. Where is the parameter? $\endgroup$ Jul 28, 2018 at 12:53
  • 1
    $\begingroup$ Plot3D[vars[p0,t],{p0,10},{t,0,25}] ? $\endgroup$ Jul 28, 2018 at 13:14
  • 1
    $\begingroup$ I updated the solution, added 3D figure. It's probably what you wanted to see. $\endgroup$ Jul 29, 2018 at 4:58

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.