LISP的积分问题
本帖最后由 ynhh 于 2019-7-7 08:06 编辑椭球体球缺的表面积计算
用积分公式用LISP程序表达并计算出来
这是一个网上的
有结果和公式
本帖最后由 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可以满足精度
)
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 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 1h (/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:26 编辑
既然你这公式都有推论了，直接利用公式计算就可以了。n取得越大，结果越精确。单次循环累加求和就是了。
如：
(setq s 0)
(repeat n
(setq si "公式")
(setq s (+ s si))
)
搜索我的帖子，我有个求积分的程序。 cad890 发表于 2019-7-1 18:24
既然你这公式都有推论了，直接利用公式计算就可以了。n取得越大，结果越精确。单次循环累加求和就是了。
...
谢谢您的指导
只上介绍上的公式太多，并且公式中的很多符号不知何意
人笨实再不好意思 本帖最后由 ynhh 于 2019-7-7 08:29 编辑
highflybird 发表于 2019-7-1 20:10
搜索我的帖子，我有个求积分的程序。
谢谢高飞大师
我也找到研究了您的积分贴子
但介绍中的公式有好几个不知用那个
主要还是公式中的很多符号都不知是什么意思
您在椭圆和积分方面是顶级高手
如您方便请帮我写个实例我也好对照学习一下
如不方便也不能麻烦您
谢谢您了 ynhh 发表于 2019-7-2 08:50
谢谢高飞大师
我也找到研究了您的积分贴子
http://bbs.mjtd.com/forum.php?mod=viewthread&tid=82704&e ...
几个求积分的结果是一致的，只不过效率有区别，一般来说用龙贝格最快。 highflybird 发表于 2019-7-2 10:43
几个求积分的结果是一致的，只不过效率有区别，一般来说用龙贝格最快。
高飞大师您好
能不能就这个实例
写个LISP的实用例子
让大家学习见识一下真实的运用
谢谢您 highflybird 发表于 2019-7-2 12:45
譬如这个求椭圆弧长的例子，就是采用了龙贝格积分。
只要你把积分的函数代入，然后代入下限上限值，就 ...
谢谢您的多次指导
我慢慢研究一下
谢谢 highflybird 发表于 2019-7-2 12:45
譬如这个求椭圆弧长的例子，就是采用了龙贝格积分。
只要你把积分的函数代入，然后代入下限上限值，就 ...
谢谢您的多次指导
我慢慢研究一下
谢谢