Improper integral at 0

lambda = 0.07;
mu_1 = 0.05;
mu_2 = 0.12;
W0 = 20000;
R1 = 3000;
R2 = 5000;
r = 0.05;

f = @(x)  ((R1-r*W0)/x)^((lambda + mu_1) / r) * (R2 - R1 + x)^(mu_2/r-1); % Example function
a = 0;         % Lower limit of integration
b = R1;         % Upper limit of integration
tol = 1e-6;    % Tolerance (optional)
result = quadadapt(f, a, b, tol);

Out of memory. The likely cause is an infinite recursion within the program.

Error in solution>quadstep (line 43)
qa = quadstep(f, a, c, tol, fa, fd, fc, varargin{:});

disp(result);
function q = quadadapt(f, a, b, tol, varargin)
% Checks if the tolerance is provided, otherwise sets a default tolerance
if nargin < 4 || isempty(tol)
    tol = 1.e-6;
end
% Initial function evaluations at the endpoints and the midpoint
c = (a + b) / 2;
fa = feval(f, a, varargin{:});
fc = feval(f, c, varargin{:});
fb = feval(f, b, varargin{:});
% Call the recursive adaptive quadrature function
q = quadstep(f, a, b, tol, fa, fc, fb, varargin{:});
end
function q = quadstep(f, a, b, tol, fa, fc, fb, varargin)
% Calculate the midpoint and the step size
h = b - a;
c = (a + b) / 2;
% Function evaluations at the quarter-points
fd = feval(f, (a + c) / 2, varargin{:});
fe = feval(f, (c + b) / 2, varargin{:});
% Simple quadrature and more refined quadrature
q1 = h / 6 * (fa + 4 * fc + fb);
q2 = h / 12 * (fa + 4 * fd + 2 * fc + 4 * fe + fb);
% Check if the refinement is within the tolerance
if abs(q2 - q1) <= tol
    q = q2 + (q2 - q1) / 15;
else
    % Recursive calls for each half of the interval
    qa = quadstep(f, a, c, tol, fa, fd, fc, varargin{:});
    qb = quadstep(f, c, b, tol, fc, fe, fb, varargin{:});
    q = qa + qb;
end
end

I am trying to approximate improper integral ((R1-r*W0)/x)^((lambda + mu_1) / r) * (R2 - R1 + x)^(mu_2/r-1) from o to R1.

Can you suggest method or or code example. It seems quadadapt doesn't work for improper integral

Answer 1

[Torsten]

Your integrand behaves like x^(-12/5) at x=0. That means that the value of the integral is + Inf.

lambda = 0.07;
mu_1 = 0.05;
mu_2 = 0.12;
W0 = 20000;
R1 = 3000;
R2 = 5000;
r = 0.05;
syms x
f = ((R1-r*W0)/x)^((lambda + mu_1) / r) * (R2 - R1 + x)^(mu_2/r-1) % Example function

f = 

[Raushan]

Can you suggest method for solving improper integrals? This method doesn't work EVEN for 1/sqrt(x) which is convergent. My integral should be convergent in the finite interval.

[Torsten]

f = @(x) 1./sqrt(x)

f = function_handle with value:
    @(x)1./sqrt(x)

integral(f,0,1)

ans = 2.0000

syms x
f = 1/sqrt(x);
int(f,0,1)

ans = 

2

[Raushan]

It can solve only simple integrals, even if I change 0 to 0.01, it gives not correct answer. Probably there are other methods?

[Raushan]

I need the method which can numerically approximate the improper integral

[Paul]

@Raushan

What is the basis for the assertion that " My integral should be convergent in the finite interval. "?

Here is the result from integral to integrate from 0.01 to R1 using default settings:

lambda = 0.07;
mu_1 = 0.05;
mu_2 = 0.12;
W0 = 20000;
R1 = 3000;
R2 = 5000;
r = 0.05;
syms x
f(x) =  ((R1-r*W0)./x).^((lambda + mu_1) ./ r) .* (R2 - R1 + x).^(mu_2/r-1); % Example function
integral(matlabFunction(f(x)),0.01,R1)

ans = 1.5769e+15

Here is the result from vpaintegral using default settings:

vpaintegral(f(x),x,0.01,R1)

ans = 

1.57687e+15

Are these not close enough to the correct answer? What is the correct answer?

Assuming vpaintegral is returning a correct result, it certainly looks like the improper integral does not converge

a = logspace(-6,-2,5);
for ii = 1:numel(a)
y(ii) = vpaintegral(f(x),x,a(ii),R1);
end
figure
loglog(a,double(y),'-o'),grid

[John D'Errico]

@Raushan

Why are you using your own code? For example, integral will survive when it is possible. NEVER write your own code to solve a problem when well written code exists, written by professionals. For example...

f = @(x) 1./sqrt(x);
integral(f,0,1)

ans = 2.0000

As you can see, integral has no problem, as long as the integral exists.

However, your problem does not have a solution, no matter what you do. And reducing the interval on an unbounded integral is not meaningful. It will generate a non-zero result that is complete garbage.

[Raushan]

@Paul @John D'ErricoThank you very much! I see your point. It is probability and when I simulated I got finite value. Probably I will re-consider my integral again. I also considered different packages in different languages and got large numbers and different ones. Therefore, I thought I wrote code incorrectly, and I assumed written codes(packages) doesn't approximate integrals with singularities.