How to solve the following four integration in MATLAB

Hello all, I am trying to solve the following expression involving four integrations in MATLAB but note getting it correctly.

This is how I had tried to code:

integrand = @(y, h, z1, z2) (1 / gamma(m_0)) * ...
    gammainc( ((u_1 ./ ...
    (a_1 * zeta_1 * y .* g_1_abs_sq .* h .* P_h1_2 - ...
    u_1 * a_2 * zeta_1 * y .*  g_1_abs_sq .* h .*P_h1_2 - ...
    u_1 * a_1 * zeta_1 *  g_1_abs_sq .* h .* z1 - ...
    u_1 * a_1 * zeta_1 *P_h1_2 * z2 .* g_1_abs_sq .* h - ...
    u_1 * a_2 * zeta_1 * g_1_abs_sq .* h .* z1 - ...
    u_1 * a_2 * zeta_1 * P_h1_2 * z2 .* g_1_abs_sq .* h)) / ...
    (Omega_0 / m_0)),m_0) .* ...
    (1/gamma(m_0))*((m_0/Omega_0)^m_0)*y^(m_0-1)*exp(-m_0*y/Omega_0).* 1/(2*pi*A_1).*lambda_1*exp(-lambda_1*z1).*lambda_2*exp(-lambda_2*z2);
outer_integral = @(z2) arrayfun(@(z2_val) integral(@(y,h,z1) integrand(y,h,z1,z2_val ), 0, y_max(z2_val)), 0:1000);

% Perform the integration
op = integral(outer_integral, 0, 1000);

I am not getting why I am getting such errors:

Not enough input arguments.

Error in Analytical_2>@(z1,z2)(P_h1_2*(a_1-u_1*a_2))./(u_1*a_1*z1+u_1*a_1*P_h1_2*z2+u_1*a_2*z1+u_1*a_2*P_h1_2*z2) (line 96)

(u_1 * a_1 * z1 + u_1 * a_1 * P_h1_2 * z2 + ...

Error in Analytical_2>@(z2_val)integral(@(y,h,z1)integrand(y,h,z1,z2_val),0,y_max(z2_val)) (line 132)

outer_integral = @(z2) arrayfun(@(z2_val) integral(@(y,h,z1) integrand(y,h,z1,z2_val ), 0, y_max(z2_val)), 0:1000);

Error in Analytical_2>@(z2)arrayfun(@(z2_val)integral(@(y,h,z1)integrand(y,h,z1,z2_val),0,y_max(z2_val)),0:1000) (line 132)

outer_integral = @(z2) arrayfun(@(z2_val) integral(@(y,h,z1) integrand(y,h,z1,z2_val ), 0, y_max(z2_val)), 0:1000);

Error in integralCalc/iterateScalarValued (line 314)

fx = FUN(t);

Error in integralCalc/vadapt (line 132)

[q,errbnd] = iterateScalarValued(u,tinterval,pathlen);

Error in integralCalc (line 75)

[q,errbnd] = vadapt(@AtoBInvTransform,interval);

Error in integral (line 87)

Q = integralCalc(fun,a,b,opstruct);

Error in Analytical_2 (line 136)

op = integral(outer_integral, 0, 1000);

Please note that I had intentionally not given other part of code which is used to produce the values of elements that are used in this integrand. Any help in this regard will be highly appreciated.

[Torsten]

Try "integralN":

https://uk.mathworks.com/matlabcentral/fileexchange/47919-integraln-m

[Sam Chak]

Hi @chaaru datta

@Mike Hosea has written a function to perform such integration. See example here:

https://www.mathworks.com/matlabcentral/answers/168330-how-to-do-numerical-multiple-integral-more-than-triple-by-using-matlab

[chaaru datta]

Thank u sir for ur response, but I am not getting why MATLAB cant solve such 4 integrals.

[Torsten]

In MATLAB, there exist functions that can directly be applied to 1d, 2d and 3d integrals (namely integral, integral2 and integral3). If you have a 4d integral, the existing MATLAB routines can be combined (e.g. 3d with 1d, 2d with 2d, 2d with two times 1d and four times 1d) to solve the higher-dimensional integral. And that's exactly what "integralN" does.

If you want to use your personal method to solve your integral (e.g. by taking "integral" four times), you can do so. But as you see from your code, you don't succeed. I recommended "integralN" because it simplifies the set-up of such calculations.

Answer 1

[Walter Roberson]

integral expects a function handle to a function that accepts a single input parameter. It normally passes a vector of values to the function handle, and expects a vector of results of the same size. The size of the vector that it passes in is variable.

You are first passing integral() a function handle to a function that accepts a single parameter, which is good. Your code does

outer_integral = @(z2) arrayfun(@(z2_val) integral(@(y,h,z1) integrand(y,h,z1,z2_val ), 0, y_max(z2_val)), 0:1000);

The arrayfun() part of it is good. But inside the arrayfun() you are passing integral() a function handle of a function that expects up to three input parameters. integral() cannot handle functions that expect three input parameters. integral3 on the other hand is able to accepts function handles that expect three input parameters. Thus, proper code would more likely be

outer_integral = @(z2) arrayfun(@(z2_val) integral3(@(y,h,z1) integrand(y,h,z1,z2_val ), 0, y_max(z2_val)), 0:1000);

[chaaru datta]

Thank u for ur response sir...but still its not clearing my doubt...In my previous research work , the analytical expression was where f_z(z) and f_Y(y) are PDF of exponential random variable and Z_lim = and I had written the MATLAB code for this double integration as follows:

% Define Z_lim a s a function of y
Z_lim = @(y) (gamma_th * ((A + beta_0^2 * y.^2) * (P_u1_lin / jj1) + 1 / jj1)) / (beta_0^2 * y * A_11);
% Define the integrand
integrand = @(y, z) (1-exp(-lambda*(gamma_th*((A + beta_0^2 * y.^2) * (P_u1_lin / jj1) + 1 / jj1)-beta_0^2 * y * z *A_11)/A_11)).* lambda .* exp(-lambda * y) .* lambda .* exp(-lambda * z);

% Define the outer integral with respect to y from 0 to infinity
outer_integral = @(y) arrayfun(@(y_val) integral(@(z) integrand(y_val, z), 0, Z_lim(y_val)), y);
% Perform the double integration
op = integral(outer_integral, 0, inf);

My query is if the above MATLAB code of double integration is working correctly and giving me proper results then why the MATLAB code for four integrals is not working.

[Torsten]

My query is if the above MATLAB code of double integration is working correctly and giving me proper results then why the MATLAB code for four integrals is not working.

Because calling integral only twice as done in your code gives you a code for two integrals, not four (1d + 1d = 2d).

As @Walter Roberson wrote: You need integral + integral3 (or integral2 + integral2 or integral2 + integral + integral or integral + integral + integral + integral).

integral(@(y,h,z1) integrand(y,h,z1,z2_val )

does not work because the function "integral" has one input argument for 1d-integration, not three.

[chaaru datta]

Can u pls suggest how should I solve this four integration?

[Walter Roberson]

outer_integral = @(z2) arrayfun(@(z2_val) integral3(@(y,h,z1) integrand(y,h,z1,z2_val ), LOWER_y, UPPER_y, LOWER_h, UPPER_h, LOWER_z1, UPPER_z1), z2);

where LOWER_* is the lowerbound for the given variable and UPPER_* is the upperbound for the given variable.

[Torsten]

@chaaru datta

For example:

f = @(x,y,z,u) x + y + z + u;
value_integral_num = integral(@(x)integral3(@(y,z,u)f(x,y,z,u),0,1,0,1,0,1),0,1,'ArrayValued',1)

value_integral_num = 2.0000

syms x y z u
f =  x + y + z + u;
value_integral_sym = int(int(int(int(f,x,0,1),y,0,1),z,0,1),u,0,1)

value_integral_sym = 

2

both compute the four-fold integral of f over [0,1]^4.

But you should try the numerical variant. The symbolic approach is slow and should only be used if there are analytical antiderivatives of your function with respect to the integration variables.

Are you sure that the expression in the denominator has no zeros in the region of integration (which would cause problems for the integration process) ?

    (a_1 * zeta_1 * y .* g_1_abs_sq .* h .* P_h1_2 - ...
    u_1 * a_2 * zeta_1 * y .*  g_1_abs_sq .* h .*P_h1_2 - ...
    u_1 * a_1 * zeta_1 *  g_1_abs_sq .* h .* z1 - ...
    u_1 * a_1 * zeta_1 *P_h1_2 * z2 .* g_1_abs_sq .* h - ...
    u_1 * a_2 * zeta_1 * g_1_abs_sq .* h .* z1 - ...
    u_1 * a_2 * zeta_1 * P_h1_2 * z2 .* g_1_abs_sq .* h)

[chaaru datta]

@Walter Roberson I had written in the following way by replacing the lower and upper limits. The limits of z2,z1 and h are from 0 to infinity but limit of y is from 0 to y_max wherein y_max depends on variables z1 and z2.

outer_integral = @(z2) arrayfun(@(z2_val) integral3(@(y,h,z1) integrand(y,h,z1,z2_val ), 0, y_max(z1,z2_val), 0, inf, 0, inf), z2);
op = integral(outer_integral, 0, inf);

here i use y_max(z1,z2_val) but it gives error as 'Unrecognized function or variable 'z1'.'

[chaaru datta]

@Torsten sir, thank u for ur well drafted example but I need to take into account limits of variables while integrating.

[Torsten]

Assuming y_max is somewhere defined as a function of z1 and z2:

ymax = @(z1,z2) ...;
outer_integral = @(h) arrayfun(@(h_val) integral3(@(z1,z2,y) integrand(y,h_val,z1,z2 ), 0, inf, 0, inf, 0, y_max), h);
op = integral(outer_integral, 0, inf);

or shorter:

y_max = @(z1,z2) ...;
op = integral(@(h)integral3(@(z1,z2,y)integrand(y,h,z1,z2),0, inf, 0, inf, 0, y_max),0, Inf, 'ArrayValued',1)

[chaaru datta]

@Torsten sir now the code is running but value of op is NaN....

I am having one more query , in the double integration, we did the following :

% Define the outer integral with respect to y from 0 to infinity

outer_integral = @(y) arrayfun(@(y_val) integral(@(z) integrand(y_val, z), 0, Z_lim(y_val)), y);

we wrote Z_lim as Z_lim(y_val), but now in case of four integation, we just wrote y_max...

[Torsten]

integral3 expects the limits of integration to be function handles, not evaluated function handles.

Take a look at the example "Integral over the Unit Sphere in Cartesian Coordinates" under

https://uk.mathworks.com/help/matlab/ref/integral3.html

[chaaru datta]

Ohh got it ...so we are using integral3 and thats why we are simply writing y_max

[chaaru datta]

@Torsten sir why u have taken variable h outside in the below statement ...

ymax = @(z1,z2) ...;
    outer_integral = @(h) arrayfun(@(h_val) integral3(@(z1,z2,y) integrand(y,h_val,z1,z2 ), 0, inf, 0, inf, 0, y_max), h);
op = integral(outer_integral, 0, inf);

Cant we take any other element ?

[Torsten]

Yes.

If you had used "integral2" instead of two times "integral" in your double integration case, you also would have used

Z_lim = @(y)...;
op = integral2(integrand,0,inf,0,Z_lim)

with Z_lim being the function handle, not the evaluated function handle.

Maybe it's interesting to see whether your code with 2-times "integral" and the above with "integral2" gives the same results.

[Torsten]

sir why u have taken variable h outside in the below statement ...

Cant we take any other element ?

No. It's necessary to use h in the 1d-integral. y_max depends on z1 and z2 in the upper limit - thus y, z1 and z2 should be used together in the 3d-integral. This leaves only h for the 1d-integral.

Further, it's necessary to use y as last argument to integral3 because only the lower and upper limits of the last dimension (here: y_max) can depend on both previous dimensions (z1 and z2):

Cited from the documentation of "integral3":

q = integral3(fun,xmin,xmax,ymin,ymax,zmin,zmax) approximates the integral of the function z = fun(x,y,z) over the region xmin ≤ x ≤ xmax, ymin(x) ≤ y ≤ ymax(x) and zmin(x,y) ≤ z ≤ zmax(x,y).

And I already see one possible reason why you get NaN: h = 0 causes the denominator of your argument for gammainc to become 0. Thus you get a division by 0.

[chaaru datta]

I am getting confused now, I will once again share my query here:

I am working on two seperate analytical expression let us call them expression A and expression B.

Expression A is: where and are PDFs of exponential random variable and .

The MATLAB code (which is working correctly) for evaluating expression (A) is as:

% Define Z_lim a s a function of y
Z_lim = @(y) (gamma_th * ((A + beta_0^2 * y.^2) * (P_u1_lin / jj1) + 1 / jj1)) / (beta_0^2 * y * A_11);
% Define the integrand
integrand = @(y, z) (1-exp(-lambda*(gamma_th*((A + beta_0^2 * y.^2) * (P_u1_lin / jj1) + 1 / jj1)-beta_0^2 * y * z *A_11)/A_11)).* lambda .* exp(-lambda * y) .* lambda .* exp(-lambda * z);  
% Define the outer integral with respect to y from 0 to infinity
outer_integral = @(y) arrayfun(@(y_val) integral(@(z) integrand(y_val, z), 0, Z_lim(y_val)), y);
% Perform the double integration
op = integral(outer_integral, 0, inf);

Next, the expression B is:

---- (B)

where and

My query is based on the style of expression (A), how to write MATLAB code for expression (B) ?

[Torsten]

My query is based on the style of expression (A), how to write MATLAB code for expression (B) ?

Why are you getting confused ?

I gave you the answer previously (some * are replaced by .* and / are replaced by ./ ):

integrand = @(y, h, z1, z2) (1 ./ gamma(m_0)) .* ...
    gammainc( ((u_1 ./ ...
    (a_1 .* zeta_1 .* y .* g_1_abs_sq .* h .* P_h1_2 - ...
    u_1 .* a_2 .* zeta_1 .* y .*  g_1_abs_sq .* h .* P_h1_2 - ...
    u_1 .* a_1 .* zeta_1 .*  g_1_abs_sq .* h .* z1 - ...
    u_1 .* a_1 .* zeta_1 .* P_h1_2 .* z2 .* g_1_abs_sq .* h - ...
    u_1 .* a_2 .* zeta_1 .* g_1_abs_sq .* h .* z1 - ...
    u_1 .* a_2 .* zeta_1 .* P_h1_2 .* z2 .* g_1_abs_sq .* h)) ./ ...
    (Omega_0 ./ m_0)),m_0) .* ...
    (1./gamma(m_0)).*((m_0./Omega_0).^m_0).*y.^(m_0-1).*exp(-m_0.*y./Omega_0).* 1./(2.*pi.*A_1).*lambda_1.*exp(-lambda_1.*z1).*lambda_2.*exp(-lambda_2.*z2);
y_max = @(z1,z2) P_h1_2 .* (a_1 - u_1 .* a_2) ./ ...
                 (u_1 .* a_1 .* z1 + u_1 .* a_1 .* P_h1_2 .* z2 + ...
                  u_1 .* a_2 .* z1 + u_1 .* a_2 .* P_h1_2 .* z2);
op = integral(@(h)integral3(@(z1,z2,y)integrand(y,h,z1,z2),0, inf, 0, inf, 0, y_max),0, Inf, 'ArrayValued',1)
% or
%outer_integral = @(h) arrayfun(@(h_val) integral3(@(z1,z2,y) integrand(y,h_val,z1,z2 ), 0, inf, 0, inf, 0, y_max), h);
%op = integral(outer_integral, 0, inf);

I assumed z1 = Z1 and z2 = Z2 in your mathematical description.

But as written, it will not work because h = 0 gives a division by zero because the denominator of the argument to the gammainc function becomes 0.

[chaaru datta]

Ok sir...