Why the NIntegrate can not give a result?

Question

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:

Clear["`*"]
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_] := 
  kc^2*vh*Total[
    Table[Exp[I*n*Total[(vKh[2] - vKh[1])*a3]], {n, 0, m}]];
f21[x_, y_] := 
  kc^2*vh*Total[
    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}];
Guu

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? — Michael E2, Aug 12 '18 at 00:59
@MichaelE2 Thank you for your reminder. I have re-edited the question. — user59546, Aug 12 '18 at 01:25
@Bill If you give the values of x and y, you can find Tuu can be carried out. — user59546, Aug 12 '18 at 01:30
@Bill But I have the answer of 1.551584515978786723934004674758063793199421595*10^10. — user59546, Aug 12 '18 at 01:35
You have constants that differ by dozens of orders of magnitude; numerical results will be catastrophically bad. Please, work in natural units, where all parameters are of order one. — AccidentalFourierTransform, Aug 12 '18 at 03:29
@AccidentalFourierTransform I have change the parameters to order one, but it does not work. — user59546, Aug 12 '18 at 04:16
Your re-edit is not really accurate: I get loads of error messages and you report none, even though my first question was whether you got errors. See http://i.stack.imgur.com/IcJWL.png -- I would say that probably the issues lie with poor programming practice. It's better to write functions in which all the dependencies appear in their arguments. For instance, LinearSolve should probably not be evaluated until the matrix-vector arguments are completely numeric. But you cannot control that with the way your code is written. — Michael E2, Aug 12 '18 at 14:27
@MichaelE2 Thank you very much! I'm sorry that the errors have not been added. — user59546, Aug 13 '18 at 04:08
@MichaelE2 I have write functions in which all the dependencies appear in their arguments. But there still appear a lot of errors. What I should do next? — user59546, Aug 13 '18 at 05:05
@MichaelE2 I have rewritten the functions in which all the dependencies appear in their arguments, and the above codes have been renewed. But there still appear a lot of errors. What I should do next? — user59546, Aug 13 '18 at 09:08
I think you missed my glancing remark about NumericQ: Try Tuu[x_?NumericQ, y_?NumericQ] := ... -- You may want to add it to any other function that might get called at the "top level" (like Tuu gets called by NIntegrate) and which requires x and y to be numbers for them to work without errors. — Michael E2, Aug 13 '18 at 11:05
@MichaelE2 Thank you very much! However, it will also takes a long time. Besides rewrite the Tuu, is there any other methods to optimize the codes? — user59546, Aug 14 '18 at 01:07
@MichaelE2 Besides the long time, the error "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." will also appear. — user59546, Aug 14 '18 at 01:37

score 4 · Answer 1 · answered Aug 12 '18 at 15:35

There are a few programming issues:

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.
You need to control the evaluation of certain parts of the code
You compute the solution the same linear problem in the code for Tuu eight times. It's not incorrect, but it wastes time.
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.
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},
   ku*Sin[x]*Tuu
   ];
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  *)

Why the NIntegrate can not give a result?

1 Answers1