2
$\begingroup$

I have the following equation

$$\pi_{xx} = \exp(-2 k_x^2 \sin^2(\frac{t}{2})) k_x^2 \theta(t) \times \bigg[ -\frac{\sin(2t)}{2} - \cos(t) \sin(t) (1 - k_x^2 \sin^2(\frac{t}{2})) + 2 \cos(t) \sin^2(\frac{t}{2}) \sin^2(t) + 2 \sin^2(t) - 4 k_x^2 \sin^2(t) \sin^2(\frac{t}{2}) - 2 \sin(2t) \cos(t) \sin^2(t) - \sin(2t) \sin^2(t) + 2 k_x^2 \sin(2t) \sin^2(t) \sin^2(\frac{t}{2})\bigg]$$

enter image description here

I want to find maximums of this function (points where the time derivative is zero). The goal here is to have a comparison with my code result, so it should be accurate. Using Table or FindPeaks will not be accurate enough. The proper way is to use D to compute the time derivative, and FindRoots to find the maximum values.

k here is the wave number. It could be {10, 20, 30, 40, ...}.

What I tried so far is as follows.

xxrot = E^(-2 k^2 Sin[t/2]^2)k^2 Cos[0.4292036732051034 + t]^2 (-Sin[t] (1 - k^2 Sin[t]^2) - Sin[2 t]) + E^(-2 k^2 Sin[t/2]^2) k^2 (-4 Sin[t/2]^2 Sin[t] - (1 - 4 k^2 Sin[t/2]^4) Sin[t]) Sin[0.4292036732051034 + t]^2 - E^(-2 k^2 Sin[t/2]^2)k^2 (2 Cos[t] Sin[t/2]^2 + Sin[t]^2 - 2 k^2 Sin[t/2]^2 Sin[t]^2) Sin[2 (0.4292036732051034 + t)]; 
    
xxrotplt = Plot[xxrot /. k -> -i2k@20 /. t -> t - 2 // Evaluate, {t, 0, 2000 0.01}, PlotRange -> All, PlotStyle -> Red];
    
i2k[i_] = 1/2 (i - 1) ;
    
pl = xxrot /. k -> -i2k[80] /. t -> t - 200*0.01;
tab = Table[pl, {t, 8.1, 8.3, 0.01}];
pp = Rest@FindPeaks[tab];
plta = ListLinePlot[tab, 
  Epilog -> {Red, PointSize[Large], Point /@ pp}, 
  GridLines -> Automatic, ImageSize -> Large, PlotRange -> All]
pts = Cases[plta, 
   Point[{x_, y_ /; y > 20}] :> {x/100, y}, \[Infinity]] // Sort

Here the wave number is 80, and the result of the above code is enter image description here

You can see how bad it is for high wave numbers and it is not possible to compare this result with the code result.

However, when I used D and FindRoots it got more complicated (it is not helpful to put my try here).

I appreciate it if you could help me.

$\endgroup$
1
  • $\begingroup$ Unclear: Are you looking for "points where the time derivative is zero" or global maximum? $\endgroup$ Oct 16, 2023 at 20:09

2 Answers 2

1
$\begingroup$

Here is your function written out as a fuction of t and k:

xxrot[t_, k_] := 
  E^(-2 k^2 Sin[t/2]^2) k^2 Cos[
      0.4292036732051034 + t]^2 (-Sin[t] (1 - k^2 Sin[t]^2) - 
      Sin[2 t]) + 
   E^(-2 k^2 Sin[t/2]^2) k^2 (-4 Sin[t/2]^2 Sin[
        t] - (1 - 4 k^2 Sin[t/2]^4) Sin[t]) Sin[
      0.4292036732051034 + t]^2 - 
   E^(-2 k^2 Sin[t/2]^2) k^2 (2 Cos[t] Sin[t/2]^2 + Sin[t]^2 - 
      2 k^2 Sin[t/2]^2 Sin[t]^2) Sin[2 (0.4292036732051034 + t)];

If you plot it, you'll see if matches your plot above:

Plot[xxrot[t, 10], {t, 0, 20}, PlotRange -> All]

Take the derivative:

dxx[t_] := D[xxrot[t, 10], t]

Now use FindRoot to find values near the peaks:

FindRoot[dxx[t], {t, #}] & /@ {Range[5, 6.5, 0.01]}

This gives out a collection of answers, some of which have failed. But among them you will find

t-> 6.207973968942875`

for which

xxrot[6.207973968942875`, 10]

(* 12.4548 *)

is the largest peak in the region. So I would just search through the returned values to find the largest.

$\endgroup$
3
  • $\begingroup$ Thanks bill. However, how can I know which t is the right one? your code produces many answers. How can I quickly find the correct one among them? $\endgroup$
    – Lohrasb
    Oct 18, 2023 at 11:21
  • $\begingroup$ And actually, what is the reason of having many answers? $\endgroup$
    – Lohrasb
    Oct 18, 2023 at 11:22
  • $\begingroup$ If you evaluate xxrot at all the candidate answers, one of them will be the largest... this is the one you are looking for. $\endgroup$
    – bill s
    Oct 18, 2023 at 16:52
3
$\begingroup$
$Version

(* "13.3.1 for Mac OS X ARM (64-bit) (July 24, 2023)" *)

Clear["Global`*"]

xxrot[t_, k_] := 
  E^(-2 k^2 Sin[t/2]^2) k^2 Cos[
       0.4292036732051034 + t]^2 (-Sin[t] (1 - k^2 Sin[t]^2) - 
       Sin[2 t]) + 
    E^(-2 k^2 Sin[t/2]^2) k^2 (-4 Sin[t/2]^2 Sin[
         t] - (1 - 4 k^2 Sin[t/2]^4) Sin[t]) Sin[
       0.4292036732051034 + t]^2 - 
    E^(-2 k^2 Sin[t/2]^2) k^2 (2 Cos[t] Sin[t/2]^2 + Sin[t]^2 - 
       2 k^2 Sin[t/2]^2 Sin[t]^2) Sin[2 (0.4292036732051034 + t)] // 
   Rationalize[#, 0] &;

sol = Table[{k, 
   Maximize[{xxrot[t, k], 0 <= t <= 20}, t] // N[#, 50] & // N}, 
   {k, 10, 80, 10}]

(* {{10, {12.4548, {t -> 18.7743}}}, {20, {24.848, {t -> 
     18.8119}}}, {30, {37.2395, {t -> 18.8244}}}, {40, {49.6305, {t ->
      18.8307}}}, {50, {62.0213, {t -> 18.8345}}}, {60, {74.412, {t ->
      18.837}}}, {70, {86.8026, {t -> 18.8388}}}, {80, {99.1932, {t ->
      18.8401}}}} *)

Plot[xxrot[t, 80], {t, 187/10, 19},
 WorkingPrecision -> 50,
 PlotRange -> All,
 MaxRecursion -> 5,
 Epilog -> {Red, AbsolutePointSize[6],
   Point[{t, xxrot[t, 80]}] /. sol[[-1, -1, -1]]}]

enter image description here

EDT: For each of the wave numbers,

Plot[Evaluate[
  xxrot[t, #] & /@ Range[10, 80, 10]],
 {t, 1874/100, 1886/100},
 Frame -> True,
 PlotRange -> {-30, All},
 Epilog -> {Red, AbsolutePointSize[4],
   Point[({t, xxrot[t, #[[1]]]} /. #[[2, 2, -1]]) & /@
     sol]},
 PlotLegends -> Placed[
   LineLegend[
    Range[10, 80, 10],
    LegendLabel -> Style["k", 14],
    LegendLayout -> {"Column", 2}],
   {.3, .7}]]

enter image description here

EDIT 2: To find all of the maxima for each wave number

The first derivative w.r.t. t is

der[t_, k_] = D[xxrot[t, k], t] //
   FullSimplify[#, k > 0 && t >= 0] &;

The second derivative w.r.t. t is

der2[t_, k_] = D[der[t, k], t] //
   FullSimplify[#, k > 0 && t >= 0] &;

The maxima occur at

Column[(sol2 = Table[{k, Solve[
      {der[t, k] == 0, der2[t, k] < 0, 0 <= t <= 20},
      t, Reals]},
    {k, 10, 80, 10}]) /. x_Root :> N[N[x, 50]]]

enter image description here

Plotting,

Manipulate[
 Plot[xxrot[t, k], {t, 0, 20},
  WorkingPrecision -> 50,
  PlotPoints -> 100,
  MaxRecursion -> 7,
  PlotRange -> {-5, All},
  Epilog -> {Red, AbsolutePointSize[4],
    Point[{t, xxrot[t, k]} /. sol2[[k/10, 2]]]}],
 {{k, 10}, Range[10, 80, 10], 
    ControlType -> RadioButtonBar},
 SynchronousUpdating -> False,
 TrackedSymbols :> {k}]

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.