1
$\begingroup$

I am trying to plot the roots of a function f(x,y) using ContourPlot as follows

Sk[j_] = Sqrt[1 - x^2 + j 2 y + y^2];
ρ[j_] = (1 - x^2 + y^2 - Sk[j]^2)/(2 y);
A[s_, j_] = ( (s I x + Sk[j]) (y + ρ[j]))/(s - s y^2);
ζ[s_, j_] = ( (s I x + Sk[j]) (1 + y ρ[j]))/(s - s y^2);

Gs1[y_] = {{A[-1, 1], A[1, 1], A[-1, -1], A[1, -1]}, {1, 1, 1, 
    1}, {ζ[-1, 1], ζ[1, 1], ζ[-1, -1], ζ[
     1, -1]}, {ρ[1], ρ[1], ρ[-1], ρ[-1]}};
Ms1[z_, y_] = {{E^(I Sk[1] z), 0, 0, 0}, {0, E^(-I Sk[1] z), 0, 
    0}, {0, 0, E^(I Sk[-1] z), 0}, {0, 0, 0, E^(-I Sk[-1] z)}};
Mxy = Inverse[Ms1[0, y]].Inverse[Gs1[y]].Gs1[y - 3].Ms1[0, 
    y - 3].Inverse[Ms1[8, y - 3]].Inverse[Gs1[y - 3]].Gs1[y].Ms1[8, 
    y];
ContourPlot[{N@Mxy[[4, 3]] == 0}, {x, -2, 2}, {y, 0.001, 3}, 
  Exclusions -> y == 1, ExclusionsStyle -> None, PlotPoints -> 30, 
  ContourStyle -> {Directive[Blue]}, 
  RegionFunction -> 
   Function[{x, y, f},  Abs[y - 1] - Abs[x] < -0.08 ]] // Timing

I had two problems:
1- In M10 (win10) it takes about 73.46 sec to get the results (see Fig. below) while in M12 (12.0.0.0 Linux x86 64-bit) I wait for more than 15 min but no results. So, why is this discrepancy?

enter image description here

2- To speed up the computation and enhance the results (smoothen the curves) I tried to define the function Mxy[[4, 3]] with numeric variables as follows

f[x_?NumericQ, y_?NumericQ] = Mxy[[4, 3]];   

ContourPlot[{f[x, y] == 0}, {x, -2, 2}, {y, 0.001, 3}, 
  Exclusions -> y == 1, ExclusionsStyle -> None, PlotPoints -> 30, 
  ContourStyle -> {Directive[Blue]}, 
  RegionFunction -> 
   Function[{x, y, f},  Abs[y - 1] - Abs[x] < -0.08 ]] // Timing    

but this takes a long time with the same output, see Fig. below! So, How can I speed up the computation and smooth the curves and why defining numeric variables did not help?

enter image description here

$\endgroup$

1 Answer 1

5
$\begingroup$
$Version

(* "12.1.1 for Mac OS X x86 (64-bit) (June 19, 2020)" *)

Clear["Global`*"]

Sk[j_] = Sqrt[1 - x^2 + j 2 y + y^2];
ρ[j_] = (1 - x^2 + y^2 - Sk[j]^2)/(2 y);
A[s_, j_] = ((s I x + Sk[j]) (y + ρ[j]))/(s - s y^2);
ζ[s_, j_] = ((s I x + Sk[j]) (1 + y ρ[j]))/(s - s y^2);

Gs1[y_] = {{A[-1, 1], A[1, 1], A[-1, -1], A[1, -1]}, {1, 1, 1, 
    1}, {ζ[-1, 1], ζ[1, 1], ζ[-1, -1], ζ[
     1, -1]}, {ρ[1], ρ[1], ρ[-1], ρ[-1]}};
Ms1[z_, y_] = {{E^(I Sk[1] z), 0, 0, 0}, {0, E^(-I Sk[1] z), 0, 0}, {0, 0, 
    E^(I Sk[-1] z), 0}, {0, 0, 0, E^(-I Sk[-1] z)}};
Mxy = Inverse[Ms1[0, y]].Inverse[Gs1[y]].Gs1[y - 3].Ms1[0, y - 3].Inverse[
    Ms1[8, y - 3]].Inverse[Gs1[y - 3]].Gs1[y].Ms1[8, y];

Simplify the expression prior to using it in ContourPlot

expr = Mxy[[4, 3]] // Simplify

(* -((3 E^(-8 I (Sqrt[-x^2 + (-4 + y)^2] - Sqrt[-x^2 + (-1 + y)^2])) (-1 + E^(
       16 I Sqrt[-x^2 + (-4 + y)^2])) x (x + 
       I Sqrt[-x^2 + (-1 + y)^2]))/(2 Sqrt[-x^2 + (-4 + y)^2]
      Sqrt[-x^2 + (-1 + y)^2] (-1 + y))) *)

ContourPlot[expr == 0, {x, -2, 2}, {y, 0.001, 3}, Exclusions -> y == 1, 
  ExclusionsStyle -> None, PlotPoints -> 30, 
  ContourStyle -> {Directive[Blue]}, 
  RegionFunction -> 
   Function[{x, y, f}, Abs[y - 1] - Abs[x] < -0.08]] // Timing

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.