0
$\begingroup$

I am trying to plot Mathieu Functions, running for an hour without a result. Here is my code:

  1. First define the relevant physical constant:
Needs["Notation`"]
Symbolize[ParsedBoxWrapper[SubscriptBox["v", "F"]]]; Symbolize[ParsedBoxWrapper[SubscriptBox["g", "2"]]]; 
  Symbolize[ParsedBoxWrapper[SubscriptBox["q", "c"]]]; 
  1. Define all the relevant functions:
g\[UnderBracket]Subscript\[UnderBracket]2[\[Nu]_, q_, t_] := (v\[UnderBracket]Subscript\[UnderBracket]F*Exp[(-(1/2))*(q/q\[UnderBracket]Subscript\[UnderBracket]c)^2]*(1 - Cos[(\[Nu]*Pi*t)/\[Tau]]))/
   (4*(2*Pi*v\[UnderBracket]Subscript\[UnderBracket]F))
w[\[Nu]_, q_] := Sqrt[(1 + g\[UnderBracket]Subscript\[UnderBracket]2[\[Nu], q, \[Tau]])^2 - g\[UnderBracket]Subscript\[UnderBracket]2[\[Nu], q, \[Tau]]^2]
s[\[Nu]_, q_] := Sqrt[(1/2)*((1 + g\[UnderBracket]Subscript\[UnderBracket]2[\[Nu], q, \[Tau]])/w[\[Nu], q] - 1)]
c[\[Nu]_, q_] := Sqrt[1 + s[\[Nu], q]^2]
y[\[Nu]_, q_, t_] = FullSimplify[
    a[t] /. First[DSolve[{Derivative[2][a][t] + v\[UnderBracket]Subscript\[UnderBracket]F^2*q^2*(1 + g\[UnderBracket]Subscript\[UnderBracket]2[\[Nu], q, t])*a[t] == 
         0, a[0] == 1, Derivative[1][a][0] == (-I)*v\[UnderBracket]Subscript\[UnderBracket]F*q}, a[t], t]]]; 
v[\[Nu]_, q_, t_] := FullSimplify[(1/2)*y[\[Nu], q, t] - (I*D[y[\[Nu], q, t], t])/(2*v\[UnderBracket]Subscript\[UnderBracket]F*q)]
u[\[Nu]_, q_, t_] := FullSimplify[(1/2)*y[\[Nu], q, t] + (I*D[y[\[Nu], q, t], t])/(2*v\[UnderBracket]Subscript\[UnderBracket]F*q)]
  1. Plot numq:
mode[q_] := FullSimplify[ComplexExpand[Abs[c[3, q*q\[UnderBracket]Subscript\[UnderBracket]c]*v[3, q*q\[UnderBracket]Subscript\[UnderBracket]c, \[Tau]] + 
       s[3, q*q\[UnderBracket]Subscript\[UnderBracket]c]*u[3, q*q\[UnderBracket]Subscript\[UnderBracket]c, \[Tau]]]^2, TargetFunctions -> {Re, Im}]]
numq[q_] = mode[q] /. {q\[UnderBracket]Subscript\[UnderBracket]c -> Exp[1.5]/(v\[UnderBracket]Subscript\[UnderBracket]F*\[Tau])}
Plot[numq[q], {q, 0, 3}, PlotPoints -> 30]

I was wondering what could be the reason for slow plotting. (or simply an error from my code) Is it because I defined the function mode[q] using

:=

? Is it still possible to speed up the plotting process? Remark: I evaluate all the relevant functions at q q_c so that everything is dimensionless.

Thanks in Advance!

$\endgroup$
3
  • $\begingroup$ I suggest displaying each time how long it takes to compute y[nu,q,t] and count how many times you must do that. Your := means that the right hand side of each of those will be calculated again every time. Timing how long FullSimplify needs for each of those will give you more information about where the hour is spending it's time. $\endgroup$
    – Bill
    Commented Oct 15, 2023 at 22:12
  • $\begingroup$ I tried doing as I suggested in the previous comment. I tried various values for q, individual values and things like Table[numq[q],{q,0,3,1/3}] and if I didn't make any mistakes then it appears that DSolve completes in a second or a few seconds for some values of q and other values take so long that I gave up waiting. So I am guessing that the reason your plot is taking more than an hour and displays nothing is that it is grinding away on some "bad" values for q. See if you can reproduce what I'm seeing to confirm this. $\endgroup$
    – Bill
    Commented Oct 16, 2023 at 6:08
  • $\begingroup$ I think the problem comes from ComplexExpand[( Abs[c[3, q Subscript[q, c]] v[3, q Subscript[q, c], t] + s[3, q Subscript[q, c]] u[3, q Subscript[q, c], t]]^2) Somehow, this is an expensive procedure, I changed CompleExpand it to Conjugate[] and the computation is much faster. $\endgroup$ Commented Oct 16, 2023 at 16:07

1 Answer 1

2
$\begingroup$

I found the equations difficult to read with the Symbolized notations, so I rewrote them. I also changed the Greek letters, since they don't paste well on web sites.

If your goal is to produce a plot, you're going to be working with numeric quantities, so running FullSimplify, ComplexExpand, Conjugate, etc. probably isn't going to help much. Also, FullSimplify is usually used in cases where you want to simplify special functions, but unless you encounter those in DSolve, it'll just slow things down as you apply lots of transformation rules which may not result in any simplification. Best to generate one result, examine it, then decide if FullSimplify might help. I removed all of these from your definitions. I also cached the values for w,s, and c since they may be reused during the computation and may not need to be recomputed.

g2[nu_, q_, 
   t_] := (vF*
     Exp[(-(1/2))*(q/qc)^2]*(1 - Cos[(nu*Pi*t)/tau]))/(4*(2*Pi*vF));
w[nu_, q_] := 
  w[nu, q] = Sqrt[(1 + g2[nu, q, tau])^2 - g2[nu, q, tau]^2];
s[nu_, q_] := 
  s[nu, q] = Sqrt[(1/2)*((1 + g2[nu, q, tau])/w[nu, q] - 1)];
c[nu_, q_] := c[nu, q] = Sqrt[1 + s[nu, q]^2];
y[nu_, q_, t_] = 
  a[t] /. First[
    DSolve[{Derivative[2][a][t] + vF^2*q^2*(1 + g2[nu, q, t])*a[t] == 
       0, a[0] == 1, Derivative[1][a][0] == (-I)*vF*q}, a[t], t]];
v[nu_, q_, t_] := (1/2)*y[nu, q, t] - (I*D[y[nu, q, t], t])/(2*vF*q);
u[nu_, q_, t_] := (1/2)*y[nu, q, t] + (I*D[y[nu, q, t], t])/(2*vF*q);

I put a Chop around numq. While running your plot, I saw a number of errors complaining about very small numbers (both real and complex) being generated. There appears to be a singularity near 0.

mode[q_] := 
  Abs[c[3, q*qc]*v[3, q*qc, tau] + s[3, q*qc]*u[3, q*qc, tau]]^2;
numq[q_] = Chop[mode[q] /. {qc -> Exp[1.5]/(vF*tau)}];

Running your plot still appears to take a long time, so I generated a set of points with Table and wrapped ListLinePlot around it. There are still error messages, but you do get an output after 23 seconds. There's a gap near zero, which is where you have a singularity, and the plot is quite flat and near zero for q > 1.5. You might be able to fine-tune your Table generation and improve the resolution (or just use Plot in those regions directly).enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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