Skip to main content
Mathematica

Return to Answer

added 57 characters in body
Source Link
Michael E2
Michael E2
  • 244.7k
  • 18
  • 351
  • 774

Here's a polynomial interpolation method, which can be be found in Chapter 5 of Boyd (2014).

nn = 64;
z0 = w1 + w2;
rr = 1.1 w1;
ff = N[WeierstrassP[z0 + rr #, inv] - L, Precision[#]] &;
wprec = MachinePrecision;
tj = 2 Pi*Range[0, nn - 1]/nn;
wj = N[Exp[I tj], wprec];
fj = ff /@ wj; (* f[zj] *)
aa = InverseFourier[fj]/Sqrt[nn];

(* Rough check of accuracy of interpolation *)
"condition"@# -> Log10@Ratios[#] &@N@MinMax@Abs@fj&@ N@ MinMax@ Abs@ fj
ip = FromDigits[Reverse@aa, (z - z0)/rr];
Max@Table[Max@ Table[(ff[z]WeierstrassP[z, inv] - L - ip)/ff[z](0WeierstrassP[z, inv] - L) /. 
    z -> z0 + r Exp[I t] // N // Abs,
  {r, (1/8 - 0.01) rr, rr, rr/8}, {t, 0., 2 Pi, 0.2/r}]
(*
  "condition"[{0.0694093, 0.75655}] -> {1.03742}
  31.98344*10^7511*10^-1112
*)

z2 = Eigenvalues@companionMatrix[aa];Eigenvalues@ companionMatrix[aa];
roots = z0 + rr*Select[z2, Abs[#] < 0.999 &]
(*  {0.776416 + 2.35619 I, 2.36518 + 2.35619 I}  *)

Graphics[{
  EdgeForm@Gray, LightBlue, Disk[ReIm@z0, rr], Gray, Point[ReIm@z0],
  First@plt,
  Black, PointSize[Medium], Point@ReIm[roots]},
 Frame -> True]

Mathematica graphics

Auxiliary code:

For the companion matrix, you can use

companionMatrix[coeffs_] := 
  Join[SparseArray[
    Band[{1, 2}] -> 1, {Length@coeffs - 2, Length@coeffs - 1}],
   {-coeffs[[;; -2]]/coeffs[[-1]]}
   ];

Or

companionMatrix = NRoots`CompanionMatrix

Here's a polynomial interpolation method, which can be be found in Chapter 5 of Boyd (2014).

nn = 64;
z0 = w1 + w2;
rr = 1.1 w1;
ff = N[WeierstrassP[z0 + rr #, inv] - L, Precision[#]] &;
wprec = MachinePrecision;
tj = 2 Pi*Range[0, nn - 1]/nn;
wj = N[Exp[I tj], wprec];
fj = ff /@ wj; (* f[zj] *)
aa = InverseFourier[fj]/Sqrt[nn];

(* Rough check of accuracy of interpolation *)
"condition"@# -> Log10@Ratios[#] &@N@MinMax@Abs@fj
ip = FromDigits[Reverse@aa, (z - z0)/rr];
Max@Table[(ff[z] - ip)/ff[z] /. z -> z0 + r Exp[I t] // N // Abs,
  {r, (1/8 - 0.01) rr, rr, rr/8}, {t, 0., 2 Pi, 0.2/r}]
(*
  "condition"[{0.0694093, 0.75655}] -> {1.03742}
  3.98344*10^-11
*)

z2 = Eigenvalues@companionMatrix[aa];
roots = z0 + rr*Select[z2, Abs[#] < 0.999 &]
(*  {0.776416 + 2.35619 I, 2.36518 + 2.35619 I}  *)

Graphics[{
  EdgeForm@Gray, LightBlue, Disk[ReIm@z0, rr], Gray, Point[ReIm@z0],
  First@plt,
  Black, PointSize[Medium], Point@ReIm[roots]},
 Frame -> True]

Mathematica graphics

Auxiliary code:

For the companion matrix, you can use

companionMatrix[coeffs_] := 
  Join[SparseArray[
    Band[{1, 2}] -> 1, {Length@coeffs - 2, Length@coeffs - 1}],
   {-coeffs[[;; -2]]/coeffs[[-1]]}
   ];

Or

companionMatrix = NRoots`CompanionMatrix

Here's a polynomial interpolation method, which can be be found in Chapter 5 of Boyd (2014).

nn = 64;
z0 = w1 + w2;
rr = 1.1 w1;
ff = N[WeierstrassP[z0 + rr #, inv] - L, Precision[#]] &;
wprec = MachinePrecision;
tj = 2 Pi*Range[0, nn - 1]/nn;
wj = N[Exp[I tj], wprec];
fj = ff /@ wj; (* f[zj] *)
aa = InverseFourier[fj]/Sqrt[nn];

(* Rough check of accuracy of interpolation *)
"condition"@# -> Log10@Ratios[#] &@ N@ MinMax@ Abs@ fj
ip = FromDigits[Reverse@aa, (z - z0)/rr];
Max@ Table[(WeierstrassP[z, inv] - L - ip)/(0WeierstrassP[z, inv] - L) /. 
    z -> z0 + r Exp[I t] // N // Abs,
  {r, (1/8 - 0.01) rr, rr, rr/8}, {t, 0., 2 Pi, 0.2/r}]
(*
  "condition"[{0.0694093, 0.75655}] -> {1.03742}
  1.7511*10^-12
*)

z2 = Eigenvalues@ companionMatrix[aa];
roots = z0 + rr*Select[z2, Abs[#] < 0.999 &]
(*  {0.776416 + 2.35619 I, 2.36518 + 2.35619 I}  *)

Graphics[{
  EdgeForm@Gray, LightBlue, Disk[ReIm@z0, rr], Gray, Point[ReIm@z0],
  First@plt,
  Black, PointSize[Medium], Point@ReIm[roots]},
 Frame -> True]

Mathematica graphics

Auxiliary code:

For the companion matrix, you can use

companionMatrix[coeffs_] := 
  Join[SparseArray[
    Band[{1, 2}] -> 1, {Length@coeffs - 2, Length@coeffs - 1}],
   {-coeffs[[;; -2]]/coeffs[[-1]]}
   ];

Or

companionMatrix = NRoots`CompanionMatrix
Source Link
Michael E2
Michael E2
  • 244.7k
  • 18
  • 351
  • 774

Here's a polynomial interpolation method, which can be be found in Chapter 5 of Boyd (2014).

nn = 64;
z0 = w1 + w2;
rr = 1.1 w1;
ff = N[WeierstrassP[z0 + rr #, inv] - L, Precision[#]] &;
wprec = MachinePrecision;
tj = 2 Pi*Range[0, nn - 1]/nn;
wj = N[Exp[I tj], wprec];
fj = ff /@ wj; (* f[zj] *)
aa = InverseFourier[fj]/Sqrt[nn];

(* Rough check of accuracy of interpolation *)
"condition"@# -> Log10@Ratios[#] &@N@MinMax@Abs@fj
ip = FromDigits[Reverse@aa, (z - z0)/rr];
Max@Table[(ff[z] - ip)/ff[z] /. z -> z0 + r Exp[I t] // N // Abs,
  {r, (1/8 - 0.01) rr, rr, rr/8}, {t, 0., 2 Pi, 0.2/r}]
(*
  "condition"[{0.0694093, 0.75655}] -> {1.03742}
  3.98344*10^-11
*)

z2 = Eigenvalues@companionMatrix[aa];
roots = z0 + rr*Select[z2, Abs[#] < 0.999 &]
(*  {0.776416 + 2.35619 I, 2.36518 + 2.35619 I}  *)

Graphics[{
  EdgeForm@Gray, LightBlue, Disk[ReIm@z0, rr], Gray, Point[ReIm@z0],
  First@plt,
  Black, PointSize[Medium], Point@ReIm[roots]},
 Frame -> True]

Mathematica graphics

Auxiliary code:

For the companion matrix, you can use

companionMatrix[coeffs_] := 
  Join[SparseArray[
    Band[{1, 2}] -> 1, {Length@coeffs - 2, Length@coeffs - 1}],
   {-coeffs[[;; -2]]/coeffs[[-1]]}
   ];

Or

companionMatrix = NRoots`CompanionMatrix