-2
$\begingroup$

I have the following code for a numerical integration. I don't get an error message but it gives me a clearly wrong value (it is not evaluation the integral). What do you think the problem is?

   Tc = 170 (* critical temprature*) ; T = 
   1.5 Tc (*temprature in that moment in terms of the critical \
   temprature*); a = 0.15  ; \[Zeta] = 1;(* \[Kappa]=1.05; *)T0 = 10.8; \
   fs = 11; b = 7; Y0 = 7.6; Z = 82 ;
   e = Sqrt[4 \[Pi]/137];
   Tf = 130 (*final temprature? (it cooles down)*);
   t0 = 0.125(* fm *);
   R = 5 (* fm *); RA = 7 (* fm *);
   \[Epsilon] = 1;
   \[Sigma] =(*0.37 (T/Tc)^2*) \[Zeta] 0.018 T /197; (* fm^-1 *)
  \[Tau] := t0/cosh[\[Eta]];
   xpf = 5 ;
  points = 2;

  i1[\[Phi]p_?NumericQ] := i2[\[Phi]p] =
   NIntegrate[ -Z*
    x ((3/2 \[Pi]*RA^3) Sqrt[
    RA^2 - x^2 + 
     b x Cos[\[Phi]p] +(*I think there was a sign problem here*)
     b^2/4])(*first part of the field*)((e^2/4 \[Pi]) Sinh[
    Y0] (xp Cos[\[Phi]] - 
     x Cos[\[Phi]p]) (\[Tau]^2 Sinh[Y0 - \[Eta]]^2 + xp^2 + x^2 - 
      2 xp x Cos[\[Phi] - \[Phi]p])^(-3/
      2) (\[Sigma] Sinh[
       Y0] Sqrt[\[Tau]^2 Sinh[Y0 - \[Eta]]^2 + xp^2 + x^2 - 
         2 xp x Cos[\[Phi] - \[Phi]p]]/2 + 
     1) Exp[\[Sigma] Sinh[Y0 - \[Eta]] Sinh[Y0] \[Tau]/
       2 - \[Sigma]/2 Sinh[
       Y0] Sqrt[\[Tau]^2 Sinh[Y0 - \[Eta]]^2 + xp^2 + x^2 - 
        2 xp x Cos[\[Phi] - \[Phi]p]]])(*second part of the \
     field*)((e^2/4 \[Pi]) Sinh[
    Y0] (xp Cos[\[Pi] - \[Phi]] - 
     x Cos[\[Phi]p]) (\[Tau]^2 Sinh[Y0 + \[Eta]]^2 + xp^2 + x^2 - 
      2 xp x Cos[\[Pi] - \[Phi] - \[Phi]p])^(-3/
      2) (\[Sigma] Sinh[
       Y0] Sqrt[\[Tau]^2 Sinh[Y0 + \[Eta]]^2 + xp^2 + x^2 - 
         2 xp x Cos[\[Pi] - \[Phi] - \[Phi]p]]/2 + 
     1) Exp[\[Sigma] Sinh[Y0 + \[Eta]] Sinh[Y0] \[Tau]/
       2 - \[Sigma]/2 Sinh[
       Y0] Sqrt[\[Tau]^2 Sinh[Y0 + \[Eta]]^2 + xp^2 + x^2 - 
        2 xp x Cos[\[Pi] - \[Phi] - \[Phi]p]]]), {x, 0, 5}];
     i3[x_?NumericQ] := 
     NIntegrate[i2[\[Phi]p], {\[Phi]p, -\[Pi]/2, \[Pi]/2}];
$\endgroup$
  • $\begingroup$ How do you know the result is wrong? What result are you expecting? $\endgroup$ – FalafelPita Aug 2 '17 at 2:50
0
$\begingroup$

That would be a lot easier to read if you were to gather repeated portions into auxiliary constants or functions. For example

  fun1[x_] := \[Tau]^2 Sinh[Y0 - \[Eta]]^2 + xp^2 + x^2 - 
     2 xp x Cos[\[Phi] - \[Phi]p];
  fun2[x_] := \[Tau]^2 Sinh[Y0 + \[Eta]]^2 + xp^2 + x^2 - 
     2 xp x Cos[\[Pi] - \[Phi] - \[Phi]p];
  const1 = Sinh[Y0] \[Tau]/2 - \[Sigma]/2 Sinh[Y0];
  const2 = (e^2/4 \[Pi]) Sinh[Y0];
  funMain[x_] := -Z*
    x ((3/2 \[Pi]*RA^3) Sqrt[
       RA^2 - x^2 + 
        b x Cos[\[Phi]p] +(*I think there was a sign problem here*)
        b^2/4]) (const2 (xp Cos[\[Phi]] - 
        x Cos[\[Phi]p]) (fun1[x])^(-3/
         2) (\[Sigma] Sinh[Y0] Sqrt[fun1[x]]/2 + 
        1) Exp[\[Sigma] Sinh[Y0 - \[Eta]] const1 Sqrt[
         fun1[x]]]) (const2 (xp Cos[\[Pi] - \[Phi]] - 
        x Cos[\[Phi]p]) (fun2[x])^(-3/
         2) (\[Sigma] Sinh[Y0] Sqrt[fun2[x]]/2 + 
        1) Exp[\[Sigma] Sinh[Y0 + \[Eta]] const1 Sqrt[fun2[x]]]), {x, 
   0, 5}];


i1[\[Phi]p_?NumericQ] := 
  i2[\[Phi]p] = NIntegrate[funMain[x], {x, 0, 5}];

i3[x_?NumericQ] := 
  NIntegrate[i2[\[Phi]p], {\[Phi]p, -\[Pi]/2, \[Pi]/2}];

Other potential problems include

  • You defined i3 as a function of x but you didn't use x on the RHS, so it is a constant with respect to x.
  • xp is used, but not given a value, nor treated as a variable (through pattern matching). Did you mean x*p? Or did you mean xpf, which does have a value?
$\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.