The integrand includes ku, Sin[x], and Tuu. ku is constant, x is one of the variables, and Tuu is functions of x. I have tested that when x and y are given, Tuu can be carried out. But why NIntegrate[ku*Sin[x]Tuu, {x, 0, Pi/2}, {y, 0, Pi/4}] give a result of itself (i.e. NIntegrate[kuSin[x]*Tuu, {x, 0, Pi/2}, {y, 0, Pi/4}])? How to solve the problem? Many thanks!

According to the suggestion of Michael E2, I have rewrite the codes. But there still appear a lot of errors.

The codes are as following:

m = 4;
vh = 4;
mu = 11;
delta = 8;
HBAR = SetPrecision[1.05457266*10^(-34), 50];
ME = SetPrecision[9.1093897*10^(-31), 50];
ELEC = SetPrecision[1.60217733*10^(-19), 50];
Kh = SetPrecision[0.211, 50];
vKh[1] = {0, 0, 0};
vKh[2] = {Kh, 0, 0};
kc = Sqrt[2*ME*ELEC/HBAR^2]*10^(-11);
ku = kc*Sqrt[mu + delta];
kd = kc*Sqrt[mu - delta];
a3 = {Pi/Kh, Pi/Kh, Sqrt[2]*Pi/Kh};
k[x_, y_] := {-ku*Sin[x]*Cos[y], -ku*Sin[x]*Sin[y], kz};
f11[x_, y_] := Total[(k[x, y] + vKh[1])^2] - ku^2;
f22[x_, y_] := Total[(k[x, y] + vKh[2])^2] - ku^2;
f12[x_, y_] := 
    Table[Exp[I*n*Total[(vKh[2] - vKh[1])*a3]], {n, 0, m}]];
f21[x_, y_] := 
    Table[Exp[I*n*Total[(vKh[1] - vKh[2])*a3]], {n, 0, m}]];
t[x_, y_] := {{f11[x, y], f12[x, y]}, {f21[x, y], f22[x, y]}};
slu[x_, y_] := 
  Select[kz /. NSolve[Det[t[x, y]] == 0, kz], 
   Re[#] >= 0 && Im[#] >= 0 &];
td1[x_, y_] := t[x, y] /. kz -> slu[x, y][[1]];
td2[x_, y_] := t[x, y] /. kz -> slu[x, y][[2]]
nu1[x_, y_] := Flatten[NullSpace[td1[x, y]]];
nu2[x_, y_] := Flatten[NullSpace[td2[x, y]]];
pi[x_, y_] := Transpose[{nu1[x, y], nu2[x, y]}];
pei := Table[If[q == 1, 1, 0], {q, 2}];
Tuu[x_, y_] := 
  Im[Conjugate[(LinearSolve[pi[x, y], pei][[1]]*nu1[x, y][[1]]*
          Exp[I*slu[x, y][[1]]*d] + 
         LinearSolve[pi[x, y], pei][[2]]*nu2[x, y][[1]]*
          Exp[I*slu[x, y][[2]]*d])]*(Dt[
        LinearSolve[pi[x, y], pei][[1]]*nu1[x, y][[1]]*
          Exp[I*slu[x, y][[1]]*d] + 
         LinearSolve[pi[x, y], pei][[2]]*nu2[x, y][[1]]*
          Exp[I*slu[x, y][[2]]*d], d]) + 
     Conjugate[(LinearSolve[pi[x, y], pei][[1]]*nu1[x, y][[2]]*
          Exp[I*slu[x, y][[1]]*d] + 
         LinearSolve[pi[x, y], pei][[2]]*nu2[x, y][[2]]*
          Exp[I*slu[x, y][[2]]*d])]*(Dt[
        LinearSolve[pi[x, y], pei][[1]]*nu1[x, y][[2]]*
          Exp[I*slu[x, y][[1]]*d] + 
         LinearSolve[pi[x, y], pei][[2]]*nu2[x, y][[2]]*
          Exp[I*slu[x, y][[2]]*d], d]) /. d -> 2*10^(-9)];
Guu := NIntegrate[ku*Sin[x]*Tuu[x, y], {x, 0, Pi/2}, {y, 0, Pi/4}];
  • $\begingroup$ Does it give an error? Does it take longer than you're willing to wait? Does it run out of memory? Does it cause a crash? $\endgroup$ – Michael E2 Aug 12 '18 at 0:59
  • $\begingroup$ @MichaelE2 Thank you for your reminder. I have re-edited the question. $\endgroup$ – user59546 Aug 12 '18 at 1:25
  • $\begingroup$ @Bill If you give the values of x and y, you can find Tuu can be carried out. $\endgroup$ – user59546 Aug 12 '18 at 1:30
  • $\begingroup$ @Bill But I have the answer of 1.551584515978786723934004674758063793199421595*10^10. $\endgroup$ – user59546 Aug 12 '18 at 1:35
  • $\begingroup$ @Bill I guess you should change yy to y. $\endgroup$ – user59546 Aug 12 '18 at 1:36

There are a few programming issues:

  1. It's better to write functions in which all the dependencies appear in their arguments. This requires a complete rewrite which I'm not going to do.

  2. You need to control the evaluation of certain parts of the code

  3. You compute the solution the same linear problem in the code for Tuu eight times. It's not incorrect, but it wastes time.

  4. The condition Im[#] >= 0 in the code for slu does not accommodate rounding error on zero imaginary part, if the error were a small negative number.

  5. Non-issue: Some will complain about the scale of the parameters (i.e. numbers of the order 10^-31). Often problems with numerics are due to the size of numbers getting out of hand and catastrophic round-off error. But not always. This has yet to prove to be an issue in this particular set of problems (see OP's other recent posts), and the scale of the numbers in the computations have fit comfortably within the range of double-precision floats.

First of all, there are all these errors that were not mentioned in the question:

Mathematica graphics

They arise because NSolve and probably LinearSolve should not be evaluated until the expressions involved are completely numeric (see advice on _?NumericQ). But that cannot be controlled with the way the code is written (see point 1 above). It can be see that the same errors come from evaluating Tuu when x and y do not have numeric values. If we give them numeric values first, then Tuu evaluates to a number without errors.

Mathematica graphics

The first error you get from NIntegrate[] originates in td[1]. This is easy to figure out in V11, because you can click on the three-dot button next to the error message and select "Show Stack Trace". It comes up formatted a little more nicely than shown below. We see in lines 2-4 some of the code for td[1] and in line five we see the function call.

Mathematica graphics

Below is one way to get NIntegrate to compute something, but given all the programming issues, I cannot guarantee the value is correct:

(* takes a long time -- didn't time it (maybe 10-20 minutes?) *)
myG[xx_?NumericQ, yy_?NumericQ] := Block[{x = xx, y = yy},
NIntegrate[myG[x, y], {x, 0, Pi/2}, {y, 0, Pi/4}]

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.

(*  1.8222*10^20  *)
| improve this answer | |

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.