1
$\begingroup$

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?

$\endgroup$
  • $\begingroup$ How are the variables r, fz, ft, fr and the function T defined? You seem to be using Octave – are you looking for Octave-compatible Matlab code or are Matlab-only options okay? $\endgroup$ – horchler Sep 11 '14 at 23:47
  • $\begingroup$ r, fz, ft, fr are constants. (Radius and forces applied in cylindrical coordinates). T is the result of solving a set of 8 ODE's with an initial condition at T(a) that is an 8x8 matrix so the solution of T at x, T(x), is also an 8x8 matrix. I am solving for T either numerically or with Eigen values and Eigen Vectors. Right now my preference is Octave as I've installed the beta version of Yosemite and having difficulties with my Matlab. $\endgroup$ – Andrew Sep 12 '14 at 0:14
0
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Yes, I should be using a typical linear solution method. I really like 2 and 3 as hopefully the set of linear equations does not have to be solved 8 times for each iteration. $\endgroup$ – Andrew Sep 12 '14 at 15:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.