parameterized function of quadgk with another array input argument

Hello all,

I have been troubled with optimizing my numerical calculation using 'quadgk', and hope to get some advice on this. I will do my best to describe my problem, so let me know if anything is unclear.

1. there is a function of 3 input variables. i.e. f(x,y,z)

2. want to integrate over only one variable. (say z)

3. want to evaluate the whole function at several discrete values of remaining variables. i.e. at (x1,y1),(x2,y2)

It works fine when both x and y are scalar, but does not work when one of them is vector. Of course I can use the nested for loop to evaluate it at several (x,y), but I was wondering if there is more efficient way of doing this.

or example, putting x and y as array and getting f(x,y) as an array as well without going through loops.

ollowing is a simplified version of working code. (MATLAB 7.11.0)

function [w] = f(xx,yy)
    myfun = @(x,y) quadgk(@(z) exp(z).*(cos(x)+sin(y)),0,1);
    w = myfun(xx,yy);
end

So I want xx,yy to be arrays and evaluate myfun(xx,yy) without loops.

Any advice will be appreciated with this issue.

--------------------------------------------------------------

Below is an addition to the original question for further clarification.

As Teja suggested, ODE45 works fine with functions with simple form. But, my actual function has dependency to another function. For example,

myfun(x,y,z) = exp[f(y+z)+f(z+x)]

f(x) values are also calculated from another numerical calculation. With this function, I would like to integrate myfun(x,y,z) over z, at several (x,y) values.

Answer 1

[Teja Muppirala]

You could try to use an ODE solver instead of just a scalar integration function. For example.

%% Just make some random arrays of 1000 elements to loop over

xx = rand(1000,1);
yy = rand(1000,1);

%% Method 1. Use a loop

tic
I = zeros(size(xx));
for n = 1:numel(I)
    x = xx(n);
    y = yy(n);
    I(n) = integral( @(z) exp(z).*(cos(x)+sin(y)),0,1 ); 
end
toc;
plot(I);

%% Method 2. Use ODE45 and calculate the integral of the whole array at once

tic
J = ode45( @(z,unused) exp(z).*(cos(xx)+sin(yy)),[0 1],zeros(size(X)) );
J = J.y(:,end);
toc;
hold on;
plot(J,'r:');

%% Verify that you get the essentially same results.

max(abs(I - J))

ans =
 4.2831e-11

Just a couple of other observations,

1. If your integrand is really just exp(z).*(cos(x)+sin(y)) and not something more complicated, you could just pull the cos(x)+sin(y) out of the integral.

2. If you have a fairly new version of MATLAB, the recommended functions for integration are INTEGRAL, INTEGRAL2, and INTEGRAL3.

[James Lee]

Hello Teja,

Thank you for the advice! Although ODE45 was very impressive with performance compared to the loop, my function is actually not that simple... I realized I oversimplified the problem. I edited the question to better reflect my actual function. I would appreciate if you could take another look at the question.

Answer 2

[Alan Weiss]

I think you just need to make an anonymous function for your integral:

 @(z)myfun(x,y,z)

If myfun accepts matrix arguments, the function should work.

Alan Weiss

MATLAB mathematical toolbox documentation

[Asier Mariscal]

Dear Alan, could you please expand your answer? I have not been able to use quadgk when the function has a vector (say xx and yy above) as input as in the above example. I get an error saying matrix dimensions must agree I guess because quadgk also treats the integrating variable as a vector but of another dimension than xx). Is there an example you could write?Many thanks in advance

Answer 3

[Mike Hosea]

The virtues of avoiding loops in MATLAB are sometimes clear, sometimes unclear. QUADGK is a heavily vectorized integrator designed to deliver relatively high accuracy. Vectorizing the calling of QUADGK is not likely to help. IME, the differences in speed versus QUADGK come from two main sources. Both exist because QUADGK was designed to be robust. One is that QUADGK has a relatively large initial mesh. Often the integration is successful without refining a single interval. This sounds good, but it may also mean that a lot more arithmetic was done than necessary to solve the problem. The other source of performance differences is conservative global error estimation and control. If an integration is not being solved in one or two iterations, loosening the tolerances may speed things up.

If you really want to avoid loops, you can do it with arrayfun. Here I assume xx and yy come from (x,y) pairs (hence are arrays of the same size rather than vectors representing a grid)

function [w] = f(xx,yy)
    myfun = @(x,y)quadgk(@(z)f(x,y,z),0,1);
    w = arrayfun(myfun,xx,yy);
end