# Fourier Analysis: How to get rid of a discontinuity

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: The correct plot is like this 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}}] 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: • Thank Chip Hurst. You give me the detailed explaination. – Flymath Jul 31 '14 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 '14 at 3:11

Please tell me if this produces what you expect:

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

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"}}] • Thanks a lot. I guess this is what I want. – Flymath Jul 31 '14 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 '14 at 3:02
• @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 Jul 31 '14 at 6:52
• @Rahul Good observation. I wasn't digging that deep. – Mr.Wizard Jul 31 '14 at 6:55
• Thank Rahul Narain a lot for finding the problem. Also thank for Mr. Wizard again. – Flymath Jul 31 '14 at 13:10