How to do a Numerical Integration of an Array Valued Function Handle with an Array as a limit?

Hi there, I have a little trouble with the following numerical integration problem. At the moment I have created two coefficient matrices A and B which are applied to the function handle fi. Both matrices A and B are of size n x n x n. The function is of the sort:

f = @(t) cos(A*t) + sin(B + sin(2*A*t))

The real function is more complex but I think this is sufficient for the problem I'm having. This function should be integrated, where the upper limit is stored in A at each corresponding position. I already tried arrayfun.

arrayfun(@(x) integral(f, 0, x, 'ArrayValued', true), A, 'UniformOutput', false)

But I don't think that it does what I want it to do, which is simply to numerically integrate for each individual instance of the function in the array with the corresponding upper limit.

It is maybe helpful to state where I'm coming from with this problem. The core problem is that I'm having a function that uses 3 for-loops to accomplish what I'm trying to do without any loops:

for i = 1:h
      for j = 1:t
          theta_s = theta(j);
          for k = 1:v,
              To            = aux(i,j) * aux2(k);
              result        = integral(f, 0, To);
              Integ2(i,j,k) = aux3(k,j) * result / Integ1(j);
          end
      end
  end

Here theta, aux, aux2, aux3 and Integ1 are arrays, which are passed to the function. To and theta_s are declared as global variables, so that the function f always uses the correct values for the evaluation of the integral.

Since the whole function is called in other for-loops, it takes forever to do it this way and I think there must be a way to solve this loop-free and therefore faster. Maybe I'm just to worked up in the problem and can't see the wood for the trees. Help is greatly appreciated.

Thanks in advance

Answer 1

[Mike Hosea]

I'm not sure it helps anything, but you don't need global variables for that. Instead, make f a function of both t and theta_s, and on To as well if you need that. The call site would look like

result        = integral(@(t)f(t,theta_s,To), 0, To);

Now if you integrand depends on To, then the only way to use 'ArrayValued' is to transform the problem so that each integral has the same upper bound. This might or might not be a good thing to do, depending on how efficient you can make the transformed integrand function.

But on the whole, "vectorizing" multiple integrations is a mixed bag. It seems like a good idea, and it is in some cases, but adaptive integration is all about putting extra evaluations of the integrand only where they are needed. If you try to compute simultaneously a bunch of integrals with things going on at different frequencies, you are at risk of just doing a lot of evaluations everywhere, which is wasted effort on most components but needed on some. It might be faster or it might be slower.

[Kevin Redeker]

Thank you for your response. The possibility to remove the global variable helps a bit actually. The only thing I can think of now, is to "hard code" the numerical integration in such a way that it fits my needs, but then again I don't know if it will be exectuted faster or not. Is there maybe any other way to solve this faster?

[Mike Hosea]

The ways of solving it faster that I can think of all involved transforming the integrand. Here's an example I just made up showing the brute force approach with the loop, the transformed approach, and the full-on arrayvalued approach.

function arrayval
f = @(x,v,t)cos(v*x + t*x^2)
v = (1:10).'
t = linspace(0,1,7)
a = pi/2;

% First way.  Integrate f from a to t(j).
tic
Q1 = zeros(length(v),length(t));
for j = 1:length(t)
    Q1(:,j) = integral(@(x)f(x,v,t(j)),a,t(j),'ArrayValued',true);
end
toc

% Intermediate approach.  Transform the problem so that the integrals are
% calculated over the interval from 0 to 1.
g = @(u,v,t,a)f(a + u*(t - a),v,t)*(t - a);
tic
Q2 = zeros(length(v),length(t));
for j = 1:length(t)
    Q2(:,j) = integral(@(u)g(u,v,t(j),a),0,1,'ArrayValued',true);
end
toc
norm(Q2(:) - Q1(:))

% Fully vectorized and transformed.
tic
Q3 = integral(@(u)h(u,v,t,a),0,1,'ArrayValued',true);
toc
norm(Q3(:) - Q1(:))

function y = h(u,v,t,a)
tma = t - a;
x = a + u*tma;
vx = bsxfun(@times,v,x);
tx2 = t.*x.^2;
fx = cos(bsxfun(@plus,vx,tx2));
y = bsxfun(@times,fx,tma);

[Kevin Redeker]

Hi, Mike

thank you very much for the help. This will be very useful for me. I came however to think of another way to solve it faster, since I noticed, that there are many duplicate entries in the boundary condition array. My quick and dirty solution is now to filter those out with the help of the "unique" command and have a much smaller problem size. But this will only work for now and I'm not sure this will be still the case, when I fit this program to the other conditions I'm investigating. So I'm very thankful, that you gave me this as another approach