Integral over an integral without analytical solution

Hello there,

I'm facing the following challenge:

I have a nested integral without analytical solution wich depends on the same variable as the outer integral. Testing the integral routine with a similar problem (but with analytical solution) yields out that the result is wrong.

The similar problem looks like: Integral dx(Integral dx x). The analytical solution gives (1/6)*x^3.

If we now assume the area of this analytical solution from 0 to 2 we expect 8/6 (round about 1.333). But what MATLAB outputs is 4.

So what happens is that first the function x is numerically integrated from 0 to 2 what gives 2 and then the result is again integrated what gives 4 (so, the integral from 0 to 2 over 2)

Here is the code for the example:

% Call of the integration

    x1 = 0;
    x2 = 2;

    OuterInt(x1,x2)

Function for outer integral:

function [Total] = OuterInt(x1,x2)

OuterFun = @(x) InnerInt(x1,x2);
Total = integral(OuterFun,x1,x2,'ArrayValued',true);

And at last the function for the inner integral:

function [ValInnerInt] = InnerInt(x1,x2) 
                      %x1 lower boarder; x2 upper boarder

    InnerFun = @(x) x;
    ValInnerInt = integral(InnerFun,x1,x2,'ArrayValued',true);

 end

For instance: the inner integral should be evaluated at maximum to the current step in the integration of the outer integral.

Thanks to everybody sharing solutions and ideas with me!

Answer 1

[Walter Roberson]

You should be using integral2() . You can use functions for the lower and upper bounds of the second variable.

[Stefan Brehm]

Thank you Walter. I tried it but it does not work. Meanwhile I've found another way to solve this problem.

I will post it later.

Nevertheless, thanks a lot.

Answer 2

[Stefan Brehm]

As I mentioned in my last post, here is a possibility to solve the problem I descibed in my question:

Calling the outer Integral (just for convenience):

%%This is how solving of convoluted integrals work:

% Call of the outer integral:

x1 = 0;
x2 = 2;

Out = OuterInt(x1,x2)

The outer Integral:

%%Outer Integral receives the upper and lower bounds (x1,x2)
%   where x1 is the lower one. But it also works reverse with the
%   conventional sign convention vor changing the bounds.

% --> OuterIntegral calls InnerInt

function [Tot] = OuterInt(x1,x2)

InnerInt1 = @(x) InnerInt(x);
Total = integral(InnerInt1,x1,x2,'ArrayValued',true,'RelTol',1e-30,'AbsTol',1e-30);
Tot = Total;

end

The inner Integral:

%%Inner Integral receives the upper and lower bounds from the
%   outer integral. So it knows nothing about the initial bounds x1
%   and x2.

% --> OuterIntegral calls InnerInt

function [ValInnerInt] = InnerInt(x) 

Func = @(x) (Fun(x));
ValInnerInt = integral(Func,0,x,'ArrayValued',true,'RelTol',1e-30,'AbsTol',1e-30);

end

The function (one example with analytical solution):

%%Function to be integrated. Since InnerInt and OuterInt expect
%   functions of x as a variable, x should be used as the
%   functional dependency. Otherwise one has to change the
%   corresponding handles.

function [FUN] =Fun(x)

FUN = x;

end

Even though this does not work for at least one function, namely the sin(x), it seems to work for many types of functions.

Best regards Stefan