2
$\begingroup$

Usually it is advised to use NumericQ to avoid symbolic evaluation in NIntegrate. Here I'm trying to use NumericQ, but the NIntegrate gives an error:

t = 0.5;
ClearAll[vec];
vec[x_?NumericQ, t_?NumericQ] := {If[x > 1, 0, x Cos[t]],If[x > 1, 0, x Sin[t]]}
NIntegrate[v = vec[x, t]; Norm[v + {1, 1}], {x, 0, 2}]
(*NIntegrate: Integrand Sqrt[2] Abs[1+vec[x,1.2854]] is not numerical at {x} = {0.0159146}.*)

But after removing the NumericQ, NIntegrate works just fine:

ClearAll[vec];
vec[x_, t_] := {If[x > 1, 0, x Cos[t]], If[x > 1, 0, x Sin[t]]}
NIntegrate[v = vec[x, t]; Norm[v + {1, 1}], {x, 0, 2}]
(*3.31443*)

This is very confusing to me, isn't it that adding NumericQ fixes the non-numerical problem? But here it is totally opposite, removing the NumericQ fixes the non-numerical problem... Can someone explain what is going on here? Thanks.

$\endgroup$
  • 3
    $\begingroup$ The problem is that you're protecting only half of the integrand from being evaluated prematurely, thus breaking the whole thing. Try to evaluate Norm[vec[x, t] + {1, 1}] with both of your definitions for vec to see the difference $\endgroup$ – Lukas Lang Aug 7 at 7:34
  • $\begingroup$ @LukasLang, I have tried that, but changing the integrand to Norm[vec[x,t]+{1,1}] doesn't help. The function with NumericQ gives the same error, and the function without NumericQ works just fine. $\endgroup$ – JNL Aug 7 at 7:46
  • $\begingroup$ Sorry, I wasn't very clear - that's not what I meant. Execute the definition for vec without NumericQ, then evaluate Norm[vec[x, t] + {1, 1}] on its own. Now clear vec and execute the other definition, and evaluate Norm[vec[x, t] + {1, 1}] again $\endgroup$ – Lukas Lang Aug 7 at 7:48
  • $\begingroup$ Ok. I see what you meant. Thanks. $\endgroup$ – JNL Aug 7 at 7:57
  • $\begingroup$ If you want to have a composite integrand like v = vec[x, t]; Norm[v + {1, 1}], it's better to just wrap the whole thing into a new function and put that into NIntegrate instead. $\endgroup$ – Sjoerd Smit Aug 7 at 8:25
7
$\begingroup$

Expanding my comment into an answer for improved clarity.

The issue is that your only protecting half of your integrand from premature evaluation. In the case at hand, this breaks things because Plus behaves differently depending on the shape of its arguments. To see this, consider the following code:

t = 0.5;
ClearAll[vec];
vec[x_, t_] := {If[x > 1, 0, x Cos[t]], If[x > 1, 0, x Sin[t]]}
vec[x, t] + {1, 1}
% /. x -> 1
(* {
     1 + If[x > 1, 0, x Cos[0.5]],
     1 + If[x > 1, 0, x Sin[0.5]]
   } *)
(* {1.87758, 1.47943} *)

As you can see, everything works as expected - vec[x, t] is evaluated to {If[…],If[…]}, then {1,1} is added to the resulting vector. If we now insert a value for x, we get a vector, as it should be.

Now consider the definition involving NumericQ:

ClearAll[vec];
vec[x_?NumericQ, t_? NumericQ] := {If[x > 1, 0, x Cos[t]], 
  If[x > 1, 0, x Sin[t]]}


vec[x, t] + {1, 1}
% /. x -> 1
(* {1 + vec[x, 0.5], 1 + vec[x, 0.5]} *)
(* {{1.87758, 1.47943}, {1.87758, 1.47943}} *)

As you can see, the evaluation of vec[x, 0.5] is successfully prevented. The problem is that now, {1,1} is still added, but since vec[x, 0.5] is not a vector, it is added to each entry separately (see also a + {1,1}). Now, the result is no longer the same after inserting a value for x.

The fix

Knowing this, the fix is easy: Just protect the rest of the integrand from evaluation as well, using another function:

t = 0.5;
ClearAll[vec];
vec[x_, t_] := {If[x > 1, 0, x Cos[t]], If[x > 1, 0, x Sin[t]]}
func[x_?NumericQ, t_?NumericQ] := Norm[vec[x, t] + {1, 1}]
NIntegrate[func[x, t], {x, 0, 2}]

Note that we can even omit the NumericQ from vec now, since that is already taken care of by func.

Specifying the shape of a function

The question whether we can tell Mathematica somehow that a function will always return a vector with 2 entries or similar came up in the comments. I will present one a approach to achieve this below. I will try to outline a way to arrive at this approach, rather than just giving the solution (which can be found at the end). Note that for the specific problem in the question, I would use the fix above. This is more intended as an academic exercise, that may or may not be useful in some cases.

Let's take as an example the function f[x], which will return a vector with two elements:

ClearAll@f
f[x_?NumericQ] := {x, x^2}

This has the same problems as the vec function from the question:

{1, 1} + f[x]
(* {1 + f[x], 1 + f[x]} *)

We could try to make f[x] return something like {f[x][[1]], f[x][[2]]} to specify the shape explicitly:

ClearAll[f, if]
f[x_] := {if[x][[1]], if[x][[2]]}
if[x_?NumericQ] := {x, x^2}

We introduce the internal function if to distinguish between the one returning the shape and the one doing the actual computation. This has still some issues however:

f[2]
(* {2, 4} *)

f[x]
(* Part::partw: Part 2 of if[x] does not exist. *)
(* {x, if[x][[2]]} *)

The issue is that Part ([[…]]) blindly extracts items without waiting for a vector. This can be solved by using Indexed instead of Part:

ClearAll[f, if]
f[x_] := {Indexed[if[x], 1], Indexed[if[x], 2]}
if[x_?NumericQ] := {x, x^2}

f[x]
(* {Indexed[if[x], {1}], Indexed[if[x], {2}]} *)

enter image description here

This already works in principle. It has the disadvantage however that if[x] is evaluated several times, once for each index. We can mitigate this somewhat by memoizing the results of if[x]:

ClearAll[f, if]
f[x_] := f[x] = {Indexed[if[x], 1], Indexed[if[x], 2]}
if[x_?NumericQ] := {x, x^2}

Applying it to the question at hand: We can now use this to fix the issue in the question:

t = 0.5;
ClearAll[vec, ivec];
ivec[x_?NumericQ, t_?NumericQ] := 
 ivec[x, t] = {If[x > 1, 0, x Cos[t]], If[x > 1, 0, x Sin[t]]}
vec[x_, t_] := Array[Indexed[ivec[x, t], ##] &, 2]
Norm[vec[x, t] + {1, 1}]
(* Sqrt[
Abs[1 + Indexed[ivec[x, 0.5], {1}]]^2 + 
 Abs[1 + Indexed[ivec[x, 0.5], {2}]]^2]

NIntegrate[Norm[vec[x, t] + {1, 1}], {x, 0, 2}] *)
(* 3.31443 *)

Here is a generalized version that will turn any definition into one that returns something of the specified shape:

Attributes[setShape] = {HoldFirst};
setShape[sym_Symbol[args___] := rhs_, shape_] := With[
  (* create new symbol and clear it (could also use Unique[sym] *)
  {isym = Symbol["i" <> SymbolName@sym]},
  ClearAll@isym;
  (* make the original function return something of the correct shape for all arguments *)
  sym[args2___] := Array[Indexed[isym[args2], {##}] &, shape];
  (* attach definition to isym and add memoization *)
  call : isym[args] := call = rhs
  ]

We can now redo the demonstration from above:

ClearAll@f
setShape[f[x_?NumericQ] := {x, x^2}, 2]
f[x]
f[2]
(* {Indexed[if[x], {1}], Indexed[if[x], {2}]} *)
(* {2, 4} *)

An example where f returns a $2\times 2$ matrix instead:

ClearAll@f
setShape[f[x_?NumericQ] := {{x, x^2}, {x^3, x^4}}, {2, 2}]
f[x]
f[2]
(* {{Indexed[if[x], {1, 1}], 
  Indexed[if[x], {1, 2}]}, {Indexed[if[x], {2, 1}], 
  Indexed[if[x], {2, 2}]}} *)
(* {{2, 4}, {8, 16}} *)
$\endgroup$
  • $\begingroup$ Is there a way to let Mathematica knows the intended shape of the output of the function? Like can we make Mathematica assumes that vec[x_?NumericQ,t_?NumericQ]:={If[x > 1, 0, x Cos[t]],If[x > 1, 0, x Sin[t]]} always output a list of 2 elements, so that vec[x, t] + {1, 1} does not expand into {1 + vec[x, 0.5], 1 + vec[x, 0.5]} but {1 + vec[x, 0.5][[1]], 1 + vec[x, 0.5][[2]]}? $\endgroup$ – JNL Aug 7 at 9:14
  • $\begingroup$ @JNL Your question above shows that you do not understand the answer by Lukas. I suggest to read "Operations on Scalars, Vectors, and Matrices". $\endgroup$ – Anton Antonov Aug 7 at 9:37
  • $\begingroup$ Thanks for the link. So there is no way for Mathematica to treat vec[x_?NumericQ,t_?NumericQ] as a list without feeding it a numerical value as argument? $\endgroup$ – JNL Aug 7 at 10:19
  • 1
    $\begingroup$ @JNL See my edit, where I present one possible approach to do something like this. Note however that it most likely comes with a slight performance penalty (depending on what you're doing), and I would try to use the other fix if at all possible $\endgroup$ – Lukas Lang Aug 7 at 11:19
5
$\begingroup$

Here a simple workaround (only to complete @LukasLang very comprehensive answer)

vec[x_, t_] := Boole[x <= 1] {x Cos[t], x Sin[t]}
NIntegrate[v = vec[x, 0.5]; Norm[v + {1, 1}], {x, 0, 2}]
(*3.31443*)
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.