0
$\begingroup$

Hello Mathematica Users,

I'm using ?NumericQ to solve an integral numerically and then fit the data to give me values for certain parameters.

The integrals involved appear in:

$$ \alpha = \frac{\int_{0}^{r_{csda}} \rho\times e^{\mu\rho-\frac{2r_g}{\lambda}sin^{-1}(\rho/2r_g)}d\rho}{\int_{0}^{r_{csda}} \rho\times e^{\mu\rho-\rho/\lambda}d\rho} $$

where $r_g = \frac{p}{qB}$. In all this, $\alpha$ is the dependent variable, $B$ is the independent variable and the parameters I need to fit for are $\lambda$ and $r_{csda}$.

The numerator can't be solved analytically and so I have written the following Mathematica code:

ClearAll["Global`*"]
E0 := 10                                                   
(* \
Electron energy in MeV*)
(*St :=2.149;                (* 10 MeV 
Stopping power in water MeV/cm*)*)
r0 := 4.975 ;                                   
(* 10 MeV csda range \
in water cm*)
mc2 := 0.51099895 ;                      
(* Electron mass*c2 MeV*)
q := 3 ;                                                
(* Electron \
Charge compatible with units of Mev, 
cm and teslas *)
mu = 0.0157948384752221;                                                
(*attenuation constant*)

Zf[l_, csda_] = l*(l + (-csda + 
    l*(-1 + mu*csda)*Exp[(mu - 1/l)*csda])/((-1 + l*mu)^2));

rg[B_] = E0/q*B; 

NZfprime[l_, rho_, B_] = rho*Exp[mu*rho - (2*rg[B]/l)*ArcSin[rho/(2*rg[B])]];

alpha[l_?NumericQ, csda_?NumericQ, B_?NumericQ] := NIntegrate[NZfprime[l, rho, B], {rho, 0, csda}]/Zf[l, csda];

data = {{0, 1}, {0.3, 0.999552770379519}, {0.4, 
0.998805851062253}, {0.5, 0.998066565041057}, {0.6, 
0.997060919995234}};

solution1 = FindFit[data, {alpha[l,csda, B], l > 0, 
csda < 10}, {{l, 0.000001}, {csda, 2}}, {B}]

Show[Plot[solution1[x], {x, 0.1, 1}], ListPlot[data]]

The first problem is the data point {0,1}. Using my code, Mathematica doesn't recognise that when $B=0$, $\alpha =1$. However, when I eliminate that data point, Mathematica still has problems, probably because $\lambda$ is very small. I keep getting the following error messages:

NIntegrate::izero: Integral and error estimates are 0 on all integration subregions. Try increasing the value of the MinRecursion option. If value of integral may be 0, specify a finite value for the AccuracyGoal option. >>    
General::stop: Further output of NIntegrate::izero will be suppressed during this calculation. >>

and

NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small. >>
NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 9 recursive bisections in rho near {rho} = {1.789}. NIntegrate obtained -1.63470314044198918853761342200621303622880363340610336427127743697*10^-3960980+2.57268143657927767936334415959354139897206405201367090837780304657*10^-3960981 I and 8.5838860214619239726652239787055605419565681362037772502604174771`65.954589770191*^-3960980 for the integral and error estimates. >>

I'm only an occasional user of Mathematica and so do not have a really clear idea of what these error messages are prompting me to do. Does anybody have any solutions/insights?

Thanks in advance, RL

$\endgroup$
  • $\begingroup$ ArcSin-argument rho/(2*rg[B])==(3 rho)/(20 B) is singular for B==0 ? In this case ArcSin isn't real! $\endgroup$ – Ulrich Neumann Jun 18 at 10:03
  • $\begingroup$ I see that ArcSin-argument rho/(2*rg[B])==(3 rho)/(20 B) is singular for B==0. However the point that Mathematica fails to notice is that the full term (2*rg[B]/l)*ArcSin[rho/(2*rg[B]) is not singular for B==0; instead it converges to rho/l. $\endgroup$ – luiz Jun 19 at 3:07
2
$\begingroup$

With corrected variable assignments

ClearAll["Global`*"]
E0 = 10(*Electron energy in MeV*)
(*St:=2.149;(*10 MeV Stopping power in water MeV/cm*)*)
r0 = 4.975;
(*10 MeV csda range in water cm*)
mc2 = 0.51099895;
(*Electron mass*c2 MeV*)
q = 3;
(*Electron Charge compatible with units of Mev,cm and teslas*)
mu = 0.0157948384752221;
(*attenuation constant*)

Zf[l_, csda_] :=l*(l + (-csda +l*(-1 + mu*csda)*Exp[(mu - 1/l)*csda])/((-1 +l*mu)^2));
rg[B_] := E0/q*B;

and changed argument inside ArcSin

NZfprime[l_, rho_, B_] :=rho*Exp[mu*rho - (2*rg[B]/l)*ArcSin[rho/(2)*(rg[B])]];

alpha[l_?NumericQ, csda_?NumericQ, B_?NumericQ] :=NIntegrate[NZfprime[l, rho, B], {rho, 0, csda}]/Zf[l, csda];

data = {{0, 1}, {0.3, 0.999552770379519}, {0.4,0.998805851062253}, {0.5, 0.998066565041057}, {0.6,0.997060919995234}};

the fit evaluates an optimum

FindFit[data, {alpha[l, csda, B](*,l>0,csda<10*)}, {l, csda}, B,Method -> "NMinimize"]
(*{csda -> 0.162831, l -> 0.934341}*)
$\endgroup$
  • $\begingroup$ Thanks, Ulrich, for taking the time to correct my code. However, there are still a few issues. As I pointed out in the above comment, your change to the argument of ArcSin is leads to a function different from the original, correct one. (2*rg[B]/l)*ArcSin[rho/(2*rg[B]) does indeed converge to rho/l though at the moment Mathematica does not realize it. How can I correct my code so that Mathematica can deal with what is only apparently a singularity? $\endgroup$ – luiz Jun 19 at 3:09
  • $\begingroup$ What a pity. I tried to substitute the intgration variable rho->rho B but the singlarity still remains. Just an observation: If I minimize Rest[data] I get a very poor fit of the data? Are you sure about the model? $\endgroup$ – Ulrich Neumann Jun 19 at 7:07
  • $\begingroup$ ..The integrand NZfprime[l, rho, B] is only real for B==0 ! Though you try to fit a complex modell to real data? $\endgroup$ – Ulrich Neumann Jun 19 at 7:23
  • $\begingroup$ My understanding is that ArcSin[x] will be real as long as -1<x<1 . In the present case it will be real for all values of rho<2*rg[B]. Unless there is a mistake that I haven't seen, I believe the integrand should be real. $\endgroup$ – luiz Jun 19 at 7:50
  • $\begingroup$ I agree, rho<20/3 B implies the restriction csda<20/3B $\endgroup$ – Ulrich Neumann Jun 19 at 8:13
0
$\begingroup$

As I mentioned in the comment above, there was a slight error in the code. I also found a way to deal with the removable singularity. Many thanks to Ulrich for guiding and spurring me on.

ClearAll["Global`*"]
E0 = 10(*Electron energy in MeV*)
(*St:=2.149;(*10 MeV Stopping power in water MeV/cm*)*)
r0 = 4.975;
(*10 MeV csda range in water cm*)
mc2 = 0.51099895;
(*Electron mass*c2 MeV*)
q = 3;
(*Electron Charge compatible with units of Mev,cm and teslas*)
mu = 0.0157948384752221;
(*attenuation constant*)

Zf[l_, csda_] :=l*(l + ((-csda +l*(-1 + mu*csda))*Exp[(mu - 1/l)*csda])/((-1 
+l*mu)^2));
arc[B_,rho_,l_] := (2*E0/(q*B*l))*ArcSin[rho/(2*E0/(q*B))];
NZfprime[l_, rho_, B_] :=rho*Exp[mu*rho - arc[B,rho,l]];
limit=rho*Exp[mu*rho - rho/l]; 

alpha[l_?NumericQ,csda_?NumericQ,B_?NumericQ]:= 
NIntegrate[If[B==0,limit,NZfprime[l, rho, B]],
{rho,0,csda},Method{"GlobalAdaptive",Method"MultipanelRule"}]/Zf[l,csda];

data = {{0, 1}, {0.3, 0.999552770379519}, {0.4,0.998805851062253}, {0.5, 
0.998066565041057}, {0.6,0.997060919995234}};

FindFit[data, {alpha[l, csda, B](*,l>0,csda<10*)}, {l, csda}, B,Method -> 
"NMinimize"]
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.