# How to make the plot smooth

There are some stranges point which cause the figure isn't smooth. I think that the reason is the method of numericial integration. Here is my code, it's pleasure that one can revise the drawbacks in the code.

zslist = Table[0.8999 - i/480 0.8999, {i, 0, 99}];
ϵ := 0.000001;
zc = 10^-2;
d = 4;
Do[z0[i_, d_] :=
zslist[[i]] - ((1 - zslist[[i]]^d) ϵ^2)/(2 zslist[[i]]);
z1[i_, d_] := -(((1 - zslist[[i]]^d) ϵ)/zslist[[i]]);
s[i, d] =
NDSolve[{z''[ρ] == -((
2 z[ρ]^3 Derivative[1][z][ρ]^2)/(
1 - z[ρ]^4)) + (-(2/z[ρ]) -
Derivative[1][z][ρ]/(ρ (1 - z[ρ]^4))) (1 -
z[ρ]^4 + Derivative[1][z][ρ]^2),
z[ϵ] == z0[i, d], z'[ϵ] == z1[i, d]},
z, {ρ, ϵ, 10}];
f[ρ_, i_, d_] := s[i, d][[1, 1]][[2]][ρ];
r[i_, d_] :=
FindRoot[f[ρ, i, d] == zc, {ρ,
s[i, d][[1, 1]][[2]][[1, 1, 2]]}][[1, 2]];,
{i, 1, Length[zslist]}]
TabzsR = Table[{r[i, d], z0[i, d]}, {i, 1, Length[zslist]}];
ListPlot[TabzsR, AxesOrigin -> {0, 0}, PlotRange -> All]
Clear[zsR]
zsR[x_] := Interpolation[TabzsR][x];
Rmax = r[3, d]
Rmin = r[99, d]
△l = (Rmax - Rmin)/360
Rlist = Table[Rmax - △l i, {i, 0, 359}];
Do[zn0[i_, d_] :=
zsR[Rlist[[i]]] - ((1 - zsR[Rlist[[i]]]^d) ϵ^2)/(
2 zsR[Rlist[[i]]]);
zn1[i_, d_] := -(((1 - zsR[Rlist[[i]]]^d) ϵ)/
zsR[Rlist[[i]]]);
sn[i, d] =
NDSolve[{
z''[ρ] == -((
2 z[ρ]^3 Derivative[1][z][ρ]^2)/(
1 - z[ρ]^4)) + (-(2/z[ρ]) -
Derivative[1][z][ρ]/(ρ (1 - z[ρ]^4))) (1 -
z[ρ]^4 + Derivative[1][z][ρ]^2),
z[ϵ] == zn0[i, d], z'[ϵ] == zn1[i, d]},
z, {ρ, ϵ, 10}];
fn[ρ_, i_, d_] := sn[i, d][[1, 1]][[2]][ρ];
,
{i, 1, Length[Rlist]}]
δ = 10^(-7);
Do[Dfn[i_,
d_] := -((fn[ρ, i,
d] /. {ρ -> Rlist[[i]] - δ}) - (
fn[ρ, i,
d] /. {ρ -> Rlist[[i]] - 1.001 δ} ))/(δ -
1.001 δ);
D2fn[i_,
d_] := -((D[
fn[ρ, i, d], ρ] /. {ρ ->
Rlist[[i]] - δ}) - (
D[fn[ρ, i, d], ρ] /. {ρ ->
Rlist[[i]] - 1.001 δ} ))/(δ - 1.001 δ);
fnp[ρ_, i_, d_] :=
zc + Dfn[i, d] (ρ - Rlist[[i]]) +
D2fn[i, d] (ρ - Rlist[[i]])^2/2 ;
rArea[i_?NumberQ, d_?NumberQ] :=
NIntegrate[ρ/fnp[ρ, i, d]^2 Sqrt[
1 + D[fnp[ρ, i, d], ρ]^2/(1 -
fnp[ρ, i, d]^4)], {ρ, Rlist[[i]] - δ,
Rlist[[i]]}, WorkingPrecision -> 30] +
NIntegrate[ρ/fn[ρ, i, d]^2 Sqrt[
1 + D[fn[ρ, i, d], ρ]^2/(1 -
fn[ρ, i, d]^4)], {ρ, 0, Rlist[[i]] - δ},
WorkingPrecision -> 30] - Rlist[[i]]/zc;, {i, 1, Length[Rlist]}]
TabrAreaR = Table[{Rlist[[i]], rArea[i, d]}, {i, 1, Length[Rlist]}];
ListPlot[TabrAreaR, AxesOrigin -> {0, 0}, PlotRange -> All,
Joined -> True]

• Your code doesn't produce a plot as written. It is also too long as it is, so it is too demanding to ask us to wade through it and find the errors for you. Please reduce it to a much smaller minimal working example so that we may be able to help you. Commented Jul 27, 2015 at 17:48
• Thanks for your attention. I think the reason is from the NIntegrate, there are many warning messages about NIntegrate. Commented Jul 28, 2015 at 4:03
• I can obtain the two plots. I recommend that you delete WorkingPrecision -> 30 from NIntegrate, because it does not increase accuracy and generates numerous warning messages. Because the integrands are smooth, it seems unlikely that the irregularities arise from NIntegrate. More likely, the irregularities are associated with the finite differences that feed into the integrands. Commented Jul 28, 2015 at 6:53
• How do you get the right plot? I delete WorkingPrecision -> 30 as you said, there are still the irregularity. And how to improve the inetgrand to decrease the feedback due to the finite differences? Commented Jul 28, 2015 at 9:33
• @amonxu It is difficult to try to improve your code without knowing what you are trying to accomplish. For instance, why are you using finite differences to construct Dfn instead of taking the second derivative? Commented Jul 28, 2015 at 12:50

In an effort to answer the question, I

• Rationalized all constants.
• Used WhenEvent to stop DSolvebefore z'[ρ] becomes singular.
• Solved for both z and z' in the second instance of DSolve to improve the accuracy of the integrand of the subsequent NIntegrate.
• Commented out the instance of NIntegrate involving fnp, which has no significant effect on the results.
• Tried Method -> "DoubleExponential" and ("LocalAdaptive" also) to improve results of NIntegrate.
• Experimented with increased WorkingPrecision.

Only the last of these had any substantive effect, and then not the expected one. Below is the code with WorkingPrecision not set.

ϵ = 10^-6; zc = 10^-2; d = 4; δ = 10^-7; imax = 100; imaxn = 360;
zslist = Table[(1 - i/480) 8999/10000, {i, 0, imax - 1}];

Do[z0[i_, d_] := zslist[[i]] - ((1 - zslist[[i]]^d) ϵ^2)/(2 zslist[[i]]);
z1[i_, d_] := -(((1 - zslist[[i]]^d) ϵ)/zslist[[i]]);
NDSolve[{z''[ρ] == -((2 z[ρ]^3 z'[ρ]^2)/(1 - z[ρ]^4)) + (-(2/z[ρ]) -
z'[ρ]/(ρ (1 - z[ρ]^4))) (1 - z[ρ]^4 + z'[ρ]^2),
z[ϵ] == z0[i, d], z'[ϵ] == z1[i, d],
WhenEvent[z[ρ] == zc, r[i, d] = ρ; "StopIntegration"]},
z, {ρ, ϵ, 10}], {i, 1, imax}];
TabzsR = Table[{r[i, d], z0[i, d]}, {i, 1, imax}];
ListPlot[TabzsR, PlotRange -> All, AxesLabel -> {"r", "z0"}]

zsR[x_] = Interpolation[TabzsR][x];
Rmax = r[3, d] ; Rmin = r[imax - 1, d] ; △l = (Rmax - Rmin)/imaxn ;
Rlist = Table[Rmax - △l i, {i, 0, imaxn - 1}];

Do[zn0[i_, d_] := SetPrecision[zsR[Rlist[[i]]] - ((1 - zsR[Rlist[[i]]]^d)
ϵ^2)/(2 zsR[Rlist[[i]]]), 15];
zn1[i_, d_] := SetPrecision[-(((1 - zsR[Rlist[[i]]]^d) ϵ)/zsR[Rlist[[i]]]), 15];
sn[i, d] = NDSolve[{z''[ρ] == -((2 z[ρ]^3 z'[ρ]^2)/(1 - z[ρ]^4)) + (-(2/z[ρ]) -
z'[ρ]/(ρ (1 - z[ρ]^4))) (1 - z[ρ]^4 + z'[ρ]^2),
z[ϵ] == zn0[i, d], z'[ϵ] == zn1[i, d],
WhenEvent[z[ρ] == zc, "StopIntegration"]},
{z, z'}, {ρ, ϵ, 10}];
fn[ρ_, i_, d_] := z[ρ] /. First@sn[i, d];
fnd[ρ_, i_, d_] := z'[ρ] /. First@sn[i, d], {i, 1, imaxn}]

Clear[rArea, fnp];
Do[(*fnp[ρ_,i_,d_]:=zc+fnd[Rlist[[i]]-δ, i, d] (ρ-Rlist[[i]])+(D[fnd[ρ, i, d],ρ]
/.ρ->Rlist[[i]]-δ) (ρ-Rlist[[i]])^2/2;*)
rArea[i_, d_] :=(* NIntegrate[ρ/fnp[ρ,i,d]^2 Sqrt[1+D[fnp[ρ,i,d],ρ]^2/
(1-fnp[ρ,i,d]^4)],{ρ, Rlist[[i]]-δ,Rlist[[i]]}]+*)
NIntegrate[ρ/fn[ρ, i, d]^2 Sqrt[1 + fnd[ρ, i, d]^2/(1 - fn[ρ, i, d]^4)],
{ρ, 0, Rlist[[i]](*-δ*)}, Method -> "DoubleExponential"]
- Rlist[[i]]/zc, {i, 1, imaxn}];
TabrAreaR = Table[{Rlist[[i]], rArea[i, d]}, {i, 1, imaxn}] ;
ListPlot[TabrAreaR, PlotRange -> All, AxesLabel -> {"Rlist", "rArea"}, Joined -> True]


This code produces the plots,

The irregularities described in the question are apparent in the second plot. Running the code with the commented items uncommented produces the same plots.

If the only issue were the noisy second plot, the irregularities could be removed by

ListPlot[MovingAverage[TabrAreaR, 20], Joined -> True,
PlotRange -> All, AxesLabel -> {"Rlist", "rArea"}]


However, while experimenting with the code, I found the second plot to be quite sensitive to seemingly minor changes. For instance, increasing the precision of zn0 and zn1 from 15 to 20 and adding WorkingPrecision ->20 to NIntegrate and the second instance of NDSolve yields for the second plot.

Thus, both instances of the second plot are not only noisy but, probably, wrong. The fundamental difficulty is this: fn approaches zero at the upper bound of NIntegrate, and as a consequence the integrand approaches infinity. A typical plot of fn is

Plot[fn[ρ, 215, 4], {ρ, 0, Rlist[[215]]}, AxesLabel -> {"ρ", "fn"}]


It is not difficult to show that fn varies as Sqrt[ρ0 - ρ] near where fn vanishes, designate ρ0. I leave it to someone more knowledgeable than I to obtain a final answer.