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

83 views (last 30 days)
Kevin Redeker on 11 Feb 2015
Commented: Kevin Redeker on 16 Feb 2015
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.

Mike Hosea on 13 Feb 2015
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 on 16 Feb 2015
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