# How to call an external function from within a numerical integration

Posted 8 years ago
6576 Views
|
8 Replies
|
0 Total Likes
|
 Dear Wofram community members,I have a problem with a code in which I call an external function from within an NIntegrate-command (numerical integration). I think that my problem has to do with me getting mixed up between local and global variables. In the attached file, I have tried to bring back the problem to its bare bones as best as I could. I'd be grateful for any suggestions what I'm doing wrong.Thanks in advance, René Samson Attachments:
8 Replies
Sort By:
Posted 8 years ago
 David,This explanation is very helpful. Thanks a heap. Sometimes Mathematica makes me think of the famous Einstein quote (in his refutation of the probabilistic interpretation of Quantum Mechanics): "God is subtle but he is not malicious". Something similar could be said about Mathematica, as illustrated by your explanation. Cheers, RenĂ©.
Posted 8 years ago
 Hi Rene,Here is what I meant by restricting the argument of the integrand to a numerical value. Now when NIntegrate evaluates it symbolically, it just gets the function back as the result, just like we would if we gave integrand[arg] a non-numerical argument. So then NIntegrate just proceeds to numerical computation with the function we've given it. It's good that NIntegrate behaves this way, because it allows it to adapt its strategy to the function, for example if it is piecewise, it knows to allow for discontinuity in choosing integration segments. But is does generate this issue when the argument it's given cannot be symbolically evaluated. In[1]:= pe = 25.; \[Epsilon] = 1.0/pe; psitrans = Sqrt[2./3.]; r[y_] := Sqrt[1 + y^2]; psifct[y_] := (y^2/2)*(1 - 3/(2 r[y]) + 1/(2 r[y]^3)); rhoapp[y_] := If[y < psitrans, (16*y/3)^(1/6), Sqrt[2 y]] // N; \[Psi][y_] := \[Epsilon]*y; In[6]:= (* integrand[arg] is only defined for numerical arg *) (* psiv,rhoappv,sol are localized within Module *) integrand[psiprime_?NumberQ] := Module[{psiv, rhoappv, sol}, psiv = \[Psi][psiprime]; rhoappv = rhoapp[psiv]; sol = FindRoot[psifct[x] - psiv == 0, {x, rhoappv}]; rhoprime = x /. sol; Exp[-rhoprime] ] In[7]:= integrand[a] Out[7]= integrand[a] In[8]:= integrand[10] Out[8]= 0.191518 In[9]:= (* So NIntegrate will not attempt to evaluate it symbolically *) (* edit -- it does attempt it, but gets the unevaluated function back, as stated above *) \ (* note that my use of psiprime for both the formal argument and the \ argument variable is for readability -- they are not related *) NIntegrate[integrand[psiprime], {psiprime, 0, 1000}] Out[9]= 11.6287 
Posted 8 years ago
 David,Thanks a lot. I didn't succeed in doing the _?NumberQ trick but turning the messages off will do as far as I'm concerned. You were very helpful. RenĂ©
Posted 8 years ago
 How to call an external function from within a numerical integration
Posted 8 years ago
 Dear David,Thanks for your response. I understand what you write. However, that was not my real problem. I just gave the simplified problem "One" in my Notebook as an example in order to understand better what I was doing. My real problem is under heading "Three" in the notebook I attached. Why doesn't that work while the syntax is identical to the code under my heading "Two" which does work? RenĂ©
Posted 8 years ago
 Sorry, Rene -- I read to the first issue. You will notice your code does produce a numerical solution, after the error messages. I suspect that FindRoot or NIntegrate is attempting to perform a symbolic evaluation before substituting the numerical value of psiprime, and objecting to tit's lack of a numerical value in the range. If this is the case, it is typically cured by specifying NumberQ as a conditional in the function definitions. But I suspect your solution is correct, in that once the symbolic evaluation fails (and generates the error) the computation proceeds with a numerical value.
Posted 8 years ago
 Indeed -- from the docs on NIntegrate: "NIntegrate symbolically analyzes its input to transform oscillatory and other integrands, subdivide piecewise functions, and select optimal algorithms."So NIntegrate is attempting to evaluate FindRoot with a symbolic value in the range specification. One way to deal with that is to live with (or turn off) the error messages. Another way is to wrap up the integrand in an external function with _?NumberQ used to restrict the arguments to numerical values, so NIntegrate knows it cannot use symbolic evaluation.Best, David
Posted 8 years ago
 In your notebook, you define f[x], but you then use f. f[x] is a definition for a function, where x is a formal parameter, so f as the function you defined must be called with a single argument. f with no argument has not been defined. In[1]:= Clear[x, y, f]; f[y_] := If[y < 1, y, y^2]; NIntegrate[f, {x, 0, 2}] During evaluation of In[1]:= NIntegrate::inumr: The integrand f has evaluated to non-numerical values for all sampling points in the region with boundaries {{0,2}}. >> Out[3]= NIntegrate[f, {x, 0, 2}] (* f is not defined *) f Out[1]= f (* f as a function with one argument is *) f[a] Out[5]= If[a < 1, a, a^2] In[4]:= NIntegrate[f[x], {x, 0, 2}] Out[4]= 2.83333