LISP的积分问题

ynhh · 发表于 2019-7-1 10:11

本帖最后由 ynhh 于 2019-7-7 08:06 编辑

椭球体球缺的表面积计算
用积分公式用LISP程序表达并计算出来
这是一个网上的
有结果和公式

highflybird · 发表于 2019-7-8 23:06

本帖最后由 highflybird 于 2019-7-8 23:40 编辑

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

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

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

highflybird · 发表于 2019-7-2 12:45

ynhh 发表于 2019-7-2 11:31
高飞大师您好
能不能就这个实例
写个LISP的实用例子

;;;=============================================================
;;; 说明: 此函数为纯数学计算，可不用图元作参，采用龙贝塔积分法。
;;; 适用于某些特殊情况。一般可用vlax-curve函数求。见样例。
;;; 功能: 获取椭圆弧的长度。
;;; 输入: 椭圆的长半轴，短半轴，参数1，参数2和精度
;;; 输出: 椭圆弧的长度
;;;=============================================================
(defun ELL:Length (la lb p1 p2 eps / Func ratio 0.5Pi)
(defun Func (x / cx ee)
(setq cx (cos x))
(setq ee (* (1+ ratio) (1- ratio)))
(sqrt (1+ (* ee cx cx)))
)
(if (<= lb la)
(progn
(setq ratio (/ lb (float la)))
(* la (Math:Romberg p1 p2 eps))
)
(progn
(setq ratio (/ la (float lb)))
(setq 0.5Pi (* pi 0.5))
(* lb (Math:Romberg (- p1 0.5pi) (- p2 0.5pi) eps))
)
)
)

譬如这个求椭圆弧长的例子，就是采用了龙贝格积分。
只要你把积分的函数代入，然后代入下限上限值，就可以得出结果。

cad890 · 发表于 2019-7-3 19:42

本帖最后由 cad890 于 2019-7-3 20:03 编辑

ynhh 发表于 2019-7-3 14:14
我这怎么研究都不对
(setq a 2.0 b 3.0 c 4.0 hc 3.0);已知
(setq hb (- c hc))

你这公式推论好像有错误，按理说F(t)应该不会小于0的，但是计算过程发现F(t)会存在小于0的情况；
如n=100时，xi≈cos(1/200 * pi)≈1,则B(t)≈a^2/b^2, F(t)≈b^2/c^2-1，如果b小于c，会计算出F(t)<0。由于存在F(t)开方，所以下方函数是无法计算下去的。请核实一下公式推论F(t)。
但如果a=b=c，hc=0时，可以计算。即半球面积

注：我中午列的公式是根据图片符号，区分大小写，lisp变量不区分大小写，所以你有很多变量重复了。

(defun c:tt()
(setq a 2.0 b 3.0 c 4.0 hc 3.0);已知
(setq hb (- c hc))
(setq n 100)
(setq i 1 h (/ hb c) s 0)
(repeat n
(setq wi (/ pi n)
ai (* pi (/ (- (* 2.0 i) 1) (* 2.0 n)))
xi (cos ai)
;;; B (…)
;;; F (…)
)
(setq Bt (- 1.0 (* (- 1.0 (/ (* a a) (* b b))) xi xi) ))
(setq Ft (- (/ (* a a) (* c c) Bt) 1))
(if (= Ft 0) (setq f (* 2 h (/ a c)))
;;;（为了避免程序又臭又长，可以增加中间计算）
;;;(setq f (1- (* (* (/ (* A A) (* C C))) (/ 1.0 B)) ))
(progn
(setq sqrtBt (sqrt Bt))
(setq sqrtFh (sqrt (1+ (* Ft h h))))
(setq hsqrtF (* h (sqrt Ft)))
(setq sinh-1 (arcsinh hsqrtF))
(setq f (* sqrtBt (+ sqrtFh (/ sinh-1 hsqrtF )))))
);end_if
(setq si (* b c h wi f))
(setq s (+ s si))
(setq i (1+ i))
)
s
)
(defun arcsinh(x);sinh-1
(log (+ x (sqrt (1+ (* x x)))));ln(x + sqrt(x^2 + 1))
)

(defun c:tt() 
(setq a 2.0 b 3.0 c 4.0 hc 3.0);已知 
(setq hb (- c hc)) 
(setq n 100) 
(setq i 1  h (/  hb c) s 0) 
(repeat n 
(setq wi (/ pi n) 
        ai   (* pi (/ (- (* 2.0 i) 1) (* 2.0 n))) 
         xi  (cos ai) 
;;;         B (…) 
;;;         F (…) 
) 
(setq Bt (- 1.0 (* (- 1.0 (/ (* a a) (* b b))) xi xi) )) 
(setq Ft (- (/ (* a a) (* c c) Bt) 1)) 
   
(if (= Ft 0) (setq f (* 2 h (/ a c))) 
;;;（为了避免程序又臭又长，可以增加中间计算） 
;;;(setq f (1- (* (* (/ (* A A) (* C C))) (/ 1.0 B)) )) 
(progn 
(setq sqrtBt (sqrt Bt)) 
(setq sqrtFh (sqrt (1+ (* Ft h h)))) 
(setq hsqrtF (* h (sqrt Ft))) 
(setq sinh-1 (arcsinh hsqrtF)) 
(setq f (* sqrtBt (+ sqrtFh (/ sinh-1 hsqrtF ))))) 
  );end_if 
(setq si (* b c h wi f)) 
(setq s (+ s si)) 
(setq i (1+ i)) 
) 
s 
) 
(defun arcsinh(x);sinh-1 
  (log (+ x (sqrt (1+ (* x x)))));ln(x + sqrt(x^2 + 1)) 
  )

cad890 · 发表于 2019-7-1 18:24

本帖最后由 cad890 于 2019-7-1 18:26 编辑

既然你这公式都有推论了，直接利用公式计算就可以了。n取得越大，结果越精确。单次循环累加求和就是了。
如：
(setq s 0)
(repeat n
(setq si "公式")
(setq s (+ s si))
)

highflybird · 发表于 2019-7-1 20:10

搜索我的帖子，我有个求积分的程序。

ynhh · 发表于 2019-7-2 08:47

cad890 发表于 2019-7-1 18:24
既然你这公式都有推论了，直接利用公式计算就可以了。n取得越大，结果越精确。单次循环累加求和就是了。
...

谢谢您的指导
只上介绍上的公式太多，并且公式中的很多符号不知何意
人笨实再不好意思

ynhh · 发表于 2019-7-2 08:50

本帖最后由 ynhh 于 2019-7-7 08:29 编辑

highflybird 发表于 2019-7-1 20:10
搜索我的帖子，我有个求积分的程序。

谢谢高飞大师
我也找到研究了您的积分贴子
但介绍中的公式有好几个不知用那个
主要还是公式中的很多符号都不知是什么意思
您在椭圆和积分方面是顶级高手
如您方便请帮我写个实例我也好对照学习一下
如不方便也不能麻烦您
谢谢您了

highflybird · 发表于 2019-7-2 10:43

ynhh 发表于 2019-7-2 08:50
谢谢高飞大师
我也找到研究了您的积分贴子
http://bbs.mjtd.com/forum.php?mod=viewthread&tid=82704&e ...

几个求积分的结果是一致的，只不过效率有区别，一般来说用龙贝格最快。

ynhh · 发表于 2019-7-2 11:31

highflybird 发表于 2019-7-2 10:43
几个求积分的结果是一致的，只不过效率有区别，一般来说用龙贝格最快。

高飞大师您好
能不能就这个实例
写个LISP的实用例子
让大家学习见识一下真实的运用
谢谢您

ynhh · 发表于 2019-7-2 13:38

highflybird 发表于 2019-7-2 12:45
譬如这个求椭圆弧长的例子，就是采用了龙贝格积分。
只要你把积分的函数代入，然后代入下限上限值，就 ...

谢谢您的多次指导
我慢慢研究一下
谢谢

ynhh · 发表于 2019-7-2 13:38

highflybird 发表于 2019-7-2 12:45
譬如这个求椭圆弧长的例子，就是采用了龙贝格积分。
只要你把积分的函数代入，然后代入下限上限值，就 ...

谢谢您的多次指导
我慢慢研究一下
谢谢

账号		自动登录	找回密码
密码			注册

[基础] LISP的积分问题

本帖子中包含更多资源