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.

I would like to ask your advice how to speed up Contour and Parametric plots in the following example.

Let as start by defining a function

A = {(8 Sqrt[1035741])/15625, 0, 0, 0, 0, Sqrt[6597877/15]/6250, 0, 0, 0, 
0, -(Sqrt[(6338501/10)]/3125), 0, 0, 0, 0, (29 Sqrt[27347/30])/3125, 0, 0, 
0, 0, -(Sqrt[(4358874/5)]/3125), 0, 0, 0, 0, -(Sqrt[(352454081/3)]/31250), 
0, 0, 0};
Clear[s]
s[m_, l_] := If[m < 0, (-1)^m s[-m, l]]
s[m_?NonNegative, l_] := If[Abs[m] <= l, A[[m + 1]], 0]
Clear[ysaf]
L = 28;
ysaf[θ_, ϕ_] = 
ExpToTrig[Sum[s[m, L] SphericalHarmonicY[L, m, θ, ϕ], {m, -5 IntegerPart[L/5], L, 5}]];

ysaf is a pretty complicated trigonometric function of two angles in spherical coordinate system.

Next, we define a few parameters for a stereographic projection. I am plotting a part of the unit sphere as seen from a point on z-axis.

θmx = ArcCos[1/Sqrt[5]];
λ = 3/4;
rmx = Sin[θmx/λ];

Following plot takes a minute or so

gr = ContourPlot[{ysaf[λ ArcSin[Sqrt[x^2 + y^2]] , ArcTan[x, y]] == 0,
Sqrt[x^2 + y^2] == rmx}, {x, -rmx, rmx}, {y, -rmx, rmx}, 
RegionFunction -> Function[{x, y, z}, Sqrt[x^2 + y^2] <= rmx], 
Axes -> False, Frame -> False]

Next, I am constructing interpolating functions for contour lines (like in this post)

lines = Cases[Normal[gr], Line[i__] -> i, Infinity];
fns = Map[Map[ListInterpolation[#, {{0, 1}}] &, Transpose[#]] &, lines];

and define a Schwarz-Christoffel map (see p. 468 of Complex Analysis with MATHEMATICA® by W. T. Shaw, isbn 0521836263):

PolyMap[z_] = Integrate[1/(1 - ξ^5)^(2/5), {ξ, 0, z}];

This final plot is the objective of my question. How can I speed it up?

ParametricPlot[
Table[PolyMap[(Through[fns[[a]][x]]/rmx) /. {List -> Complex}] /. 
Complex[i_, j_] :> {i, j}, {a, 1, First[Dimensions[fns]]}], {x, 0, 1}, Axes -> False, Frame -> False]

plot1

plot2

Any help is appreciated!

share|improve this question
    
Extremely nice plots. Despite some efforts I couldn't make them faster. Will try tomorrow again. –  eldo Aug 29 at 22:05
    
@eldo Thanks a lot for comments! –  yarchik Sep 1 at 8:07

4 Answers 4

up vote 5 down vote accepted

You can get a significant boost by compiling ysaf:

ysafc = Compile[{x, y}, Evaluate @ N @ ysaf[λ ArcSin[Sqrt[x^2 + y^2]], ArcTan[x, y]]];

func[x_?NumericQ, y_?NumericQ] := ysafc[x, y]

gr = ContourPlot[{func[x, y] == 0, Sqrt[x^2 + y^2] == rmx}, {x, -rmx, rmx}, {y, -rmx, rmx}, 
  RegionFunction -> Function[{x, y, z}, Sqrt[x^2 + y^2] <= rmx], 
  Axes -> False, Frame -> False]

For the mapping I would suggest transforming the coordinate data in gr directly. I have also used GenerateConditions -> False in the integral to remove the ConditionalExpression:

PolyMap[z_] = Integrate[1/(1 - ξ^5)^(2/5), {ξ, 0, z}, GenerateConditions -> False];

pm[{x_, y_}] := Through[{Re, Im} @ PolyMap[x + I y]]

Show[Normal[gr] /. Line[data_] :> Line[pm /@ data], PlotRange -> All]

enter image description here

share|improve this answer
    
I learnt a lot from your answer. Even though your code was not the fastest for the ContourPlot I am accepting this answer because it gives an incredible boost was the second part -the Schwarz-Cristoffel transform. This was the slowest part of the code and now it runs in less than 1 sec! –  yarchik Sep 1 at 8:11

Some codes are from the answer of PhilChang

Remove["Global`*"];
A = {(8 Sqrt[1035741])/15625, 0, 0, 0, 0, Sqrt[6597877/15]/6250, 0, 0,
     0, 0, -(Sqrt[(6338501/10)]/3125), 0, 0, 0, 0, (29 Sqrt[27347/30])/3125, 0, 0, 0, 
    0, -(Sqrt[(4358874/5)]/3125), 0, 0, 0, 0, -(Sqrt[(352454081/3)]/31250), 0, 0, 0} // N;
s[m_, l_] := If[Abs[m] > l, 0, If[Negative[m], (-1)^m A[[(-m) + 1]], A[[m + 1]]]]
L = 28;
mlist = Table[m, {m, -5 IntegerPart[L/5], L, 5}];
slist = s[#, L] & /@ mlist;
ysaf[θ_, ϕ_] = Dot[s[#, L] & /@ mlist, SphericalHarmonicY[L, mlist, θ, ϕ]];
θmx = ArcCos[1.0/Sqrt[5.0]];
λ = 3.0/4.0;
rmx = Sin[θmx/λ];

ysafNew = 
 Compile[{{x, _Real}, {y, _Real}}, 
  Evaluate@With[{θ = λ ArcSin[
        Sqrt[x^2 + y^2]], ϕ = ArcTan[x, y]}, 
    Re@Dot[slist, SphericalHarmonicY[L, mlist, θ, ϕ]]]];

ContourPlot[{ysafNew[x, y] == 0.0, Sqrt[x^2 + y^2] == rmx}, {x, -rmx, 
    rmx}, {y, -rmx, rmx}, 
   RegionFunction -> Function[{x, y, z}, Sqrt[x^2 + y^2] <= rmx], 
   Axes -> False, Frame -> False] // Quiet // AbsoluteTiming

Mathematica graphics

share|improve this answer
    
hi,do you know how to optimize fns and ParametricPlot? i can not understand what it means –  PhilChang Aug 30 at 9:45
    
@ Chenminqi When I tried your solution on my machine I got incredible 1.4 sec! Combining your answer with the solution of @Simon Woods for the Schwarz-Cristoffel mapping I am getting 1.4 + 0.8 sec in total! –  yarchik Sep 1 at 8:20

One can gain a considerable speed gain by omitting RegionFunction and setting PlotPoints to 10. The following only needs 8 seconds (it uses some ideas of @PhilChang):

A = {
   (8 Sqrt[1035741])/15625, 0, 0, 0, 0,
   Sqrt[6597877/15]/6250, 0, 0, 0, 0,
   -(Sqrt[(6338501/10)]/3125), 0, 0, 0, 0,
   (29 Sqrt[27347/30])/3125, 0, 0, 0, 0,
   -(Sqrt[(4358874/5)]/3125), 0, 0, 0, 0,
   -(Sqrt[(352454081/3)]/31250), 0, 0, 0
   };

s[m_, l_] := If[Abs[m] > l, 0, If[Negative[m], (-1)^m A[[(-m) + 1]], A[[m + 1]]]]

L = 16;
mx = ArcCos[1/Sqrt[5]];
lambda = 3/4;
rmx = Sin[mx/lambda];
mlist = Table[m, {m, -5 IntegerPart[L/5], L, 5}];
slist = s[#, L] & /@ mlist;
hlist = SphericalHarmonicY[L, #, a, b] & /@ mlist;

ysaf[a_, b_] = Dot[slist, hlist];

ContourPlot[
  ysaf[lambda ArcSin[Sqrt[x^2 + y^2]], ArcTan[x, y]],
  {x, -rmx, rmx}, {y, -rmx, rmx},
  ColorFunction -> "DarkRainbow",
  Contours -> 6,
  PlotPoints -> 10,
  PlotLegends -> Automatic]

enter image description here

Plot3D[
  ysaf[lambda ArcSin[Sqrt[x^2 + y^2]], ArcTan[x, y]],
  {x, -rmx, rmx}, {y, -rmx, rmx},
  ColorFunction -> "DarkRainbow",
  PlotPoints -> 50,
  ImageSize -> 500,
  PlotLegends -> Automatic]

enter image description here

share|improve this answer

The orginal program takes 78s on my computer. First, Make vector A a real vector will make the program a little faster.

 A = {(8 Sqrt[1035741])/15625, 0, 0, 0, 0, Sqrt[6597877/15]/6250, 0, 0,
     0, 0, -(Sqrt[(6338501/10)]/3125), 0, 0, 0, 
    0, (29 Sqrt[27347/30])/3125, 0, 0, 0, 
    0, -(Sqrt[(4358874/5)]/3125), 0, 0, 0, 
    0, -(Sqrt[(352454081/3)]/31250), 0, 0, 0} // N;

Second, avoiding recursive in function s also help.

s[m_, l_] := If[Abs[m] > l, 0,
    If[Negative[m], (-1)^m A[[(-m) + 1]], A[[m + 1]] ]

Third, Sum is a little bit slow. Use Dot instead. And ExpToTrig is not necessary.

mlist = Table[m, {m, -5 IntegerPart[L/5], L, 5}];
ysaf[θ_, ϕ_] = 

   Dot[s[#, L] & /@ mlist, 
    SphericalHarmonicY[L, #, θ, ϕ] & /@ mlist];

Fourth,Evaluate

gr = ContourPlot[
  Evaluate@{  
    ysaf[λ ArcSin[Sqrt[x^2 + y^2]], ArcTan[x, y]] == 0.0, 
    Sqrt[x^2 + y^2] == rmx}, {x, -rmx, rmx}, {y, -rmx, rmx}, 
  RegionFunction -> Function[{x, y, z}, Sqrt[x^2 + y^2] <= rmx], 
  Axes -> False, Frame -> False]

In this way, the contour plot gr cost 20s compared with original 78s on my computer.

Besides, if you define your function ysaf as

ysafNew[x_, y_] := 
 With[{  θ = λ ArcSin[Sqrt[x^2 + y^2]], ϕ = 
    ArcTan[x, y]},
                Re@Dot[slist, SphericalHarmonicY[L, mlist, θ, ϕ]]

the time will be 9s.

I hope someone could improve more.

In ParametricPlot, Use MaxRecursion and Plotpoints

ParametricPlot[
Table[PolyMap[(Through[fns[[a]][x]]/rmx) /. {List -> Complex}] /. 
Complex[i_, j_] :> {i, j}, {a, 1, First[Dimensions[fns]]}],
{x, 0, 1}, Axes -> False, Frame -> False, 
MaxRecursion -> 2, 
PlotPoints -> 30]

============================================================================================== Here is all the code:

    Remove["Global`*"];
A = {(8 Sqrt[1035741])/15625, 0, 0, 0, 0, Sqrt[6597877/15]/6250, 0, 0,
     0, 0, -(Sqrt[(6338501/10)]/3125), 0, 0, 0, 
    0, (29 Sqrt[27347/30])/3125, 0, 0, 0, 
    0, -(Sqrt[(4358874/5)]/3125), 0, 0, 0, 
    0, -(Sqrt[(352454081/3)]/31250), 0, 0, 0} // N;
Clear[s]
(*s[m_,l_]:=If[m<0,(-1)^m s[-m,l]]
s[m_?NonNegative,l_]:=If[Abs[m]≤l,A[[m+1]],0]*)
s[m_, l_] := If[Abs[m] > l, 0,
    If[Negative[m], (-1)^m A[[(-m) + 1]], A[[m + 1]] ]
  ]
Clear[ysaf]
L = 28;
(*ysaf[θ_,ϕ_]=ExpToTrig[Sum[s[m,L] \
SphericalHarmonicY[L,m,θ,ϕ],{m,-5 \
IntegerPart[L/5],L,5}]];*)
mlist = Table[m, {m, -5 IntegerPart[L/5], L, 5}];
slist = s[#, L] & /@ mlist; 
ysaf[θ_, ϕ_] = 
  Dot[s[#, L] & /@ mlist, 
   SphericalHarmonicY[L, mlist, θ, ϕ]] ;
θmx = ArcCos[1.0/Sqrt[5.0]];
λ = 3.0/4.0;
rmx = Sin[θmx/λ];
ysafNew[x_, y_] := 
  With[{  θ = λ ArcSin[Sqrt[x^2 + y^2]], ϕ = 
     ArcTan[x, y]},
                Re@Dot[slist, SphericalHarmonicY[L, mlist, θ, ϕ]]
        ];
{time, gr} = 
 ContourPlot[
   Evaluate@{ysafNew[x, y] == 0.0, Sqrt[x^2 + y^2] == rmx}, {x, -rmx, 
    rmx}, {y, -rmx, rmx}, 
   RegionFunction -> Function[{x, y, z}, Sqrt[x^2 + y^2] <= rmx], 
   Axes -> False, Frame -> False] // AbsoluteTiming

enter image description here enter image description here

share|improve this answer
    
Aside from the fact that your answer contains several typos I cannot confirm at all your timings. Any proof? –  eldo Aug 30 at 8:20
    
@eldo, I pasted all my source code in my answer, if you can not run the code, please let me know. I did not understand how to get fns and the function inside the ParametricPlot very much, so I can not optimize the performance. –  PhilChang Aug 30 at 8:47

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.