Message Boards Message Boards

0
|
7397 Views
|
8 Replies
|
0 Total Likes
View groups...
Share
Share this post:

How to call an external function from within a numerical integration

Posted 9 years ago

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:
POSTED BY: Rene Samson
8 Replies
Posted 9 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 BY: Rene Samson
Posted 9 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 BY: David Keith
Posted 9 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 BY: Rene Samson
Posted 9 years ago

How to call an external function from within a numerical integration

POSTED BY: Rene Samson
Posted 9 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 BY: Rene Samson
Posted 9 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 BY: David Keith
Posted 9 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 BY: David Keith
Posted 9 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
POSTED BY: David Keith
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract