Numerical Integration with functions as borders

Hi. I am new to this Forum and I have a question to which I did yet not find an appropriate answer. Maybe someone can help.

I would like to solve an triple integral which has a function as border (I am not sure if this is the right word for the boundaries of an integral, but let's assume it is). The border-function itself depends on the variable which will be integrated with the 2nd integral.

For instance, look at this integral:

My problem is not analyticly solvable (using symbolic variables ans int()-function), so I need a numeric solution. Working with quad() and function-handles as borders doesn't work either. So how to integrate?

Thank you for your help,

Manu

[Sean de Wolski]

Should it be h1(phi) to h2(phi)? Above looks like you're only integrating at a specific point.

[Manuel]

Hi Sean,

it should be from h1(phi) to h2(phi). Sorry, I made typo. I will update the graphic.

Thank you for that remark.

Answer 1

[Mike Hosea]

QUAD2D accepts function handles for limits. Take care that you make them "vectorized" (usually just use .*, .+, and ./ instead of *, +, and /).

>> h1 = @(phi)-2 + sin(phi); % or whatever
>> h2 = @(phi)2 - cos(phi); % or whatever
>> B = @(xi,b)quad2d(@(phi,x)(1./(h1(phi) - h2(phi))).*(1./(b-(xi-x).^2).^(3/2)),-pi,pi,h1,h2);
>> B(5,100)
ans =
  -0.010761692235804

-- Mike

function q = quad3d(f,x1,x2,y1x,y2x,z1xy,z2xy)
q = quadgk(@inner,x1,x2);
    function q = inner(x)
        q = arrayfun( ...
            @(xs,y1s,y2s) ...
            quad2d(@(y,z)f(xs*ones(size(y)),y,z),y1s,y2s, ...
                @(y)z1xy(xs*ones(size(y)),y), ...
                @(y)z2xy(xs*ones(size(y)),y)), ...
                x,y1x(x),y2x(x));
    end
end

[Matt Tearle]

I *knew* I had seen that in the doc before. Why couldn't I find it just now? (Obvious answer: I'm an idiot.)

[Manuel]

Thank you Mike for your answer.

quad2d let's me integrate with function handles over two integrals. It works like a charm. But I need to integrate over three integrals. I tried the following two things:

1. Your function quad3d which gave me an error. I used for testing purposes a simple function.

My Code is (I hope it is readable since it are only two lines):

I3 = @(xi,yi,zi)quad3d(@(r,phi,z)(xi+yi+zi+r+phi+z), 0,1, 0,1, 0,1); %no error

I3(0,0,0) % Error: ??? Subscript indices must either be real positive integers or logicals.

2. Integrate over the inner two integrals and then integrate with quad() over the last one (which has fixed limits). My Code for that is:

I2 = @(xi,yi,zi,r)quad2d(@(z,phi)(xi+yi+zi+r+phi+z), 0,1, 0,1);

I = @(xi,yi,zi)quad(@(r)I2, 0,1);

I(0,0,0) %Gives error: ??? Error using ==> quad at 80

% The integrand function must return an output vector of the same length as the input vector.

Trying to vectorice I2 with

I2v=@(xi,yi,zi,r) arrayfun(I2, xi,yi,zi,r);

leads to the same error.

BTW: I did take a look at your function quad3d and did not find out which value "y" has when you use it in "ones(size(y))" for example.

Manuel

[Mike Hosea]

You said the limits had to be functions, so I didn't code it to accept scalars. The call for your example should be

I3 = @(xi,yi,zi)quad3d(@(r,phi,z)(xi+yi+zi+r+phi+z),0,1,@(r)zeros(size(r)),@(r)ones(size(r)),@(r,phi)zeros(size(r)),@(r,phi)ones(size(r))); I3(0,0,0)

This returns 1.5 on my computer.

Answer 2

[the cyclist]

I think you should be able to do this with dblquad(). To take care of the variable limits, define a rectangular area (with fixed limits) that you know contains the entire range. Then, include an extra multiplicative term in your integrand that becomes zero when you are outside the true limits of your inner integral.

I'm pretty sure that a similar question was asked and answered here before, so maybe do a search on "dblquad", "limit", etc. [Have to admit I did not find it myself, but I did not search thoroughly.]

[Manuel]

Thank you cyclist for the answer.

I tried today to follow your suggestions but did not jet succeed. I also found another post from you: http://www.mathworks.de/matlabcentral/answers/15260-calculating-triple-integral-between-two-surfaces-gained-with-triscatteredinterp-or-griddata

But I find it hard apply on my problem because it is more complex (three integrals).

If you have a more precise adivice feel free telling me :)

Answer 3

[Matt Tearle]

You could do this as an iterated integral using quad (twice). Hooray for Fubini's Theorem! Here's the volume of a hemisphere, calculated as a double integral over the unit circle:

% Define upper and lower limits
limupr = @(x) sqrt(1-x.^2);
limlwr = @(x) -sqrt(1-x.^2);

% Define integrand
f = @(x,y) sqrt(1-x.^2-y.^2);

% Do inner integral w.r.t. y:
% I(x) = integrate_l(x)^u(x) f(x,y) dy
int1 = @(x) quad(@(y) f(x,y),limlwr(x),limupr(x));

% Make sure we can evaluate int1 at a vector of x values
g = @(x) arrayfun(int1,x);
% Integrate w.r.t. x
I = quad(g,-1,1)

[Mike Hosea]

>> quad2d(f,-1,1,limlwr,limupr)

ans =

2.094395102392511

[Manuel]

Hi Matt Tearle and thank you for your suggestion. I feel like it brings me very near to the solution of my Problem. Anyway, if i try to add a third integral it fails.

The error i get is:

??? Error using ==> arrayfun

All of the input arguments must be of the same size and shape.

Previous inputs had size 7 in dimension 2. Input #3 has size 1.

What I am doing is the following:

h1=@(phi)phi.^2;

h2=@(phi)phi.^2-phi;

n=1;

ra=10e-3;

ri=4e-3;

% Limits inner integral (Integral 1)

lim1upr = h1;

lim1lwr = h2;

% Limits Integral 2

lim2upr = n*pi;

lim2lwr = -n*pi;

% Limits Integral 3

lim3upr = ra;

lim3lwr = ri;

% Function

b=@(z,phi,r)1./(h1(phi)-h2(phi)) * (r-(yi.*sin(phi)+xi.*cos(phi))) ./ ((xi-r.*cos(phi)).^2+(yi+r.*sin(phi)).^2+(zi-z(phi)).^2).^(3/2);

% Calculate Integral 1

% int1 = integrate_l1l^l1u...

int1 = @(phi,r) quad(@(z) b(z,phi,r), lim1upr, lim1lwr);

% modify int1 to take vectors (Did not understand that part fully)

int1v = @(phi,r) arrayfun(int1,phi,r);

% Calculate Integral 2

% int2 = integrate_l2l^l2u...

int2 = @(r) quad(@(phi) int1v(phi,r),lim2upr,lim2lwr);

% modify int2 to take vectors (Did not understand that part fully)

int2v = @(r) arrayfun(int2,r);

% Calculate Integral 3

% int3 = integrate_l3l^l3u...

int3 = quad(int2v,lim3upr,lim3lwr); %Does not work!

Do you have any advice?

@Mike Hosea:

In my case I would need quad3d() and I did not find it in Matlab. Is there something like triplequad which takes function-handles?

Regards,

Manu

[Mike Hosea]

Unfortunately, there is not. It's a little tricky to extend QUAD2D, so I'll add the code to my original answer (so that it will be readably formatted).

Answer 4

[Manuel]

@Mike Hosea

Thank you for you reply. I have one integral with functions as limits and two with scalar values. I apologize for being unclear about that. In my original Post there are one integral with scalar limits and one with function limits. I thought it would then be easy to add another scalar integral.

Anyway, I did find an similar solution:

f=@(xi, yi, zi, z, phi, r)(xi+yi+zi+r+phi+z); %Just any function

for testing

I2  = @(xi,yi,zi,r)quad2d(@(z,phi)f(xi,yi,zi,r,phi,z), 0,1, 0,1); % since the use of quad2d (thanks to Mike) one can set a function as boundary (like @(phi)sin(phi))
I2v = @(xi,yi,zi,r) arrayfun(I2, xi,yi,zi,r); %Vectorize
I   = @(xi,yi,zi)quad(@(r)I2v(xi+0.*r,yi+0.*r,zi+0.*r,r), 0,1); %Solve the last integral but add "additional fields" to xi,yi,zi (so they have the same size as r)
Iv  = @(xi,yi,zi) arrayfun(I, xi,yi,zi); %Vectorize

Hm, which Answer do I mark as the right one?

Manu

[Matt Tearle]

So, just to be clear: you're all set now? You got it to work?

If so, it seems like the key point was using quad2d to have function handle limits on the inner integral, so please accept Mike's answer. Thx. (And if not, let us know what's still outstanding.)

[Manuel]

Yes, using quad2d() which accepts function handles in addition with the variable.*r was it. Sometimes I get Errors which are because quad2d can not solve the integral with some specific parameters. But I will work on that.

Thank you all.