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.

When I compute the phase error of a spatial series data using Fourier analysis in Mathematica there's a discontinuity @ parameter c1 = 1.35. However @ c1 = 0.5 produces the correct result.

Code:

    Clear[G, σ, ϕ];
G = -σ/2*((1 - Cos[ϕ])^2 + I*(3 - Cos[ϕ])*Sin[ϕ]);
Ztri = (1 + G + 1/2*G^2 + 1/6*G^3 + 1/24*G^4);
g[σ_Real, ϕ_Real] = -ArcTan[Re[Ztri], Im[Ztri]]/(σ*ϕ);
linecolors=Blue;
framecolors=Black;
c1 = 1.35
gp1 = Plot[g[σ, ϕ] /. {σ -> c1}, {ϕ, 0, Pi}, 
   PlotRange -> {-2, 2.}, PlotStyle -> {linecolors, Thickness[0.006]},
    PlotLegends -> Placed[{"CFL 1.35"}, {0.2, 0.4}], 
   AspectRatio -> Automatic];
c1 = 0.5;
gp2 = Plot[g[σ, ϕ] /. {σ -> c1}, {ϕ, 0, Pi}, 
   PlotRange -> {-2, 2.}, PlotStyle -> {linecolors,Dotted, Thickness[0.006]},
    PlotLegends -> Placed[{"CFL 0.5"}, {0.2, 0.4}], 
   AspectRatio -> Automatic];
BB = Show[gp1, gp2, Axes -> False, Frame -> True, 
  FrameStyle -> Directive[Thick, framecolors, 15], 
  FrameLabel -> {{"Phase error", ""}, {ω, 
     "Numerical dispersion"}}]

Plot after running code: enter image description here

The correct plot is like this Correct plot for c1=CFL=1.35

share|improve this question
1  
Related: (5782), (11714) –  Mr.Wizard Jul 31 at 0:28

2 Answers 2

Please tell me if this produces what you expect:

link[a : {{_, _} ..}, b : {{_, _} ..}] := 
  Join[a, (b\[Transpose] - First[b] + Last[a])\[Transpose]]

link[x__] := Fold[link, {x}]

gp1fix = gp1 /. Line[x_] :> Line[link @@ Split[x, EuclideanDistance[##] < 1 &]];

Show[gp1fix, gp2, Axes -> False, Frame -> True, FrameStyle -> Directive[Thick, 15], 
 FrameLabel -> {{"Phase error", ""}, {\[Omega], "Numerical dispersion"}}]

enter image description here

share|improve this answer
    
Thanks a lot. I guess this is what I want. –  Flymath Jul 31 at 2:57
    
Now it seem that another question appear: the linked solution is not so smooth compared with correct plot( cfl=1.35 should be corresponding to clf=1) I give. In addition, let gp1=..., PlotRange -> {-0, 2.}..., the resulted plot is still discontious. –  Flymath Jul 31 at 3:02
3  
@Flymath: Ah, the problem is that one needs to "unwrap" the result of ArcTan before dividing by $\sigma\phi$. In this case it is apparently enough to just define g[σ_Real, ϕ_Real] = (If[# < 0, # + 2 π, #] &)@-ArcTan[Re[Ztri], Im[Ztri]]/(σ*ϕ). –  Rahul Narain Jul 31 at 6:52
    
@Rahul Good observation. I wasn't digging that deep. –  Mr.Wizard Jul 31 at 6:55
    
Thank Rahul Narain a lot for finding the problem. Also thank for Mr. Wizard again. –  Flymath Jul 31 at 13:10

This comes from the jump disontinuity of two argument ArcTan. My solution can be automated, but I'll leave that to you for the time being.

My strategy is to find where the discontinuity is, then lift the right part of the graph up by a constant.

The jump of two argument ArcTan occurs when $x < 0$ and $y = 0$:

Plot3D[ArcTan[x, y], {x, -π, π}, {y, -π, π}, AxesLabel -> {x, y}, Exclusions -> {{y == 0, x <= 0}}]

enter image description here

You have

g[σ, φ] == -(1/(σ φ))ArcTan[384 + Re[σ ((-1 + Cos[φ])^2 - 
    I (-3 + Cos[φ]) Sin[φ]) (-192 + 
    48 σ ((-1 + Cos[φ])^2 - 
      I (-3 + Cos[φ]) Sin[φ]) - 
    8 σ^2 ((-1 + Cos[φ])^2 - 
      I (-3 + Cos[φ]) Sin[φ])^2 + σ^3 ((-1 + 
        Cos[φ])^2 - I (-3 + Cos[φ]) Sin[φ])^3)], 
   Im[σ ((-1 + Cos[φ])^2 - 
    I (-3 + Cos[φ]) Sin[φ]) (-192 + 
    48 σ ((-1 + Cos[φ])^2 - 
     I (-3 + Cos[φ]) Sin[φ]) - 
    8 σ^2 ((-1 + Cos[φ])^2 - 
     I (-3 + Cos[φ]) Sin[φ])^2 + σ^3 ((-1 + 
       Cos[φ])^2 - I (-3 + Cos[φ]) Sin[φ])^3)]]

so we need to solve Im[__] == 0, when σ == 1.35.

Reduce[Im[σ ((-1 + Cos[φ])^2 - 
    I (-3 + Cos[φ]) Sin[φ]) (-192 + 
    48 σ ((-1 + Cos[φ])^2 - 
     I (-3 + Cos[φ]) Sin[φ]) - 
    8 σ^2 ((-1 + Cos[φ])^2 - 
     I (-3 + Cos[φ]) Sin[φ])^2 + σ^3 ((-1 + 
       Cos[φ])^2 - I (-3 + Cos[φ]) Sin[φ])^3)] == 0 && 
   2 < φ < 3 /. σ -> Rationalize[1.35], φ] // RootReduce

(* φ == 2 ArcTan[Root[-2000 - 9570 #1^2 - 10020 #1^4 + 6629 #1^6 - 
            1974 #1^8 - 21504 #1^10 + 2671 #1^12 &, 2]] *)

This is the location of the discontinuity. From here you can also find the difference between the left and right limits to get the offset too:

loc = 2 ArcTan[Root[-2000 - 9570 #1^2 - 10020 #1^4 + 6629 #1^6 - 
          1974 #1^8 - 21504 #1^10 + 2671 #1^12 &, 2]];

offset = -2 g[Rationalize[1.35], loc];

Now add this only when x > loc. Here are your first four lines now:

Clear[G, Ztri, σ, φ];
G[σ_, φ_] := -σ/2*((1 - Cos[φ])^2 + I*(3 - Cos[φ])*Sin[φ]);
Ztri[σ_, φ_] := (1 + G[σ, φ] + 1/2*G[σ, φ]^2 + 1/6*G[σ, φ]^3 + 1/24*G[σ, φ]^4);
g[σ_, φ_] := -ArcTan[Re[Ztri[σ, φ]], Im[Ztri[σ, φ]]]/(σ*φ) + offset Boole[σ == 1.35 && φ > loc];

Again, this is not a general solution and only works for σ == 1.35. Here is the plot now: enter image description here

share|improve this answer
    
Thank Chip Hurst. You give me the detailed explaination. –  Flymath Jul 31 at 3:03
    
In fact, now the line of cfl=1.35 should be that of cfl=1 in the given correct plot. Why the line not discontinuous? And the max-values also different? I think that it is sure that the expressions of G and Ztri are correct. –  Flymath Jul 31 at 3:11

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.