Numerical Integration matlab Help

1 次查看(过去 30 天)
mutah
mutah 2014-5-2
评论: mutah 2014-5-5
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);

回答(1 个)

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

类别

Help CenterFile 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