definite 4D integration in octave, with functions as limits

Question

I need to do 4D integration in octave. My function is f(x,y,phi,theta) and some of the integration limits are function of the outer limits.

0 < theta < pi
t1(x,y) < phi < t2(x,y)
h1 < y < h2
w1 < x < w2

I wrote in octave like this (generalization):

[q1(i)] = integral( @(x) (integral3( @(y, phi, theta) f3(x, y, phi, theta), h1 , h2 , @(x,y) t1(x,y), @(x,y) t2(x,y), 0, pi)), w1, w2, 'ArrayValued',true);

my actual code:

clear all;
clc;
rho_bulk = 2.44; # rho_bulk = 2.44 uOhm.cm
h = 20e-9;
p = 0.5;
lambda = 40e-9;
n = 10;
w = linspace(20e-9,80e-9,n);

  for i = 1:n
  # limit for theta
  p2    = pi;
  p1    = 0;

  # limit for phi
  p4    = @(x,y) atan(x/(h-y)) + (pi/2);
  p3    = @(x,y) -atan((h-y)/(w(i)-x));

  # limit for y
  p6    = h;
  p5    = 0;

  # limit for x
  p8(i) = w(i);
  p7    = 0;

  #   f(x, y, phi, theta); outer --> inner
  #   limits;              inner --> outer

  f1 = @(x, y, phi, theta) exp(-(h-y)/(lambda *sin(theta) *sin(phi)));  
  f3 = @(x, y, phi, theta) sin(theta).*cos(theta).^2 .* f1(x, y, phi, theta);

  [q1(i)] = integral( @(x) (integral3( @(y, phi, theta) f3(x, y, phi, theta), p5, p6, @(x,y) p3(x,y), @(x,y) p4(x,y), p1, p2)), p7, p8(i), 'ArrayValued',true);

I have error from the integration line

error: 'y' undefined near line 51 column 98

I learned about the integration by following these:

https://www.mathworks.com/matlabcentral/answers/77571-how-to-perform-4d-integral-in-matlab

Quadruple Integral Using Nested Integral2 in Matlab

Answer 1

I believe the problem is in your definition of the limits in the call to integral3:

integral3( @(y, phi, theta) f3(x, y, phi, theta), p5, p6, @(x,y) p3(x,y), @(x,y) p4(x,y), p1, p2)

You're trying to integrate in y, phi and theta from p5 to p6, p3(x,y) to p4(x,y) and p1 to p2, respectively. integrate3 allows function-valued limits, but only in a very specific way:

q = integral3 (f, xa, xb, ya, yb, za, zb, prop, val, …)

    Numerically evaluate the three-dimensional integral of f using adaptive quadrature
    over the three-dimensional domain defined by xa,
    xb, ya, yb, za, zb (scalars may be finite or infinite). Additionally,
    ya and yb may be scalar functions of x and za, and zb maybe be scalar
    functions of x and y, allowing for integration over non-rectangular
    domains.

This actually reflects how integration would work on paper. So what you have is:

p5, p6,
@(x,y) p3(x,y), @(x,y) p4(x,y),
p1, p2

Your second limits try to depend on two variables, whereas they may only depend on one: the first integration variable. But this is fine, since x is not an integration variable, it's merely a parameter. So I believe the following should work:

integral3(@(y, phi, theta) f3(x, y, phi, theta), p5, p6, @(y) p3(x,y), @(y) p4(x,y), p1, p2)

By defining the limits as single-valued functions we basically curry your p3 and p4 functions to have a single-valued function that depends on the y integration variable.