0
$\begingroup$

I have a function to numerically integrate with NIntegrate. This function can be written mathematically as

$$ F(a,b,c,y,z) = \sum_i a^i b^j c^k f_{ijk}(y,z)$$

where $a,b,c$ are undefined variable (i.e. constant ) and $f_{ijk}(y,z)$ are function depending on the integration variables $y,z$.

I want to integrate $F(a,b,c,y,z)$ doing something like

NIntegrate[F[a,b,c,y,z],{y,0,1},{z,0,1}]

without assign any numerical value to $a,b,c$ ( I know that NIntegrate gives an error if $a,b,c$ have not numerical values) . I did it manually, that is I just integrated separately the functions $f_{ijk}(y,z)$ and then I reconstructed the final function of $a,b,c$ but I would like to do it quicker.

$\endgroup$
2
$\begingroup$

It is always useful to provide a simple example. Here is an example of your function:

SeedRandom[0];
F = Total[Times @@@ (Power[{a,b,c,y,z}^#]& /@ RandomInteger[5, {10,5}])]

a^4 b c^2 y+a^5 b^3 c^4 y+a^3 b^3 c^4 y^2+a^5 b y^4+a^5 y z+a^5 b^4 c^4 y z+a^2 b c^2 y^3 z+a^4 b^5 c^5 y^2 z^2+a^4 b c^2 y^5 z^2+a^4 b c^3 z^3

One applicable function is Collect:

Collect[F,{a,b,c},NIntegrate[#,{y,0,1},{z,0,1}]&]
(* expand if you desire *)
Expand @ %

0.125 a^2 b c^2+0.333333 a^3 b^3 c^4+a^5 (0.25 +0.2 b+0.5 b^3 c^4+0.25 b^4 c^4)+a^4 (0.111111 b^5 c^5+b (0.555556 c^2+0.25 c^3))

0.25 a^5+0.2 a^5 b+0.125 a^2 b c^2+0.555556 a^4 b c^2+0.25 a^4 b c^3+0.333333 a^3 b^3 c^4+0.5 a^5 b^3 c^4+0.25 a^5 b^4 c^4+0.111111 a^4 b^5 c^5

Another idea is to use CoefficientRules and FromCoefficientRules:

(* get exponent -> coefficient rules *)
rules = CoefficientRules[F, {a, b, c}];

(* use NIntegrate on coefficients *)
rules[[All, 2]] = NIntegrate[rules[[All, 2]], {y, 0, 1}, {z, 0, 1}];

(* reconstruct polynomial *)
FromCoefficientRules[rules, {a,b,c}]

0.25 a^5+0.2 a^5 b+0.125 a^2 b c^2+0.555556 a^4 b c^2+0.25 a^4 b c^3+0.333333 a^3 b^3 c^4+0.5 a^5 b^3 c^4+0.25 a^5 b^4 c^4+0.111111 a^4 b^5 c^5

$\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.