2
$\begingroup$

Code:

(*defining constants*)

neff = 1.98;
a        = 1/3;
λ        = 1.55;
β        = (2 π*neff)/λ + I*a/2;
k21    = I;
k12    = I;

R1[θ] = a/(2 π)*θ; 
R2[θ] = a + a/(2 π)*θ;
(*defining the two ODEs*)
ode1 = u1'[θ] == 
   I*β*R1[θ]*u1[θ] + k21 * u2[θ];
ode2 = u2'[θ] == 
   I*β*R2[θ]*u2[θ] + k12 * u1[θ];

(*numerically solve via Runge-Kutta*)
sol = NDSolve[
  
  {ode1, ode2, 
   u1[0] == 1, 
   u2[0] == E^(I*β*2 Sqrt[1 + 1/(36 π^2)] π)},  
  
  {u1[θ], u2[θ]}, {θ, 0, 2 π}, 
  Method -> "ExplicitRungeKutta", "StartingStepSize" -> 1/10]

I try to plot the output, but fail:

Plot[Evaluate[{u1[θ], u2[θ]} /. sol], {θ, 0, 
  2 π}, Frame -> True]

as this plots the frame of the graph, but not $u_1$ or $u_2$. Why?

$\endgroup$
6
  • 1
    $\begingroup$ The values of the functions you are trying to plot are complex. Whenever you have trouble plotting, try looking at the values resulting from your calculation, for instance by tabulating a few samples: Table[{u1[\[Theta]], u2[\[Theta]]} /. sol, {\[Theta], 0, 2 Pi, Pi/2}]. $\endgroup$
    – MarcoB
    Jul 28 at 21:43
  • $\begingroup$ I wasn't aware we could do that tabulating, Thank you very much. I'll check how to plot the real/imaginary $\endgroup$
    – Nalt
    Jul 28 at 21:45
  • $\begingroup$ Try replacing your Evaluate expression in Plot with Evaluate[ReIm@{u1[\[Theta]], u2[\[Theta]]} /. sol]. $\endgroup$
    – MarcoB
    Jul 28 at 21:46
  • $\begingroup$ There's also ReImPlot[]. $\endgroup$
    – Michael E2
    Jul 28 at 23:04
  • 1
    $\begingroup$ See also the "complex-valued nonlinear reaction equation" here: reference.wolfram.com/language/ref/NDSolve.html#1648865610 $\endgroup$
    – Michael E2
    Jul 28 at 23:09

1 Answer 1

2
$\begingroup$
Clear["Global`*"]

neff = 198/100;
a = 1/3;
λ = 155/100;
β = (2 π*neff)/λ + I*a/2;
k21 = I;
k12 = I;

R1[θ] = a/(2 π)*θ;
R2[θ] = a + a/(2 π)*θ;

ode1 = u1'[θ] == I*β*R1[θ]*u1[θ] + k21*u2[θ];
ode2 = u2'[θ] == 
   I*β*R2[θ]*u2[θ] + k12*u1[θ];

eqns = {ode1, ode2, u1[0] == 1, 
    u2[0] == E^(I*β*2 Sqrt[1 + 1/(36 π^2)] π)} // Simplify;

The system can be solved exactly using DSolve:

sol = DSolve[{ode1, ode2, u1[0] == 1, 
     u2[0] == E^(I*β*2 Sqrt[1 + 1/(36 π^2)] π)}, {u1, 
     u2}, θ][[1]];

Verifying the solution,

And @@ (eqns /. sol // Simplify)

(* True *)

Plot[
 Evaluate[Simplify[ReIm[{u1[θ], u2[θ]} /. sol]]],
 {θ, 0, 2 π},
 Frame -> True,
 PlotStyle -> {Automatic, Dashed, Automatic, Dotted},
 PlotLegends -> {Re[u1], Im[u1], Re[u2], Im[u2]}]

enter image description here

$\endgroup$

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