Numerically Integrate Matrix Equation in Matlab / Octave

Question

How do I integrate a matrix numerically in Octave / Matlab?

I am trying to do the following integration numerically in Matlab

$\int_a^b T(x)^{-1} B dx$

where $T(x)$ returns an 8x8 matrix. The vector B is a constant.

B = [0 0 0 0 -2*pi*r*fz -2*pi*r*ft -2*pi*r*fr 0]';
a = 5;
b = 20;

function y = f(x, a, B)
   y = T(a,x)^(-1)*B
endfunction

[q, ier, nfun, err] = quad(@(x) f(x,a,B), a, b)

I can see that y is being calculated as a 8x1 vector but q is returned as a scalar value.

Is there a suggested way for integrating an array in Octave or Matlab?

Answer 1

You essentially have eight independent quadrature problems (the rows of T). You have several options. Since you didn't provide runnable code, some of the code below may need to be tweaked. And if other issues crop up you'll have to deal with those.

1. Break up the problem and iterate over the rows using a for loop and quad (you might also try quadgk which performs better in many cases):

function y = f1(x,idx,a,B)
    % Integration function
    Tinv = T(a,x)^-1;
    y = Tinv(idx,:)*B; % Extract row
end

...
B = [0 0 0 0 -2*pi*r*fz -2*pi*r*ft -2*pi*r*fr 0]';
a = 5; b = 20;
fun = @(x,idx)f1(x,idx,a,B); % Anonymous function, pass a and B as parameters
q = zeros(8,1)
for i = 1:8
    % Iterate quadrature over each row
    [q(i),ier,nfun,err] = quad(@(x)fun(x,i), a, b);
end

2. Use vectorized quadrature via the quadv function:

function y = f2(x,a,B)
    % Integration function
    y = T(a,x)^-1*B;
end

...
B = [0 0 0 0 -2*pi*r*fz -2*pi*r*ft -2*pi*r*fr 0]';
a = 5; b = 20;
fun = @(x)f2(x,a,B); % Anonymous function, pass a and B as parameters
[q,nfun] = quadv(fun, a, b);

Note that the Matlab quadv function is scheduled to be removed in a future release.

3. Use Matlab and the newer integral function (uses adaptive Gauss-Kronrod quadrature like quadgk) with the 'ArrayValued' option:

...
B = [0 0 0 0 -2*pi*r*fz -2*pi*r*ft -2*pi*r*fr 0]';
a = 5; b = 20;
fun = @(x)f2(x,a,B); % Anonymous function, pass a and B as parameters
q = integral(fun, a, b, 'ArrayValued', true);

One final thing. Do you really need to be calculating an explicit inverse of a matrix? This is very rarely necessary. It leads to less numerical precision and, in some cases, numerical instability. You might look at you can modify your problem to use typical linear solution methods, e.g., \ (mldivide) or linsolve.