Improper integral at 0

2 次查看（过去 30 天）
Raushan
Raushan 2024-2-1
编辑： Raushan 2024-2-2
Ran in:
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

请先登录，再回答此问题。

回答（1 个）

Torsten
Torsten 2024-2-1
编辑：Torsten 2024-2-1
Ran in:
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 = 
  7 个评论
显示 5更早的评论
John D'Errico
John D'Errico 2024-2-2
Ran in:
@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
Raushan 2024-2-2
编辑：Raushan 2024-2-2
@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.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Calculus 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by