Integral of matrix determinant

I found that it seems the

integral

fucntion cannot integrate determinant. For example

f=@(x)det([1,x;0,2]);
integral(f,0,1)

does not return a proper result. Of course we can mannually take

f=@(x)(2*1-0*x);

but this can be hardly applied to larger matrices. Is there a solution to this problem?

Answer 1

[Torsten]

f=@(x)det([1,x;0,2]);
integral(f,0,1,'ArrayValued',true)

ans = 2.0000

[Zhenghao Yang]

I see. What about 2 variables? Seems not directly applicable to integral2.

[Torsten]

1d:

f = @(x) det([1,x;0,2]);
value = integral(@(X)arrayfun(@(x)f(x),X),0,1)

value = 2.0000

2d:

f = @(x,y)det([x 2*y^2;0 3*sin(x*y)]);
value = integral2(@(X,Y)arrayfun(@(x,y)f(x,y),X,Y),0,1,0,1)

value = 0.4756

3d:

f = @(x,y,z) det([x 2*y^2 z;0 3*cos(x*z) log(z+1);4*exp(y+z) cos(z) 1/(x+1)]);
value = integral3(@(X,Y,Z)arrayfun(@(x,y,z)f(x,y,z),X,Y,Z),0,1,0,1,0,1)

value = -13.6865

Answer 2

[Askic V]

I would do it symbolically:

syms x
f = det([1,x;0,2]);
Q = int(f,0,1)

Q = 

2

[Zhenghao Yang]

Will this method run quick enough when the matrix is considerably large?

[Askic V]

For sure this is slower approach, but you need to test it and see if this would be sufficient for your particular application.

Answer 3

[John D'Errico]

For this specific problem, I might just suggest that the determinant of an upper triangular matrix is just the product of the diagonal elements.

syms x
A = [1 x;0 2]

A = 

As such, the use of det or integral is wild overkill here, since the determinant is independent of the value of x.

det(A)

ans = 

2

Why does integral fail? Because it tries to pass in a list of points to evaluate the function at. And det is not vectorized. So this is not a problem of integral in reality, but of the interaction between integral and det. Integral insists on the function being vectorized. Det insists is it NOT vectorized. Failure!

The solution is to evaluate the determinant as a polynomial in symbolic form, and then integrate, or do as @Torsten has shown, to use the 'arrayvalued' option in integral.

help integral

 INTEGRAL  Numerically evaluate integral.
    Q = INTEGRAL(FUN,A,B) approximates the integral of function FUN from A
    to B using global adaptive quadrature and default error tolerances.
 
    FUN must be a function handle. A and B can be -Inf or Inf. If both are
    finite, they can be complex. If at least one is complex, INTEGRAL
    approximates the path integral from A to B over a straight line path.
 
    For scalar-valued problems the function Y = FUN(X) must accept a vector
    argument X and return a vector result Y, the integrand function
    evaluated at each element of X. For array-valued problems (see the
    'ArrayValued' option below) FUN must accept a scalar and return an
    array of values.
 
    Q = INTEGRAL(FUN,A,B,PARAM1,VAL1,PARAM2,VAL2,...) performs the
    integration with specified values of optional parameters. The available
    parameters are
 
    'AbsTol', absolute error tolerance
    'RelTol', relative error tolerance
 
        INTEGRAL attempts to satisfy |Q - I| <= max(AbsTol,RelTol*|Q|),
        where I denotes the exact value of the integral. Usually RelTol
        determines the accuracy of the integration. However, if |Q| is
        sufficiently small, AbsTol determines the accuracy of the
        integration, instead. The default value of AbsTol is 1.e-10, and
        the default value of RelTol is 1.e-6. Single precision integrations
        may require larger tolerances.
 
    'ArrayValued', FUN is an array-valued function when the input is scalar
 
        When 'ArrayValued' is true, FUN is only called with scalar X, and
        if FUN returns an array, INTEGRAL computes a corresponding array of
        outputs Q. The default value is false.
 
    'Waypoints', vector of integration waypoints
 
        If FUN(X) has discontinuities in the interval of integration, the
        locations should be supplied as a 'Waypoints' vector. Waypoints
        should not be used for singularities in FUN(X). Instead, split the
        interval and add the results from separate integrations with
        singularities at the endpoints. If A, B, or any entry of the
        waypoints vector is complex, the integration is performed over a
        sequence of straight line paths in the complex plane, from A to the
        first waypoint, from the first waypoint to the second, and so
        forth, and finally from the last waypoint to B.
 
    Examples:
        % Integrate f(x) = exp(-x^2)*log(x)^2 from 0 to infinity:
        f = @(x) exp(-x.^2).*log(x).^2
        Q = integral(f,0,Inf)
 
        % To use a parameter in the integrand:
        f = @(x,c) 1./(x.^3-2*x-c)
        Q = integral(@(x)f(x,5),0,2)
 
        % Specify tolerances:
        Q = integral(@(x)log(x),0,1,'AbsTol',1e-6,'RelTol',1e-3)
 
        % Integrate f(z) = 1/(2z-1) in the complex plane over the
        % triangular path from 0 to 1+1i to 1-1i to 0:
        Q = integral(@(z)1./(2*z-1),0,0,'Waypoints',[1+1i,1-1i])
 
        % Integrate the vector-valued function sin((1:5)*x) from 0 to 1:
        Q = integral(@(x)sin((1:5)*x),0,1,'ArrayValued',true)
 
    Class support for inputs A, B, and the output of FUN:
       float: double, single
 
    See also INTEGRAL2, INTEGRAL3, FUNCTION_HANDLE

    Documentation for integral
       doc integral

    Other uses of integral

       gpuArray/integral

And, to be honest, the help does not really state that this is how you can avoid the need for your objective function to be vectorized.

f=@(x) det([1 x;0 2]);
integral(f,0,1,'ArrayValued',true)

ans = 2.0000

[Zhenghao Yang]

Very clear explanations, appreciations!