求一复杂函数的数值积分

hankeyi

请用数学软件计算附件算式的数值积分，谢谢指点！


I0是零阶第一类修正贝塞尔函数，
m=8
0.03<t<1，当t =0.625145时，D(t)=0.391687就对了！

hankeyi

大家千万不要小看这个问题呀！

本人计算多次都提示溢出，不知何故。把积分下限调到0.00002左右能算，结果与文献不一样。当积分下限为0.01时与文献结果接近，真恼人啊！请高手指点。

wanglu

Forcal求解（高斯方法）：

1. f(s:ms,k,l,ss:m)=
   
2. {
   
3.   ms=m*m*s*s, ss=0,
   
4.   k=1,(k<=m).while{
   
5.     l=1,(l<=m).while{
   
6.       ss=ss+k*l*exp[-(k*k+l*l)/(4*ms)]*XSLSF::bessel3[0,k*l/(2*ms)],
   
7.       l++
   
8.     },
   
9.     k++
   
10.   },
    
11.   ss/(s*s)
    
12. };
    
13. yx(j,x1,y0,y1::tao)=
    
14. {
    
15.     which{
    
16.         j==0:[y0=0,y1=tao]
    
17.     }
    
18. };
    
19. D(_tao,_m::m,tao)= tao=_tao, m=_m, XSLSF::igaus[HFor("f"),HFor("yx"),1]/[m*(m+1)]^2;  //只能取1
    
20. D(0.625145,8);

复制代码

结果：

0.3971074065976576

这个问题的确很怪，徐士良算法中似乎只有高斯方法可解，且只能取一个区间，精度比较低，但有结果。若取大于1的区间数，就求不出结果。区间数越大越精确。

wanglu

是不是随m的增加，函数值（t =0.625145时）收敛到一个常数，该常数大约是：0.4167066699301382

hankeyi

应该是的，这是有物理意义的一个问题，但是我也没搞清约束条件和实际参数间的关系，所以还不能给您提供线索。您确实是帮了大忙了，非常感谢！！！不知这个问题用mathematica算需要加什么条件？

wanglu

楼主这个问题，用双精度实数运算似乎解不了。
bessel3[0,600]＝6.14630437520941e+258
exp[-600]＝2.650396553004311e-261
数再大就超出双精度实数的表示范围了。

当s很小时的值不能忽略，否则误差太大。故用双精度实数运算似乎解不了。

wanglu

用下面的程序只能近似计算，比实际值偏小：

1. f(s:ms,k,l,ss:m)=
   
2. {
   
3.   if[s<2.65413062591e-2, s=2.65413062591e-2],
   
4.   ms=m*m*s*s, ss=0,
   
5.   k=1,(k<=m).while{
   
6.     l=1,(l<=m).while{
   
7.       ss=ss+k*l*exp[-(k*k+l*l)/(4*ms)]*XSLSF::bessel3[0,k*l/(2*ms)],
   
8.       l++
   
9.     },
   
10.     k++
    
11.   },
    
12.   ss/(s*s)
    
13. };
    
14. yx(j,x1,y0,y1::tao)=
    
15. {
    
16.     which{
    
17.         j==0:[y0=0,y1=tao]
    
18.     }
    
19. };
    
20. D(_tao,_m::m,tao)= tao=_tao, m=_m, XSLSF::igaus[HFor("f"),HFor("yx"),100]/[m*(m+1)]^2;
    
21. D(0.625145,8);

复制代码

结果：
0.3911937913874259

hankeyi

那怎么办呢！我要一次计算200个不同Tau值的积分！

hankeyi

也有用laplace变换法算的，搞不懂。

hankeyi

谢谢您的帮助！现在的相对误差在1.2%，看看精度能否再提高一点？！

hankeyi

搞错了，误差0.12%，勉强可用了，再次感谢Wanglu！

pasuka

楼主这个问题，用双精度实数运算似乎解不了。
bessel3[0,600]＝6.14630437520941e+258
exp[-600]＝2.650396553004311e-261
数再大就超出双精度实数的表示范围了。

当s很小时的值不能忽略，否则误差太大。故用双 ...
wanglu 发表于 2010-1-18 20:15

1、用Real*16，看看有没有超出范围；
2、建议lz找本计算方法的书先看看，一般开头都会将到大数乘小数如何避免的问题，然后再修改算法

wanglu

1、用Real*16，看看有没有超出范围；
2、建议lz找本计算方法的书先看看，一般开头都会将到大数乘小数如何避免的问题，然后再修改算法
pasuka 发表于 2010-1-19 15:13

指数函数exp拆开算很容易，但贝塞尔函数不知道如何拆开算，如果有类似如下公式：
bessel3[0,x]＝f(bessel3[0,x/2])   //f是一个函数
问题就好解决了。

hankeyi

查到
对于充分大的x有   
  B(n，x)=sqrt(2/(Pi*x))*cos(x-Pi/4-n*Pi/2)

不知可用否？

hankeyi

对大自变量x>>ABS(n^2-1/4),修正贝塞尔函数的渐近形式为：
In(x)=1/(2sqrt(Pix)) e^x

hankeyi

wanglu

好几个式子，应该用哪一个？
B(n，x)=sqrt(2/(Pi*x))*cos(x-Pi/4-n*Pi/2)  //误差大
In(x)=1/(2sqrt(Pix)) e^x 是不是：B(n，x)=1/(2sqrt(Pix)) e^x

第二个式子误差小，但带入程序计算结果误差较大：0.4967803013415087

wanglu

以下程序验证了第二个式子：
ff(s,m,k,l:ms)= ms=m*m*s*s,k*l*exp[-(k*k+l*l)/(4*ms)]*XSLSF::bessel3[0,k*l/(2*ms)];
gg(s,m,k,l:ms,pi)= ms=m*m*s*s,pi=3.1415926,sqrt[k*l*ms/2/pi]*exp[-((k-l)^2)/(4*ms)];
ff[0.006,8,1,2];
gg[0.006,8,1,2];

结果：
2.87964667680323e-049
2.035631132829438e-049

wanglu

Forcal代码如下，选一个合适的式子，自己再调试一下吧：

1. 
   
2. f(s:ms,k,l,ss,pi:m)=
   
3. {
   
4.   ms=m*m*s*s, ss=0, pi=3.1415926,
   
5.   k=1,(k<=m).while{
   
6.     l=1,(l<=m).while{
   
7.       ss=ss+which{k*l/(2*ms)<700, k*l*exp[-(k*k+l*l)/(4*ms)]*XSLSF::bessel3[0,k*l/(2*ms)], sqrt[k*l*ms/2/pi]*exp[-((k-l)^2)/(4*ms)]},
   
8.       l++
   
9.     },
   
10.     k++
    
11.   },
    
12.   ss/(s*s)
    
13. };
    
14. yx(j,x1,y0,y1::tao)=
    
15. {
    
16.     which{
    
17.         j==0:[y0=0,y1=tao]
    
18.     }
    
19. };
    
20. D(_tao,_m::m,tao)= tao=_tao, m=_m, XSLSF::igaus[HFor("f"),HFor("yx"),100]/[m*(m+1)]^2;
    
21. D(0.625145,8);
    

复制代码

hankeyi

辛苦了，多谢！！！