Triple Numerical Integration with function boundaries

Hello,

I am having an issue quite similar to another question that has been previously addressed:

http://www.mathworks.com/matlabcentral/answers/17517-numerical-integration-with-functions-as-borders

However, what I need differs slightly.

Say instead of phi, we have phi_1 and phi_2, and say that in the integral from the original question, x ranges between h_2(phi_1,phi_2) and h_1(phi_1,phi_2), and then replace the single integration with respect to phi with a double integral over phi_1 and phi_2 where the ranges over phi_1 and phi_2 are expressed by scalars, how would we compute such messy triple integral?

Any guidance would be much appreciated.

Answer 1

[Andrew Newell]

Here is a sketch of the approach:

First, define your function to be integrated, e.g.,

f0 = @(x,phi1,phi2) 1./(h1(phi1,phi2)-h2(phi1,phi2))./(b-(xi-x).^2).^1.5;

(xi, h1, h2 to be specified by you). Define the inner integral:

f1 = @(phi1,phi2) quadl(f0,h1(phi1,phi2),h2(phi1,phi2));

And then integrate over the other two variables:

f2 = dblquad(f1,phi1_lb,phi1_ub,phi2_lb,phi2_ub);

[James]

Thank you Andrew for your insight. I've tried it out but I'm getting an error that one of my variables are undefined.

Here's what I tried:

f0 = @(g_1,g_0,g_2)(g_0.^0.1).*exp(-0.1.*g_0)...

.*(g_1.^0.2).*exp(-0.2.*g_1)...

.*(g_2.^0.3).*exp(-0.3.*g_2);

f1 = @(g_0,g_2)quadl(f0,g_0,(1+g_0).*(1+g_2)-1);

f2 = dblquad(f1,0,1000,0,1000)

Basically I want to find the triple integral of the function defined by f0, over the range g_1 = [g_0 (1+g_0)*(1+g_2)-1], followed by integrating g_0 and g_2 over [0,inf].

The return states that g_0 is undefined.

[Andrew Newell]

I didn't test my "sketch" very carefully because I knew it wasn't the function you really wanted to calculate. I think this works:

f0 = @(g_1,g_0,g_2)(g_0.^0.1).*exp(-0.1.*g_0)...

.*(g_1.^0.2).*exp(-0.2.*g_1)...

.*(g_2.^0.3).*exp(-0.3.*g_2);

f1 = @(g_0,g_2)quadv(@(x)f0(x,g_0,g_2),g_0,(1+g_0).*(1+g_2)-1);

f2 = dblquad(f1,0,1000,0,1000);

[James]

Thank you so much. Would you mind explaining why we needed the extra @(x) in front? I've been trying to read some manuals regarding how this works. My main confusion is that x does not actually appear in any of the expressions, yet we declare it as a variable within the quadv function. Somehow I assume this is also related to switching from quadl to quadv.

Sorry, semi-matlab newbie.

[Andrew Newell]

It's because you need a function of x only to integrate but you can't just do something like

f_onevar = @(x) f0(x,g_0,g_2);

f1 = (g_0,g_2)quadv(f_onevar,g_0,(1+g_0).*(1+g_2)-1);

because you don't know g_0 and g_2 in advance. So you have to define f_onevar *inside* f1. This would probably not look so surprising if you formulated the problem in terms of functions in separate files (but to do that you'd need nested functions).