Asked by XIONG Qingxiang
on 12 Jul 2019

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.

Answer by Jon
on 12 Jul 2019

Accepted Answer

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.

XIONG Qingxiang
on 15 Jul 2019

%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
on 15 Jul 2019

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
on 16 Jul 2019

Sign in to comment.

Answer by Matt J
on 12 Jul 2019

XIONG Qingxiang
on 15 Jul 2019

Sign in to comment.

Opportunities for recent engineering grads.

Apply Today
## 0 Comments

Sign in to comment.