I try to calculate a numberical integration in 2 layer for loops(outer 21 times, inner 128 times):
evalf(int(erfc(sqrt(B*xgms1*x/(x*xgmi1+xgmn1)))/gamav*exp(-1*x/gamav),x=0..infinity)), i.e., the kenerl is erfc(ax/(cx+d))×exp(-ax). But when "gamav" > 100000, the calculation time is too long to bear. So I try to do some approximation and change the upper bound into 10*gamav, 30*gamav, 50*gamav and 100*gamav. But I got the following results:
The trend of results is bigger the upper bound smaller the results. Obviously, the results is wrong. When "gamav" is a smaller number, such as 1000, the four integration results is the same as when the upper bound is "infinity".