Numerical Integration matlab Help

Hi dears

this is my code please help me and check it

close ; clear ; clc; g=4.6;

k=0.0000000000001;

m=1;n=100;

rou=3.5/sqrt(2);

 p=0.00001:0.01:0.81;

c=(rou*sqrt(pi*m)/sqrt(3*g*(1+k)))*(1/(m*exp(m*k)))*(gamma(m+0.5)/gamma(m))*kummer(m+0.5,m,m*k);

f =@(r)((1-erfc(n*p*c/r)).*((2*m.*(1+k).^((m+1)./2)./k.^((m-1)./2)./exp(m.*k)).* r.^(m).*exp(-m.*(1+k).*r.^2).*besseli(m-1, 2*m.*sqrt(k.*(k+1)).*r)));

F=integral(f,0,inf);

[José-Luis]

kummer is not a Matlab function.

[mutah]

kummer defined as in this link

http://www.mathworks.com/matlabcentral/fileexchange/29766-confluent-hypergeometric-function/content/kummer.m

Answer 1

[Roger Stafford]

To compute this integration numerically you will have to do some kind of simplification in expressing the integrand and in calculating the value of the variable 'c'.

In the expression for 'c', at one point you divide by exp(m*k) where m*k = 1e-13. Using 'double' numbers this would certainly underflow to zero, so you would be dividing by zero at this point, which produces an infinity. The third argument of 'kummer' is also 1e-13. I do not know what this produces but I would suspect it also underflows to zero. If so, the net effect would be infinity times zero for which matlab would produce a NaN.

That same kind of problem is present in calculating the values of the anonymous function 'f'. There are tiny 'k' values within 'exp' and the bessel function in its definition, which will undoubtedly create the same deleterious effect.

As a further difficulty, in the anonymous function 'f' you have included the vector 'p' which has 81 elements. That is not the proper way to define a function which is to serve as an integrand. The integrand should produce a scalar value for 'f' for each individual value of 'r'.

[mutah]

Hi dear

the k=0 and m=1 so exp(k*m)=exp(0*1)=exp(0)=1

also kummer is define by this

http://www.mathworks.com/matlabcentral/fileexchange/29766-confluent-hypergeometric-function/content/kummer.m