1
$\begingroup$

I want to plot the solution a coupled ODEs as a function of the end point only. Mathematica code:

tmax = 5;
rr = 1;
p = 0;
x[r_, t_, p_] = r Cos[(2 \[Pi])/tmax t + p];
z[r_, t_, p_] = r Sin[(2 \[Pi])/tmax t + p];
g[r_, t_, 
   p_] = ((1 + z[r, t, p]) D[x[r, t, p], t] - (I + x[r, t, p]) D[
     z[r, t, p], t])/(2 ((1 + z[r, t, p]^2) + (I + x[r, t, p]^2)));
b[r_, t_] = Sqrt[(1 + I) + r^2] t;
sol1 = NDSolve[{c1'[t] == -g[rr, t, p] Exp[-2 I b[rr, t] ] c2[t], 
    c2'[t] == g[rr, t, p] Exp[2 I b[rr, t]] c1[t], c1[0] == 1, 
    c2[0] == 0}, {c1, c2}, {t, 0, tmax}];
sol2 = NDSolve[{c1'[t] == -g[rr, t, p] Exp[-2 I b[rr, t] ] c2[t], 
    c2'[t] == g[rr, t, p] Exp[2 I b[rr, t]] c1[t], c1[0] == 0, 
    c2[0] == 1}, {c1, c2}, {t, 0, tmax}];

Then want to repeat this procedure for different tmax and plot LogPlot[Evaluate[{{Abs[c2[t]/c1[t]]} /. sol1, {Abs[c1[t]/c2[t]]} /. sol2}], {t, 0, tmax}, PlotRange -> Full]

Many thanks to the helpers

EDIT:

In block diagram this is what I want:

Given coupled ODEs that are dependent on tmax both as the end-point and as a parameter, find the solution for each tmax, c1_tmax[t] c2_tmax[t]

define the functions:

c1[tmax]=c1_tmax[tmax]

c2[tmax]=c2_tmax[tmax]

plot c2[tmax]/c1[tmax] for {tmax,0,20}

Here's what I finally did (I hope it actually does it):

rr = 0.8;
p = 0;
x[r_, t_, p_] = r Cos[(2 \[Pi])/tmax t + p];
z[r_, t_, p_] = r Sin[(2 \[Pi])/tmax t + p];
g[r_, t_, 
   p_] = ((1 + z[r, t, p]) D[x[r, t, p], t] - (I + x[r, t, p]) D[
     z[r, t, p], t])/(2 ((1 + z[r, t, p]^2) + (I + x[r, t, p]^2)));
b[r_, t_] = Sqrt[(1 + I) + r^2] t;
sol1 = ParametricNDSolve[{c1'[
      t] == -g[rr, t, p] Exp[-2 I b[rr, t] ] c2[t], 
    c2'[t] == g[rr, t, p] Exp[2 I b[rr, t]] c1[t], c1[0] == 1, 
    c2[0] == 0}, {c1, c2}, {t, 0, tmax}, {tmax}];
sol2 = ParametricNDSolve[{c1'[
      t] == -g[rr, t, p] Exp[-2 I b[rr, t] ] c2[t], 
    c2'[t] == g[rr, t, p] Exp[2 I b[rr, t]] c1[t], c1[0] == 0, 
    c2[0] == 1}, {c1, c2}, {t, 0, tmax}, {tmax}];
LogPlot[Evaluate[{{Abs[c2[tmax][tmax]/c1[tmax][tmax]]} /. 
    sol1, {Abs[c1[tmax][tmax]/c2[tmax][tmax]]} /. sol2}], {tmax, 0, 
  20}, PlotRange -> {{0, 20}, {10^-3, 10^3}}]
$\endgroup$
1
$\begingroup$

Ok, so as mentioned in the comments, you can go either of two ways. With your definitions (note I omitted tmax):

rr = 1;
p = 0;
x[r_, t_, p_] = r Cos[(2 π)/tmax t + p];
z[r_, t_, p_] = r Sin[(2 π)/tmax t + p];
g[r_, t_, p_] = ((1 + z[r, t, p]) D[x[r, t, p], t] - (I + x[r, t, p]) D[
       z[r, t, p], t])/(2 ((1 + z[r, t, p]^2) + (I + x[r, t, p]^2)));
b[r_, t_] = Sqrt[(1 + I) + r^2] t;

Using NDSolve

sol1[tmax_] /; tmax > 0 := 
  Evaluate@NDSolve[{c1'[t] == -g[rr, t, p] Exp[-2 I b[rr, t]] c2[t], 
     c2'[t] == g[rr, t, p] Exp[2 I b[rr, t]] c1[t], c1[0] == 1, 
     c2[0] == 0}, {c1, c2}, {t, 0, tmax}];
sol2[tmax_] /; tmax > 0 := 
  Evaluate@NDSolve[{c1'[t] == -g[rr, t, p] Exp[-2 I b[rr, t]] c2[t], 
     c2'[t] == g[rr, t, p] Exp[2 I b[rr, t]] c1[t], c1[0] == 0, 
     c2[0] == 1}, {c1, c2}, {t, 0, tmax}];

you can now solve for every initial condition

f1[t_, tmax_] := (Abs[c2[t]/c1[t]] /. Flatten@sol1[tmax]);
f2[t_, tmax_] := (Abs[c2[t]/c1[t]] /. Flatten@sol2[tmax]);

and create a plot for a range of tmax from 1 to 5:

AbsoluteTiming[Table[
    LogPlot[
     Evaluate[{f1[t, tmax], f2[t, tmax]}], {t, 0, tmax},
     PlotRange -> {{0, 5}, All},
     PlotLabel -> "tmax = " <> ToString@tmax], 
    {tmax, 0.5, 5, .1}] // Export["test.gif", #] &;]

{5.88, Null}

and here's how it looks:

enter image description here

Using ParametricNDSolve

sol11 = ParametricNDSolve[{c1'[
      t] == -g[rr, t, p] Exp[-2 I b[rr, t]] c2[t], 
    c2'[t] == g[rr, t, p] Exp[2 I b[rr, t]] c1[t], c1[0] == 1, 
    c2[0] == 0}, {c1, c2}, {t, 0, tmax}, {tmax}];
sol12 = ParametricNDSolve[{c1'[
      t] == -g[rr, t, p] Exp[-2 I b[rr, t]] c2[t], 
    c2'[t] == g[rr, t, p] Exp[2 I b[rr, t]] c1[t], c1[0] == 0, 
    c2[0] == 1}, {c1, c2}, {t, 0, tmax}, {tmax}];

now you define the functions somewhat differently:

g1[t_, tmax_] := Abs[Evaluate[c2[tmax][t]/c1[tmax][t]]] /. sol11;
g2[t_, tmax_] := Abs[Evaluate[c2[tmax][t]/c1[tmax][t]]] /. sol12;

and again do the same:

AbsoluteTiming[Table[
    LogPlot[
     Evaluate[{g1[t, tmax], g2[t, tmax]}], {t, 0, tmax},
     PlotRange -> {{0, 5}, All},
     PlotLabel -> "tmax = " <> ToString@tmax], 
    {tmax, 0.5, 5, .1}] // Export["test2.gif", #] &;]

{6.17, Null}

and you get the same plot:

enter image description here

$\endgroup$
  • $\begingroup$ Many thanks,but this plot the solution upto tmax whereas I wanted the solutions as function of tmax. I have clarified my request (and proposed solution) in the original question. $\endgroup$ – PhysicistRRE Mar 12 '14 at 12:56
0
$\begingroup$

Is this what you want ?

 Manipulate[
  tmax = endpoint;

   sol1 = NDSolve[{c1'[t] == -g[rr, t, p] Exp[-2 I b[rr, t]] c2[t], 
     c2'[t] == g[rr, t, p] Exp[2 I b[rr, t]] c1[t], c1[0] == 1, 
     c2[0] == 0}, {c1, c2}, {t, 0, tmax}];
  sol2 = NDSolve[{c1'[t] == -g[rr, t, p] Exp[-2 I b[rr, t]] c2[t], 
     c2'[t] == g[rr, t, p] Exp[2 I b[rr, t]] c1[t], c1[0] == 0, 
     c2[0] == 1}, {c1, c2}, {t, 0, tmax}];

  LogPlot[Evaluate[{{Abs[c2[t]/c1[t]]} /. sol1, {Abs[c1[t]/c2[t]]} /.sol2}],
   {t, 0, tmax},PlotRange -> Full],

 {{endpoint, 5}, $MachineEpsilon, 50, Appearance -> "Labeled"},

  Initialization :> {rr = 1,
    p = 0,
    x[r_, t_, p_] = r Cos[(2 \[Pi])/tmax t + p],
    z[r_, t_, p_] = r Sin[(2 \[Pi])/tmax t + p],
    g[r_, t_, p_] = ((1 + z[r, t, p]) D[x[r, t, p], t] - (I + x[r, t, p]) D[
          z[r, t, p], t])/(2 ((1 + z[r, t, p]^2) + (I + x[r, t, p]^2))),
    b[r_, t_] = Sqrt[(1 + I) + r^2] t
    }
  ]
$\endgroup$
  • $\begingroup$ Thank, but not quite. I want to solves the coupled ODEs which are dependent on tmax both in the equations themselves and in the end point of the integration. Then I'd like to take the different solution to different tmax's and present it as a function of tmax only. $\endgroup$ – PhysicistRRE Mar 12 '14 at 12:38
  • $\begingroup$ That's not what you said in your original post though:"I want to repeat this procedure for different tmax and plot LogPlot[CODE, {t, 0, tmax}, PlotRange -> Full] ." $\endgroup$ – gpap Mar 12 '14 at 14:16
  • $\begingroup$ Noted, I was a bit unclear. I've rephrased the question. $\endgroup$ – PhysicistRRE Mar 12 '14 at 20: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.