0
$\begingroup$

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]
$\endgroup$
7
  • $\begingroup$ 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. $\endgroup$
    – MarcoB
    Commented Jul 27, 2015 at 17:48
  • $\begingroup$ Thanks for your attention. I think the reason is from the NIntegrate, there are many warning messages about NIntegrate. $\endgroup$
    – amon xu
    Commented Jul 28, 2015 at 4:03
  • $\begingroup$ 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. $\endgroup$
    – bbgodfrey
    Commented Jul 28, 2015 at 6:53
  • $\begingroup$ 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? $\endgroup$
    – amon xu
    Commented Jul 28, 2015 at 9:33
  • $\begingroup$ @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? $\endgroup$
    – bbgodfrey
    Commented Jul 28, 2015 at 12:50

1 Answer 1

1
$\begingroup$

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,

enter image description here

enter image description here

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.

enter image description here

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"}]

enter image description here

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.

$\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.