How to integrate an unbounded function?

2019 7 12
2 个回答
回答已采纳
更新时间：2019 7 16
6 次查看（30 天）

0 个投票

Hello every one,
since the normal numerical integration cannnot work to an unbounded function, is there any other methods for the integration of unbounded function?.
Thank you.

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

 采纳的回答

Jon
Jon 2019-7-12

1 个投票

According to the MATLAB documentation https://www.mathworks.com/help/matlab/ref/integral.html MATLAB's integral function can handle singlularities on the boundary (upper or lower limit) of the domain to be integrated. If the singularity is not on the boundary, you could always break up the integral into the sum of two integrals, one with the singularity at the upper bound, the other with the singularity at the lower bound. If you have a finite number of singularities in your domain, you can extend this idea, and break up the overall integration into pieces each with singularities at their endpoints.
I am assuming the function you are integrating is one dimensional. MATLAB also has routines, integral2 and integral3 for two and three dimensional cases. I think these can also handle singularities at the boundaries, but you can check the documentation for details.

7 个评论
显示 5更早的评论 隐藏 5更早的评论

XIONG Qingxiang
XIONG Qingxiang 2019-7-15
Thank you very much. I'm integrating on one-dimensional function. And I've separated the integral interval at the two singularity points, but it still doesn't work. It seems that my function is an exception.
Bruno Luong
Bruno Luong 2019-7-15
Dedends on the singularity, if it converges then it give the right anwser
>> quad(@(t) 1./sqrt(t),0,1)
ans =
2.0000
It if doesn't then it likely that your integral is not defined (diverge).
XIONG Qingxiang
XIONG Qingxiang 2019-7-15
Yes, that's might be the answer for improper integration. When the limitation exists, the integration convergence, otherwise it should be divergence.
Actually, I have a loop in my code, when the integrand is non-monotonic, the syntax can calculate only one time, but for the second time, it fails. So I wonder if it can work for one time, why it doesn't work for the second loop?
Bruno Luong
Bruno Luong 2019-7-15
You probably need to tell us more about the type of singularies you are dealing with and what do you change between loop iteration.
XIONG Qingxiang
XIONG Qingxiang 2019-7-15
编辑：XIONG Qingxiang 2019-7-15
%Here is the code,and I made some comments for you understanding. Thank you very much.
syms r
phi0=0.22;
rmin = 1.4824*10^-9;
rmax = 1.6720*10^-4;
sigma = 0.31;
mu = -17.1043;
z = 10^6;
r0 = linspace(rmin,rmax,z);
P0 = lognpdf(r0,mu,sigma);
rc =3e-8;
beta = 3.52e-9;
epsmin = -2*beta/3;
nloop =3;
r0ave=exp(mu+(sigma^2)/2);
mystruct(nloop) = struct('Matphid',[],'MatPHID',[],'Mateps',[]);
k =1;
N = zeros(1,nloop);
feval(symengine,'_assign','upper1',rc);
feval(symengine,'_assign','lower1',rc-beta);
feval(symengine,'_assign','upper5',rc+beta);
feval(symengine,'_assign','lower5',rc);
feval(symengine,'_assign','upper2',inf);
feval(symengine,'_assign','lower2',rc+beta);
feval(symengine,'_assign','upper3',rc-beta);
feval(symengine,'_assign','lower3',0);
feval(symengine,'_assign','upper4',rc);
feval(symengine,'_assign','lower4',0);
H = (3/4)*((2*r+beta-2*rc)/beta - (1/3)*((2*r+beta-2*rc)/beta)^3) + 1/2;
dH = (3/2)*(1/beta - ((2*r+beta-2*rc)^2)/(beta^3));
M = (-3/4)*((2*r-beta-2*rc)/beta - (1/3)*((2*r-beta-2*rc)/beta)^3) + 1/2;
dM = (-3/2)*(1/beta - ((2*r-beta-2*rc)^2)/(beta^3));
P = (1/(sigma*sqrt(2*pi)))*(1/r)*exp(-((log(r) - mu)^2)/(2*sigma^2));
while (k <= nloop)
n = k;
N(k) = n;
eps = zeros(1,n+1);
F = zeros(1,n+1);
G = zeros(1,n+1);
phid = zeros(1,n+1);
S3 = zeros(1,n+1);
S4 = zeros(1,n+1);
phid(1) = phi0;
PHID = zeros(1,n);
rn = zeros(n+1,z);
rn(1,:) = r0;
Pn = zeros(n+1,z);
Pn(1,:) = P0;
PHid = linspace(0.219,217,n);
for i = 1:n
if i == 1
F =1/(1+eps(i)*dH);
G =1/(1+eps(i)*dM);
else
F = F/(1+eps(i)*dH);% (1+eps(i)*dH) is the derivative, it could be zero at two staionary points
G = G/(1+eps(i)*dM);% (1+eps(i)*dH) is the derivative.
end
f1 = P;
f2 = r*F*P;
f3 = r*G*P;
f4 = r*P;
f5 = H*F*P;
f6 = M*G*P;
Pd1 = P;
Pd2 = F*P;
Pd3 = M*P;
% J1 and J2 are the integrations
J1 = eval(feval(symengine,'numeric::int',f4,'r=lower3..upper3')) + ...
eval(feval(symengine,'numeric::int',f2,'r=lower1..upper1')) + ...
eval(feval(symengine,'numeric::int',f3,'r=lower5..upper5')) + ...
eval(feval(symengine,'numeric::int',f4,'r=lower2..upper2'));
J2 = eval(feval(symengine,'numeric::int',f5,'r=lower1..upper1')) + ...
eval(feval(symengine,'numeric::int',f6,'r=lower5..upper5'));
PHID(i) = phi0*((epsmin*J2 + J1)/r0ave)^2;
phid(i+1) = PHid(i);
eps(i+1) =(sqrt(phid(i+1)/phi0)-(1/r0ave)*J1)/((1/r0ave)*J2);
% S3 and S4 are the function value of two stationary points of function (r+eps(i)*H)
S3(i+1)=(-sqrt(3*eps(i+1)*(2*beta+3*eps(i+1)))*(2*beta+3*eps(i+1))-9*...
eps(i+1)*(beta-2*rc-eps(i+1)))/(18*eps(i+1));
S4(i+1)=(sqrt(3*eps(i+1)*(2*beta+3*eps(i+1)))*(2*beta+3*eps(i+1))-9*...
eps(i+1)*(beta-2*rc-eps(i+1)))/(18*eps(i+1));
end
mystruct(k).Matphid = phid;
mystruct(k).Mateps = eps;
k = k + 1;
end
disp(eps)
Bruno Luong
Bruno Luong 2019-7-15
Not sure with those complicated formula (and a lot of them is not important for the discussion) but if I understand
P(r) = (1/(sigma*sqrt(2*pi)))*(1/r)*exp(-((log(r) - mu)^2)/(2*sigma^2))
which is equivalent to O(1/r) * O(1/r^(1/(2*sigma^2))) for r ~0.
The first term already induce non-defined integral. So there is no point trying to compute such thing.
Back to the blackboard.
PS: it looks like you are trying to integrate some EM potential layers.
XIONG Qingxiang
XIONG Qingxiang 2019-7-16
In fact, it's a log-normal distribution and this code is for the propose of solving the parameter epsilon. What I'm working on is porous media. Anyway, thank you for your warming help and I will vote your response as the answer.

请先登录，再进行评论。

更多回答（1 个）

Matt J
Matt J 2019-7-12

0 个投票

You cantry use the symbolic toolbox, or you could just take the numerical integral over a sufficiently large finite interval. A numerical integration is an approximation anyway, so who cares about the extra epsilon of error that you get when you truncate the interval to something finite.

1 个评论
显示 -1更早的评论 隐藏 -1更早的评论

XIONG Qingxiang
XIONG Qingxiang 2019-7-15
Thank you for your idea. But could you please give me more details that how to use symbolic toolbox to solve this problem?

请先登录，再进行评论。

类别

在 帮助中心 和 File Exchange 中查找有关 Numerical Integration and Differentiation 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by