# 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

• Related: (5782), (11714) Jul 31, 2014 at 0:28

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. Jul 31, 2014 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. Jul 31, 2014 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. Jul 31, 2014 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. Jul 31, 2014 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]]/(σ*ϕ).
– user484
Jul 31, 2014 at 6:52
• @Rahul Good observation. I wasn't digging that deep. Jul 31, 2014 at 6:55
• Thank Rahul Narain a lot for finding the problem. Also thank for Mr. Wizard again. Jul 31, 2014 at 13:10