# Numerical Integration

I have used the following code to evaluate an integral (val) numerically

mu = 0.0173262004;(*attenuation coefficient for E = ?*)
k = (mu - 0.00324543007)/(0.00324543007);
h = 150.;
sz = 50.;
eg = 0.7;
qx = 1.0;
u = 1.0;
the = (22.5 Pi)/180;
f1 = 1/(l^2 + z^2);
f2 = 1 + k mu (l^2 + z^2)^0.5;
f3 = Exp[-(z - h)^2/(2 sz sz)];
f4 = Exp[-(z + h)^2/(2 sz sz)];
f5 = Exp[-mu (l^2 + z^2)];
val = NIntegrate[
2 Pi l f1 (f2 (f3 + f4) f5), {l, 0, Infinity}, {z, 0, Infinity}]
dr[r_] := 0.0404 0.00324543007 eg qx val/((2 Pi)^0.5 u r sz the)
LogPlot[dr[r], {r, 0, 10000}]


Now, i will be using the parameter (sz) as a function of r. For example, sz=0.26*r^0.69; How to do this? Thanking you in advance

• Do you mean you wish to define sz as a function of r instead of a constant (as it is now)? – honeste_vivere Sep 29 '14 at 13:43
• Yes, exactly... – Anand Sep 30 '14 at 12:40
• Then define it as sz[r_]:=0.26*r^0.69 and when calling it within something else, use sz[r], replacing r with whatever variable is relevant (r_ is a pattern so the input need not be explicitly the letter r). – honeste_vivere Sep 30 '14 at 13:20

You can define your val now as a function of r, rather than as a constant. Then you will be able to operate with this function. I will show it within a short example, and you can then implement into your code by analogy. Consider an integral

sz=5;
int=NIntegrate[x^2*Exp[-sz*x^2], {x, 0, Infinity}]
(*   0.0396333  *)


Now let us replace this sonstant sz by a function sz = 0.25*r^0.69:

int[r_] := NIntegrate[x^2*Exp[-0.25*r^0.69 x^2], {x, 0, Infinity}];


Now one can do with this function whatever he needs. Let us plot it, for example:

Plot[int[r], {r, 1, 5}, AxesLabel -> {r, int}]


giving this:

Have fun!

Corrected per input from Michael E2

mu = 0.0173262004;
(*attenuation coefficient for E=?*)
k = (mu - 0.00324543007)/(0.00324543007);
h = 150;
sz[r_] = 0.26*r^0.69;
eg = .7;
qx = 1.;
u = 1.;
the = (22.5 Pi)/180;
f1 = 1/(l^2 + z^2);
f2 = 1 + k mu (l^2 + z^2)^0.5;
f3[r_] = Exp[-(z - h)^2/(2 sz[r] sz[r])];
f4[r_] = Exp[-(z + h)^2/(2 sz[r] sz[r])];
f5 = Exp[-mu (l^2 + z^2)];

val[r_?NumericQ] := NIntegrate[
Evaluate[2 Pi l f1 (f2 (f3[r] + f4[r]) f5)],
{l, 0, Infinity}, {z, 0, Infinity}]

dr[r_?NumericQ] :=
0.0404*0.00324543007 eg qx val[r]/((2 Pi)^(1/2) u r sz[r] the) //
Evaluate

logPlot = LogPlot[dr[r], {r, 200, 10000},
PlotRange -> All,
PlotPoints -> 51]


dr[200]


8.80707*10^-45

• dr[100] gives errors for me. Note that LogPlot effectively uses Block[{r = 100},..] etc. – Michael E2 Oct 1 '14 at 1:56
• The problem with dr[100] is that the parameter r_?NumericQ is localized and does not refer to the global variable r that appears in sz and its dependents. LogPlot uses Block to set r, which sets the global r. In dr[100], the local function parameter r is 100, but the global r remains an undefined variable. The proper fix is to make sz, f3, and f4 functions of r. There are shorter workarounds, but they're kludgy (e.g., dr[r0_?NumericQ] := Block[{r = r0}, <code>] and val := NIntegrate[...]). – Michael E2 Oct 4 '14 at 12:34