0
$\begingroup$

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

Improve this question
$\endgroup$
3
  • $\begingroup$ Do you mean you wish to define sz as a function of r instead of a constant (as it is now)? $\endgroup$ Commented Sep 29, 2014 at 13:43
  • $\begingroup$ Yes, exactly... $\endgroup$ Commented Sep 30, 2014 at 12:40
  • $\begingroup$ 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). $\endgroup$ Commented Sep 30, 2014 at 13:20

2 Answers 2

Reset to default
1
$\begingroup$

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:

enter image description here Have fun!

Improve this answer
$\endgroup$
1
$\begingroup$

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]

enter image description here

dr[200]

8.80707*10^-45

Improve this answer
$\endgroup$
2
  • $\begingroup$ dr[100] gives errors for me. Note that LogPlot effectively uses Block[{r = 100},..] etc. $\endgroup$ Commented Oct 1, 2014 at 1:56
  • $\begingroup$ 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[...]). $\endgroup$ Commented Oct 4, 2014 at 12:34

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.