Numerical Integration matlab Help

2 views (last 30 days)
mutah
mutah on 2 May 2014
Commented: mutah on 5 May 2014
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);

Answers (1)

Roger Stafford
Roger Stafford on 2 May 2014
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'.

Community Treasure Hunt

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

Start Hunting!