Numerically integrate a function f(x) over x using MATLAB where f(x) has another argument y which is a vector

I would like to numerically integrate a vector which represents a function f(x) over the range of x specified by bounds x0 and x1 in Matlab. I would like to check that the output of the integration is correct and that it converges.

There are the quad and quadl functions that serve well in identifying the required error tolerance, but they need the input argument to be a function and not the resulting vector of the function. There is also the trapz function where we can enter the two vectors x and f(x), but then it computes the integral of f(x) with respect to x depending on the spacing used by vector x. However, there is no given way using trapz to adjust the tolerance as in quad and quadl and make sure the answer is converging.

The main problem why I can't use quad and quadl functions is that f(x) is the following equation: f(x) = sum(exp(-1/2 *(x-y))), the summation is over y, where y is a vector of length n and x is an element that is given each time to the function f(x). Therefore, all elements in vector y are subtracted from element x and then the summation over y is calculated to give us the value f(x). This is done for m values of x, where m is not equal to n.

When I use quadl as explained in the Matlab manual, where f(x) is defined in a separate function .m file and then in the main calling file, I use Q = quadl(@f,x0,x1,tolerance,X,Y); here X is a vector of length m and Y is a vector of length L. Matlab gives an error: "??? Error using ==> minus Matrix dimensions must agree." at the line where I define the function f(x) in the .m function file. f(x) = sum(exp(-1/2 *(x-y)))

I assume the problem is that Matlab treats x and y as vectors that should be of the same length when they are subtracted from each other, whereas what's needed is to subtract the vector Y each time from a single element from the vector X.

Would you please recommend a way to solve this problem and successfully numerically integrate f(x) versus x with a method to control the tolerance?

Answer 1

[Mike Hosea]

If you have R2012a,

integral(@(x)sum(exp(-0.5*(x-y))),x0,x1,'ArrayValued',true)

For example,

    >> y = 0:0.5:4;
    >> x0 = 1; x1 = 7;
    >> integral(@(x)sum(exp(-0.5*(x-y))),x0,x1,'ArrayValued',true,'AbsTol',1e-12,'RelTol',1e-10)
  ans =
    34.445963744387342

Even though this is scalar-valued, the ArrayValued option is useful because it prevents the integrator from trying to evaluate more than one x value at a time. This might be somewhat faster, and the technique works with QUADGK (and QUAD or QUADL, but don't use them).

    >> y = 0:0.5:4;
    >> x0 = 1; x1 = 7;
    >> f = @(x)arrayfun(@(t)sum(exp(-0.5*(t-y))),x);
    >> integral(f,x0,x1,'AbsTol',1e-12,'RelTol',1e-10)
  ans =
    34.445963744387342

[Matt Fig]

Mike, why "don't use" quadl?

[Mike Hosea]

QUAD, QUADL, QUADV, DBLQUAD, and TRIPLEQUAD are far less reliable and quite often slower and less accurate than the newer functions INTEGRAL, INTEGRAL2, and INTEGRAL3. In numerical integration nothing is always true, if you catch my meaning, but the new functions are usually better. Often they are substantially better, and rarely are the older functions better. I especially don't like the error control of the old functions. It's absolute error control, but it's a local error tolerance. There's no theoretical reason to expect that the global error is approximately bounded by the tolerance. The newer functions accept global tolerances. It's one reason, BTW, why it's sometimes hard to compare apples-to-apples. A given input tolerance is, in effect, taken more seriously by INTEGRAL or QUADGK than by QUAD, QUADL, or QUADV.

[Matt Fig]

Thanks Mike!

[Dina]

Thank you all very much. Unfortunately I don't have R2012a, I have R2011a, so I can't use the integral function.

However, I have posted this question on another forum and someone kindly provided me with the answer to my question using quadl.

The website is: stackoverflow.com/questions The person who provided the following answer is: macduff

From the documentationon quad it says:

The function y = fun(x) should accept a vector argument x and return a vector result y, the integrand evaluated at each element of x.

So every time we call the function, we need to evaluate the integrand at each given x.

Also, to parameterize the function call with the constant vector Y, I recommend an anonymous function call. There's a reasonable demo here. Here's how I implemented your problem in Matlab:

function Q = test_num_int(x0,x1,Y) Q = quad(@(x) myFun(x,Y),x0,x1); end

function fx = myFun(x,Y) fy = zeros(size(Y)); fx = zeros(size(x)); for jj=1:length(fx) for ii=1:length(Y) fy(ii) = exp(-1/2 *(x(jj)-Y(ii))); end fx(jj) = sum(fy); end end Then I called the function and got the following output:

Y = 0:0.1:1; x0 = 0; x1 = 1; Q = test_num_int(x0,x1,Y)

Q =

   11.2544
The inputs for the lower and upper bound and the constant array are obviously just dummy values, but the integral converges very quickly, almost immediately. Hope this helps!

Thank you all for your help. I really appreciate it.

Answer 2

[Matt Fig]

Mike's solution is essentially the same, but just to show you how to do what you want in your function M-file:

function [Y] = func_x(x)
% help here
y  = 0:0.5:4;  % Or however you were defining y.
Y = zeros(size(x));
for ii = 1:length(x)
    Y(ii) = sum(exp(-1/2 *(x(ii)-y)));
end

Then from the command line you could call:

>> quadl(@func_x,1, 7)

ans =

            34.4459638717806

[Dina]

Thanks a lot Matt, that's really helpful.

Best regards,