# trying to contourplot deviation but it is not giving output or error

I am trying to have contourplots of deviation dc for variables a and b.

In[27]:= \[Alpha]1[l_, b_, a_, \[Theta]_, M_,
r_] := -(a^2*(1 + l)*(2*M + 2*r + b) +
r*(2*r^2 + 3*b*r + b^2 - 2*M*(3*r + b)))/((a*
Sqrt[1 + l]*(2*M - 2*r - b))*Sin[\[Theta]])

In[28]:= \[Beta]1[l_, b_, a_, \[Theta]_, M_, r_] :=
Sqrt[(r^2*(8*a^2*(1 + l)*
M*(2*r + b) - (2*r^2 + 3*b*r + b^2 - 2*M*(3*r + b))^2)/((1 +
l)*a^2*(2*M - 2*r - b)^2)) +
a^2*Cos[\[Theta]]^2*(1 + l) - \[Alpha]1[l, b, a, \[Theta], M, r]^2*
Cos[\[Theta]]^2]

In[29]:= re1[l_, b_, a_, \[Theta]_, M_] :=
RankedMax[
r /. NSolve[{\[Beta]1[l, b, a, \[Theta], M, r] == 0, r > 0}, r,
Reals], 1]

In[30]:= pr1[l_, b_, a_, \[Theta]_, M_] :=
RankedMax[
r /. NSolve[{\[Beta]1[l, b, a, \[Theta], M, r] == 0, r > 0}, r,
Reals], 2]

In[31]:= xc1[l_, b_, a_, \[Theta]_, M_] :=
NIntegrate[\[Alpha]1[l, b, a, \[Theta], M, r]*\[Beta]1[l, b,
a, \[Theta], M, r]*D[\[Alpha]1[l, b, a, \[Theta], M, r], r], {r,
pr1[l, b, a, \[Theta], M], re1[l, b, a, \[Theta], M]}]/
NIntegrate[\[Beta]1[l, b, a, \[Theta], M, r]*
D[\[Alpha]1[l, b, a, \[Theta], M, r], r], {r,
pr1[l, b, a, \[Theta], M], re1[l, b, a, \[Theta], M]}]

In[32]:= R[l_, b_, a_, \[Theta]_, M_,
r_] := (\[Alpha]1[l, b, a, \[Theta], M, r] -
xc1[l, b, a, \[Theta], M])*
D[\[Beta]1[l, b, a, \[Theta], M, r], r] - \[Beta]1[l, b,
a, \[Theta], M, r]*D[\[Alpha]1[l, b, a, \[Theta], M, r], r]

In[33]:= raverage[l_, b_, a_, \[Theta]_, M_] :=
Sqrt[(1/\[Pi])*
NIntegrate[
R[l, b, a, \[Theta], M, r], {r, re1[l, b, a, \[Theta], M],
pr1[l, b, a, \[Theta], M]}]]

In[34]:= R1[l_, b_, a_, \[Theta]_, M_,
r_] := ((Sqrt[(\[Alpha]1[l, b, a, \[Theta], M, r] -
xc1[l, b, a, \[Theta], M])^2 + \[Beta]1[l, b, a, \[Theta],
M, r]^2] -
raverage[l, b, a, \[Theta],
M])^2)*(((\[Alpha]1[l, b, a, \[Theta], M, r] -
xc1[l, b, a, \[Theta], M])*
D[\[Beta]1[l, b, a, \[Theta], M, r], r] - \[Beta]1[l, b,
a, \[Theta], M, r]*
D[\[Alpha]1[l, b, a, \[Theta], M, r],
r])/((\[Alpha]1[l, b, a, \[Theta], M, r] -
xc1[l, b, a, \[Theta], M])^2 + \[Beta]1[l, b, a, \[Theta], M,
r]^2))

In[41]:= dc1[l_, b_, a_, \[Theta]_,
M_] := (1/raverage[l, b, a, \[Theta], M])*
Sqrt[(1/\[Pi])*
NIntegrate[
R1[l, b, a, \[Theta], M, r], {r, re1[l, b, a, \[Theta], M],
pr1[l, b, a, \[Theta], M]}]];

In[43]:= ContourPlot[
ConditionalExpression[dc1[0, b, a, 17*\[Pi]/180, 1],
Im[dc1[0, b, a, 17*\[Pi]/180, 1]] == 0], {b, 0, .42^2}, {a, .001, 1}]

Out[43]= \$Aborted


The code calculates the average radius of a parametric curve and then find the deviation of the curve from circularity in terms of rms value. the last line is not giving any output neither giving any error.

• Welcome to Mathematica.SE! I hope you will become a regular contributor. To get started, 1) take the introductory tour now, 2) when you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge, 3) remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign, and 4) give help too, by answering questions in your areas of expertise. – bbgodfrey Jun 13 at 0:09
• Your code simply may be slow. Have you tried evaluating dc1 for typical values of a and b to see whether the results are reasonable and the computations not too slow? Repeated calls to NIntegrate could consume a lot of time. – bbgodfrey Jun 13 at 0:13

Clear["Global*"]

α1[l_, b_, a_, θ_, M_,
r_] := -(a^2*(1 + l)*(2*M + 2*r + b) +
r*(2*r^2 + 3*b*r + b^2 - 2*M*(3*r + b)))/((a*
Sqrt[1 + l]*(2*M - 2*r - b))*Sin[θ])

β1[l_, b_, a_, θ_, M_, r_] :=
Sqrt[(r^2*(8*a^2*(1 + l)*
M*(2*r + b) - (2*r^2 + 3*b*r + b^2 - 2*M*(3*r + b))^2)/((1 + l)*
a^2*(2*M - 2*r - b)^2)) +
a^2*Cos[θ]^2*(1 + l) - α1[l, b, a, θ, M, r]^2*
Cos[θ]^2]

re1[l_, b_, a_, θ_, M_] :=
RankedMax[
r /. NSolve[{β1[l, b, a, θ, M, r] == 0, r > 0}, r, Reals], 1]

pr1[l_, b_, a_, θ_, M_] :=
RankedMax[
r /. NSolve[{β1[l, b, a, θ, M, r] == 0, r > 0}, r, Reals], 2]


Since xc1 uses a numeric technique, restrict its arguments to being numeric.

xc1[l_?NumericQ, b_?NumericQ, a_?NumericQ, θ_?NumericQ, M_?NumericQ] :=
NIntegrate[α1[l, b, a, θ, M, r]*β1[l, b, a, θ, M,
r]*D[α1[l, b, a, θ, M, r], r], {r,
pr1[l, b, a, θ, M], re1[l, b, a, θ, M]}]/
NIntegrate[β1[l, b, a, θ, M, r]*
D[α1[l, b, a, θ, M, r], r], {r, pr1[l, b, a, θ, M],
re1[l, b, a, θ, M]}]


Define R using Set rather than SetDelayed so that the derivatives are only done once.

R[l_, b_, a_, θ_, M_,
r_] = (α1[l, b, a, θ, M, r] - xc1[l, b, a, θ, M])*
D[β1[l, b, a, θ, M, r], r] - β1[l, b, a, θ, M,
r]*D[α1[l, b, a, θ, M, r], r];


Since raverage uses a numeric technique, restrict its arguments to being numeric.

raverage[l_?NumericQ, b_?NumericQ, a_?NumericQ, θ_?NumericQ,
M_?NumericQ] :=
Sqrt[(1/π)*
NIntegrate[
R[l, b, a, θ, M, r], {r, re1[l, b, a, θ, M],
pr1[l, b, a, θ, M]}]]

R1[l_?NumericQ, b_?NumericQ, a_?NumericQ, θ_?NumericQ, M_?NumericQ,
r_] := ((Sqrt[(α1[l, b, a, θ, M, r] -
xc1[l, b, a, θ, M])^2 + β1[l, b, a, θ, M,
r]^2] -
raverage[l, b, a, θ,
M])^2)*(((α1[l, b, a, θ, M, r] -
xc1[l, b, a, θ, M])*
D[β1[l, b, a, θ, M, r], r] - β1[l, b, a, θ,
M, r]*
D[α1[l, b, a, θ, M, r],
r])/((α1[l, b, a, θ, M, r] -
xc1[l, b, a, θ, M])^2 + β1[l, b, a, θ, M, r]^2))

dc1[l_?NumericQ, b_?NumericQ, a_?NumericQ, θ_?NumericQ,
M_?NumericQ] := (1/raverage[l, b, a, θ, M])*
Sqrt[(1/π)*
NIntegrate[
R1[l, b, a, θ, M, r], {r, re1[l, b, a, θ, M],
pr1[l, b, a, θ, M]}]]


It takes more than 5 seconds to calculate a single value for the plot.

ConditionalExpression[dc1[0, b, a, 17*π/180, 1],
Im[dc1[0, b, a, 17*π/180, 1]] == 0] /. {b -> 0,
a -> 0.001} // AbsoluteTiming

(* {5.47811, 1.23194*10^-7} *)


To control the number of calculations, generate a Table of data and use ListContourPlot

n = 9;

(plotData = Flatten[Table[
{b, a, ConditionalExpression[
dc1[0, b, a, 17*π/180, 1],
Im[dc1[0, b, a, 17*π/180, 1]] == 0]},
{b, 0, (42/100)^2, (42/100)^2/n},
{a, 1/1000, 1, (999/1000)/n}], 1];) // AbsoluteTiming

(* {635.234, Null} *)

ListContourPlot[plotData,
FrameLabel -> (Style[#, 14, Bold] & /@ {b, a}),
PlotLegends -> Automatic]


EDIT: Re "Can't I use ContourPlot so that the plot becomes continuous coloured instead of lines separating them? " Look at the documentation for ContourPlot. ContourPlot looks just like ListContourPlot, i.e., contour lines with ContourShading. Perhaps you meant to ask about DensityPlot or ListDensityPlot.

ListDensityPlot[
plotData,
FrameLabel -> (Style[#, 14, Bold] & /@ {b, a}),
PlotLegends -> Automatic]


• When I define R using set , the script is running . It doesn't end – Skj Jun 13 at 21:14
• Can't I use ContourPlot so that the plot becomes continuous coloured instead of lines separating them? – Skj Jun 13 at 21:31
• thanks for your answer and reply. But when I define R using set the script is running. It doesn't stop – Skj Jun 14 at 0:01
• You say you changed R`; did you run all of the code I posted above? What version are you running? I used v12.3 to produce the posted results. – Bob Hanlon Jun 14 at 0:06
• Now the code is running fine. Thank you Bob – Skj Jun 14 at 8:28