LISP的积分问题

ynhh

cad890 发表于 2019-7-4 18:18
应该是图片化简过程出错了，绿色线框处积分化简使用条件为a>0，即F>0；
而F不能保证大于0（如φ=0时，B= ...

之前说不是这个图
也不知我说的图对不对
如不对就算了
太麻烦您真是不好意思
那就算了吧
您也不用再费心了
谢谢您的多次费心指导
谢谢

highflybird

cad890 发表于 2019-7-4 18:18
应该是图片化简过程出错了，绿色线框处积分化简使用条件为a>0，即F>0；
而F不能保证大于0（如φ=0时，B= ...

CAD890很有见识，的确是要考虑小于0 的情况，如果考虑小于0 的情况，则是正确的，用LISP得到的结果也正确了。

highflybird

已经用龙贝格积分法解决，需要lisp的移步到http://bbs.mjtd.com/forum.php?mo ... mp;extra=#pid844670查看。2#楼贴出了源代码。

ynhh

highflybird 发表于 2019-7-5 20:03
已经用龙贝格积分法解决，需要lisp的移步到http://bbs.mjtd.com/forum.php?mod=viewthread&tid=179815&page ...

谢谢高大师
已选择您为最佳答案
应是悬赏已兑现
我好好学习一下
谢谢您的耐心指导

highflybird

在此，我还是 贴出源代码，可以看出，高斯-切比雪夫法在求此类型积分上有很高的效率。

1. ;;;=============================================================
   
2. ;;; 高斯-切比雪夫积分                                          
   
3. ;;; 此方法用于针对 sqrt(1-x^2)*f(x)型的积分有很高的效率。      
   
4. ;;; 输入: foo 函数名,arg 除自变量x外的参数列表, n迭代次数。     
   
5. ;;; 输出: 此类积分的数值.                                       
   
6. ;;; 说明: 此积分法效率较高，n取值一般10次左右就达到有效浮点精度.
   
7. ;;; http://mathworld.wolfram.com/Chebyshev-GaussQuadrature.html
   
8. ;;;=============================================================
   
9. (defun Math:Gauss-Chebyshev (foo args n)
   
10.   (setq wi (/ pi n))
    
11.   (setq sx 0)
    
12.   (setq i 1)
    
13.   (repeat n
    
14.     (setq xi (cos (/ (* pi (- i 0.5)) n)))
    
15.     (setq fi (apply foo (cons xi args)))
    
16.     (setq sx (+ sx fi))
    
17.     (setq i  (1+ i))
    
18.   )
    
19.   (* wi sx)
    
20. )
    
21. 
    
22. ;;;=============================================================
    
23. ;;; 用高斯-切比雪夫积分法求椭球带面积                           
    
24. ;;; 输入: a,b,c，椭球的x,y,z半轴长。h是其高(Z方向)与c的比值(-1<=h<=1)
    
25. ;;; 输出: 所求椭球带的表面积.                                   
    
26. ;;;=============================================================
    
27. (defun ELL:GetSurfaceArea3 (a b c h / func v)
    
28.   (defun Func (x a b c h / F G K R Y Z)
    
29.     (setq k (* a a))
    
30.     (setq g (- 1 (* (- 1 (/ k b b)) x x)))
    
31.     (setq f (1- (/ k c c g)))
    
32.     (setq y (1+ (* f h h)))
    
33.     (if        (zerop f)
    
34.       (* (sqrt g) (1+ (sqrt y)))
    
35.       (if (< f 0)
    
36.         (setq z        (* h (sqrt (- f)))
    
37.               r        (* (sqrt g) (+ (sqrt y) (/ (asin z) z)))
    
38.         )
    
39.         (setq z        (* h (sqrt f))
    
40.               r        (* (sqrt g) (+ (sqrt y) (/ (asinh z) z)))
    
41.         )
    
42.       )
    
43.     )
    
44.   )
    
45.   (setq v (/ (+ a b c) 3.0))                                ;设置与三个半轴长的平均值，防止误差增大
    
46.   (setq a (/ a v))                                                ;X 半轴长与平均值的比值
    
47.   (setq b (/ b v))                                                ;Y 半轴长与平均值的比值
    
48.   (setq c (/ c v))                                                ;Z 半轴长与平均值的比值
    
49.   (* b c h v v (Math:Gauss-Chebyshev 'func (list a b c h) 20))  ;在此例程中设置20可以满足精度
    
50. )
    

;;;=============================================================<br /> ;;; 高斯-切比雪夫积分                                           <br /> ;;; 此方法用于针对 sqrt(1-x^2)*f(x)型的积分有很高的效率。       <br /> ;;; 输入: foo 函数名,arg 除自变量x外的参数列表, n迭代次数。     <br /> ;;; 输出: 此类积分的数值.                                       <br /> ;;; 说明: 此积分法效率较高，n取值一般10次左右就达到有效浮点精度.<br /> ;;; <a href="http://mathworld.wolfram.com/Chebyshev-GaussQuadrature.html" target="_blank">http://mathworld.wolfram.com/Chebyshev-GaussQuadrature.html</a> <br /> ;;;=============================================================<br /> (defun Math:Gauss-Chebyshev (foo args n)<br />   (setq wi (/ pi n))<br />   (setq sx 0)<br />   (setq i 1)<br />   (repeat n<br />     (setq xi (cos (/ (* pi (- i 0.5)) n)))<br />     (setq fi (apply foo (cons xi args)))<br />     (setq sx (+ sx fi))<br />     (setq i  (1+ i))<br />   )<br />   (* wi sx)<br /> )<br /> <br /> ;;;=============================================================<br /> ;;; 用高斯-切比雪夫积分法求椭球带面积                           <br /> ;;; 输入: a,b,c，椭球的x,y,z半轴长。h是其高(Z方向)与c的比值(-1<=h<=1)<br /> ;;; 输出: 所求椭球带的表面积.                                   <br /> ;;;=============================================================<br /> (defun ELL:GetSurfaceArea3 (a b c h / func v)<br />   (defun Func (x a b c h / F G K R Y Z)<br />     (setq k (* a a))<br />     (setq g (- 1 (* (- 1 (/ k b b)) x x)))<br />     (setq f (1- (/ k c c g)))<br />     (setq y (1+ (* f h h)))<br />     (if        (zerop f)<br />       (* (sqrt g) (1+ (sqrt y)))<br />       (if (< f 0)<br />         (setq z        (* h (sqrt (- f)))<br />               r        (* (sqrt g) (+ (sqrt y) (/ (asin z) z)))<br />         )<br />         (setq z        (* h (sqrt f))<br />               r        (* (sqrt g) (+ (sqrt y) (/ (asinh z) z)))<br />         )<br />       )<br />     )<br />   )<br />   (setq v (/ (+ a b c) 3.0))                                ;设置与三个半轴长的平均值，防止误差增大<br />   (setq a (/ a v))                                                ;X 半轴长与平均值的比值<br />   (setq b (/ b v))                                                ;Y 半轴长与平均值的比值<br />   (setq c (/ c v))                                                ;Z 半轴长与平均值的比值<br />   (* b c h v v (Math:Gauss-Chebyshev 'func (list a b c h) 20))  ;在此例程中设置20可以满足精度<br /> )<br />

highflybird

highflybird 发表于 2019-7-8 23:06
在此，我还是 贴出源代码，可以看出，高斯-切比雪夫法在求此类型积分上有很高的效率。

此类求法与龙贝格求得的一样的结果，不过迭代次数更少。

ynhh

highflybird 发表于 2019-7-8 23:07
此类求法与龙贝格求得的一样的结果，不过迭代次数更少。

感谢高飞大师
您的精神是全国科技工作者学习的榜样