Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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}}]
share|improve this question
1  
Please have a look at this. –  b.gatessucks Mar 12 at 10:17
    
Also this Wolfram training video has a bit on FunctionInterpolation that is very relevant –  gpap Mar 12 at 10:38

2 Answers 2

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

share|improve this answer
    
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. –  PhysicistRRE Mar 12 at 12:56

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
    }
  ]
share|improve this answer
    
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. –  PhysicistRRE Mar 12 at 12:38
    
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] ." –  gpap Mar 12 at 14:16
    
Noted, I was a bit unclear. I've rephrased the question. –  PhysicistRRE Mar 12 at 20:35

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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