积分函数中数值太大怎么办

bedynamic

1. v = 2;
   
2. r = 0;
   
3. for l = 10:50:1000;
   
4.       a = factorial(l-1);
   
5.         R = @(x)x.^(l-1).* exp(-x)./a.* log2 ( 1 + x./( ( l / v ) - 1 ) );   
   
6.         r = integral(R,0,Inf);
   
7.         
   
8.         plot(l,r,'rs');
   
9.         axis([10 1000 0.5 5]);
   
10.         hold on
    
11. end

复制代码

50以后的阶乘就是INF了怎么办？我前面数试了可以积出来，后面数太大了就不行了，有什么办法可以解决？求指点！！！

amnesia

Matlab中的数一律用双精度格式存储，其范围一般而言已经足够了，但也有例外的时候，比如说计算200以上大数的阶乘。
Mathematics擅长符号计算，Matlab擅长数值计算。下面给出Matlab计算阶乘的几个方法：
1： factorial( )函数
>> factorial(4)

ans = 24

4!=4*3*2*1=24

2：借用Gamma函数的性质，当n为正整数时，gamma(n+1)=n！
>> gamma(5)

ans = 24

gamma(N)=(N-1)*(N-2)*...*2*1

3：prod( )函数 prod函数用于求数组元素的乘积
>> A=1:1:100;
>> prod(A)

ans =

9.3326e+157
还可以用它求双阶乘：
>> prod(1:2:7)

ans =

105
4.编写M文件也是可以的...
BUT，以上方法只适合小数阶乘，碰到大数阶乘就inf歇菜了...So,Matlab里面有个Maple引擎用来符号计算的，可以借用符号计算实现大数阶乘：

>> s='4!'
s =4!
>> vpa(s)
ans =24.

>> factorial(365)

ans =

Inf

>> simple(sym('365!'))

simplify:

25104128675558732292929443748812027705165520269876079766872595193901106138220937419666018009000254169376172314360982328660708071123369979853445367910653872383599704355532740937678091491429440864316046925074510134847025546014098005907965541041195496105311886173373435145517193282760847755882291690213539123479186274701519396808504940722607033001246328398800550487427999876690416973437861078185344667966871511049653888130136836199010529180056125844549488648617682915826347564148990984138067809999604687488146734837340699359838791124995957584538873616661533093253551256845056046388738129702951381151861413688922986510005440943943014699244112555755279140760492764253740250410391056421979003289600000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
这里省略了一部分输出..

simple( )这种表现形式很不好，改用极限的形式试一下：
>> syms x;
>> limit('x+365!',x,0)

ans =

25104128675558732292929443748812027705165520269876079766872595193901106138220937419666018009000254169376172314360982328660708071123369979853445367910653872383599704355532740937678091491429440864316046925074510134847025546014098005907965541041195496105311886173373435145517193282760847755882291690213539123479186274701519396808504940722607033001246328398800550487427999876690416973437861078185344667966871511049653888130136836199010529180056125844549488648617682915826347564148990984138067809999604687488146734837340699359838791124995957584538873616661533093253551256845056046388738129702951381151861413688922986510005440943943014699244112555755279140760492764253740250410391056421979003289600000000000000000000000000000000000000000000000000000000000000000000000000000000000000000


>> whos
Name Size Bytes Class Attributes

ans 1x1 1682 sym
x 1x1 126 sym
可以看出，这里ans是一个sym型值，且并不能转换为double值：

>> str2double(ans)

ans =

NaN

还是借助Maple引擎计算，这么表示也可以：
>> s='365!';
>> vpa(s)

ans =

2.5104128675558732292929443748812*10^778

>> str2double(ans)

ans =

NaN

以上三种表示方式核心都是借用Maple引擎计算求解，但是这里ans是一个sym型值，且并不能转换为double值

作者：Eric Huang
链接：http://www.zhihu.com/question/21753495/answer/19241348
来源：知乎
著作权归作者所有。商业转载请联系作者获得授权，非商业转载请注明出处。

bedynamic

amnesia 发表于 2016-3-29 10:53
Matlab中的数一律用双精度格式存储，其范围一般而言已经足够了，但也有例外的时候，比如说计算200以上大数 ...

那上面的没法算了么？？？

Lorraine

bedynamic 发表于 2016-3-30 19:04
那上面的没法算了么？？？

应该是，超出了能力范围