一个博士朋友的儿子刚上初中不久,他就提出了下面的一个问题,我感觉很有深度,很难得,现在贴出来说说。
此处用比较纯LISP的方法来解非线性方程。
主要采用Van Wijingaarden-Dekker-Brent方法求解。
继续阅读
对一元三次或者四次方程,是有数学公式求精确解的,可以不用迭代法。参考了维基的上的方法,现在我贴出一元二次、三次或者四次方程的LISP求解方法。使得在求解效率可以得到极大提高。
注明: 因为这几个方程的解有可能是复数,所以我对每个解都用表的形式来列出。
如果这个表的第二项为0,那么这个解是实数,否则是复数。
譬如 :[latex]{x^4+3x^3+7x^2+2x-5 = 0}[/latex]
(Math:Quartic_Equation 1 3 7 2 -5)
==》((-1.19281 -2.21406) (-1.19281 2.21406) (-1.24789 0) (0.633498 0))
意味着这个方程有两个实数解:-1.24789 , 0.633498
两个虚数解:-1.19281-2.21406 i ,-1.19281+2.21406i
另外在末尾附上验算测试函数。
解一元二次方程的LISP源码:
;;;=============================================================
;;;一元二次方程的解
;;;f(x) = a*x^2+b*x+c = 0
;;;Input: the coefficients a, b, c are real numbers
;;;Output: when a /= 0,one or Two solutions,all of them like
;;; this: (x y),means: x + y * i,if it's a real number,
;;; then y = 0.Otherwise , return a real number or nil
;;;Ref: http://en.wikipedia.org/wiki/Quadratic_equation
;;;=============================================================
(defun Math:Quadratic_Equation (a b c / d e g)
(if (zerop a)
(if (not (zerop b))
(list (list (/ (- c) (float b)) 0))
)
(progn
(setq a (float a))
(if (not (equal a 1 1e-14))
(setq b (/ b a)
c (/ c a)
)
)
(setq d (- (* b b) (* 4 c)))
(setq e (* b -0.5))
(cond
( (equal d 0 1e-14)
(list (list e 0) (list e 0))
)
( (> d 0)
(setq g (* (sqrt d) -0.5))
(list (list (- e g) 0) (list (+ e g) 0))
)
( (< d 0)
(setq g (* (sqrt (- d)) -0.5))
(list (list e (- g)) (list e g))
)
)
)
)
)
解一元三次方程的LISP源码:
;;;=============================================================
;;;一元三次方程的解
;;;f(x) = a*x^3+b*x^2+c*x+d = 0
;;;Input: the coefficients a, b, c, d are real numbers
;;;Output: when a /= 0,Three solutions,all of them like this:
;;; (x y),means: x+y*i,if it's a real number,then y = 0;
;;; otherwise goto "Math:Quadratic_Equation"
;;;Ref: http://en.wikipedia.org/wiki/Cubic_function
;;;=============================================================
(defun Math:Cubic_Equation (a b c d / e f g h u w x y)
(cond
( (zerop a)
(Math:Quadratic_Equation b c d)
)
( (zerop d)
(cons '(0 0) (Math:Quadratic_Equation a b c))
)
(t
(setq b (/ b a 3.0))
(setq c (/ c a 6.0))
(setq d (/ d a 2.0))
(setq e (- (* b (- (+ c c c) (* b b))) d)) ;Alpha
(setq f (- (* b b) c c)) ;Beta
(setq g (- (* e e) (* f f f))) ;Delta,The nature of the roots
(cond
( (equal g 0 1e-14)
(setq u (MATH:CubeRoot e))
(list (list (- (+ u u) b) 0.0)
(list (- (+ b u)) 0.0)
(list (- (+ b u)) 0.0)
)
)
( (> g 0)
(setq h (sqrt g))
(setq u (MATH:CubeRoot (+ e h)))
(setq w (MATH:CubeRoot (- e h)))
(setq x (- (+ b (* (+ u w) 0.5))))
(setq y (* (sqrt 3) (- u w) 0.5))
(list (list (+ u w (- b)) 0)
(list x y)
(list x (- y))
)
)
( (< g 0)
(setq x (/ e f (sqrt f)))
(setq y (sqrt (abs (- 1 (* x x)))))
(setq u (/ (atan y x) 3))
(setq w (/ (+ pi pi) 3))
(setq h (* 2 (sqrt f)))
(list (list (- (* (cos u) h) b) 0)
(list (- (* (cos (+ u w)) h) b) 0)
(list (- (* (cos (- u w)) h) b) 0)
)
)
)
)
)
)
解一元四次方程的解:
;;;=============================================================
;;;一元四次方程的解
;;;f(x) = a*x^4+b*x^3+c*x^2+d*x+e= 0
;;;Input: the coefficients a,b,c,d,e are real numbers.
;;;Output: when a /= 0,Three solutions,all of them like this:
;;; (x y),means: x+y*i,if it's a real number,then y = 0,
;;; otherwise goto "Math:Quadratic_Equation".
;;;Ref: http://en.wikipedia.org/wiki/Quartic_function
;;;=============================================================
(defun Math:Quartic_Equation (a b c d e / B2 B3 B4 EPS F G H P Q R V W X Y Z S)
(setq eps 1e-8)
(cond
( (equal a 0 eps)
(Math:Cubic_Equation b c d e)
)
( (equal e 0 eps)
(cons '(0 0) (Math:Cubic_Equation a b c d))
)
( (and (equal b 0 eps) (equal d 0 eps))
(foreach x (Math:Quadratic_Equation a c e)
(foreach y (Math:CSqrt x)
(setq s (cons y s))
)
)
)
( (setq a (float a))
(setq b (/ b a))
(setq c (/ c a))
(setq d (/ d a))
(setq e (/ e a))
(setq b2 (* b b))
(setq b3 (* b2 b))
(setq b4 (* b3 b))
(setq f (+ (* b2 -0.375) c)) ;alpha
(setq g (+ (* b3 0.125) (* b c -0.5) d)) ;beta
(setq h (+ (* -0.01171875 b4) (* 0.0625 c b2) (* -0.25 b d) e))
(if (equal g 0 eps)
(progn
(setq x (* b -0.25))
(mapcar
(function (lambda (e) (list (+ (car e) x) (cadr e))))
(Math:Quartic_Equation 1 0 f 0 h) ;Recursion
)
)
(progn
(setq p (- (* f f 2) h))
(setq q (- (* f f f 0.5) (* f h 0.5) (* g g 0.125)))
(setq r (Math:Cubic_Equation 1 (* 2.5 f) p q))
(while (not (equal (cadar r) 0 eps))
(setq r (cdr r))
)
(setq r (caar r))
(setq w (sqrt (abs (+ f r r))))
(foreach i (list + -)
(foreach j (list + -)
(setq x (i (* b -0.25) (* w 0.5)))
(setq y (- (+ f f f r r (i (/ g 0.5 w)))))
(if (< y 0)
(setq z (list x (j (* (sqrt (- y)) 0.5))))
(setq z (list (+ x (j (* (sqrt y) 0.5))) 0))
)
(setq S (cons z S))
)
)
)
)
)
)
)
复数相关代码:
;;;=============================================================
;;;立方根(为解三次和四次方程编写,因expt函数第一个参数不能为负)
;;;输入:x 实数
;;;输出:x 的立方根
;;;=============================================================
(defun MATH:CubeRoot (x)
(if (< x 0)
(- (expt (- x) 0.33333333333333333333333))
(expt x 0.33333333333333333333333)
)
)
;;;=============================================================
;;;复数加复数
;;;输入:Z1--复数,Z2--复数
;;;输出:复数跟实数相加的结果
;;;=============================================================
(defun Math:C+C (z1 z2)
(mapcar '+ z1 z2)
)
;;;=============================================================
;;;复数减复数
;;;输入:Z1--复数,Z2--复数
;;;输出:复数跟实数相减的结果
;;;=============================================================
(defun math:C-C (z1 z2)
(mapcar '- Z1 Z2)
)
;;;=============================================================
;;;复数乘复数
;;;输入:Z1--复数,Z2--复数
;;;输出:复数跟实数相乘的结果
;;;=============================================================
(defun Math:C*C (Z1 Z2)
(list
(- (* (car Z1) (car z2)) (* (cadr z1) (cadr z2)))
(+ (* (car Z1) (cadr Z2)) (* (cadr z1) (car Z2)))
)
)
;;;=============================================================
;;;复数除以复数
;;;输入:Z1--复数,Z2--复数
;;;输出:复数跟实数相除的结果
;;;=============================================================
(defun Math:C/C (Z1 Z2 / a b c d N)
(setq a (car z1))
(setq b (cadr z1))
(setq c (car z2))
(setq d (cadr z2))
(setq N (float (+ (* c c) (* d d))))
(if (not (zerop N))
(list
(/ (+ (* a c) (* b d)) N)
(/ (- (* b c) (* a d)) N)
)
)
)
;;;=============================================================
;;;复数加实数
;;;输入:Z --复数,R --实数
;;;输出:复数跟实数相加的结果
;;;=============================================================
(defun Math:C+R (Z R)
(list (+ (car z) R) (cadr z))
)
;;;=============================================================
;;;复数减实数
;;;输入:Z --复数,R --实数
;;;输出:复数跟实数相减的结果
;;;=============================================================
(defun Math:C-R (Z R)
(list (- (car z) R) (cadr z))
)
;;;=============================================================
;;;复数乘实数
;;;输入:Z --复数,R --实数
;;;输出:复数跟实数相乘的结果
;;;=============================================================
(defun Math:C*R (Z R)
(list (* (car z) R) (* (cadr z) R))
)
;;;=============================================================
;;;复数除以实数
;;;输入:Z --复数,R --实数
;;;输出:复数跟实数相除的结果
;;;=============================================================
(defun Math:C/R (Z R)
(if (not (zerop R))
(list (/ (car z) R) (/ (cadr z) R))
)
)
;;;=============================================================
;;;复数的模
;;;输入:Z --复数
;;;输出:复数的模,即复数代表的矢量长度
;;;=============================================================
(defun MATH:CNormal (Z)
(distance '(0 0) Z)
)
;;;=============================================================
;;;复数的平方根
;;;输入:Z --复数
;;;输出:复数的平方根,以复数形式表达,有两个解
;;;=============================================================
(defun Math:CSqrt (Z / x y a r)
(setq x (car z))
(setq y (cadr z))
(if (equal y 0 1e-14)
(if (> x 0)
(list (list (setq x (sqrt x)) 0) (list (- x) 0))
(list (list 0 (setq y (sqrt (- x)))) (list 0 (- y)))
)
(progn
(setq a (* (atan y x) 0.5))
(setq r (sqrt (distance '(0 0) Z)))
(setq x (* r (cos a)))
(setq y (* r (sin a)))
(list (list x y) (list (- x) (- y)))
)
)
)
;;;=============================================================
;;;复数的方根
;;;输入:Z --复数, n --正整数
;;;输出:复数的n次方根,以复数形式表达,有n个解。
;;;=============================================================
(defun MATH:CRoot (Z n / r a b i s 2Pi)
(if (and (= (type n) 'INT) (> n 0))
(progn
(setq r (expt (distance '(0 0) z) (/ 1. n)))
(setq a (atan (cadr z) (car z)))
(setq i 0)
(setq 2Pi (+ pi pi))
(repeat n
(setq b (/ (+ a (* i 2Pi)) n))
(setq s (cons (list (* r (cos b)) (* r (sin b))) s))
(setq i (1+ i))
)
(reverse s)
)
)
)
;;;=============================================================
;;;复数的实幂指数
;;;输入:Z --复数, x -- 实数
;;;输出:复数Z 的x次幂,以复数形式表达。
;;;=============================================================
(defun MATH:CRPower (Z x / a r)
(setq a (atan (cadr Z) (car Z)))
(setq r (expt (distance '(0 0) Z) x))
(list (* r (cos (* a x))) (* r (sin (* a x))))
)
;;;=============================================================
;;;复数的复幂指数
;;;输入:Z1 --复数, Z2 -- 复数
;;;输出:复数Z1的Z2次幂,以复数形式表达。
;;;=============================================================
(defun MATH:CZPower (Z1 Z2 / a r x y k)
(if (> (setq r (distance '(0 0) Z1)) 1e-20)
(progn
(setq a (atan (cadr Z1) (car Z1)))
(setq x (car z2))
(setq y (cadr z2))
(setq k (log r))
(setq r (exp (- (* k x) (* y a))))
(setq a (+ (* k y) (* x a)))
(list (* r (cos a)) (* r (sin a)))
)
)
)
;;;=============================================================
;;;复数的自然对数
;;;输入:Z --复数
;;;输出:复数Z 的自然对数,以复数形式表达。
;;;=============================================================
(defun MATH:CExp (Z / r)
(if (> (setq r (distance '(0 0) Z)) 1e-20)
(list (log r) (atan (cadr Z) (car Z)))
)
)
;;;=============================================================
;;;复数的余弦
;;;输入:Z --复数
;;;输出:复数Z 的余弦,以复数形式表达。
;;;=============================================================
(defun MATH:CCOS (Z / x y c s u v)
(setq x (car z))
(setq y (cadr z))
(if (equal y 0 1e-20)
(list (cos x) 0)
(progn
(setq c (* 0.5 (cos x)))
(setq s (* 0.5 (sin x)))
(setq u (exp y))
(setq v (exp (- y)))
(list (* c (+ v u)) (* s (- v u)))
)
)
)
;;;=============================================================
;;;复数的正弦
;;;输入:Z --复数
;;;输出:复数Z 的正弦,以复数形式表达。
;;;=============================================================
(defun MATH:CSIN (Z / x y c s u v)
(setq x (car z))
(setq y (cadr z))
(if (equal y 0 1e-20)
(list (sin x) 0)
(progn
(setq s (* 0.5 (sin x)))
(setq c (* 0.5 (cos x)))
(setq u (exp (- y)))
(setq v (exp y))
(list (* s (+ v u)) (* c (- v u)))
)
)
)
;;;=============================================================
;;;复数的正切
;;;输入:Z --复数
;;;输出:复数Z 的正切,以复数形式表达。
;;;=============================================================
(defun MATH:CTAN (Z / x y c s u v)
(MATH:C/C (MATH:CSIN Z) (MATH:CCOS Z))
)
;;;=============================================================
;;;实系数的复数多项式运算
;;;输入:Z --复数, RealNumbers --实系数列表
;;;输出:根据实系数多项式复数方程值,用复数表示。
;;;=============================================================
(defun MATH:CReal_Polynomial (Z RealNumbers / f)
(cond
( (cdr RealNumbers)
(setq f (MATH:C*R Z (car RealNumbers)))
(repeat (- (length RealNumbers) 2)
(setq RealNumbers (cdr RealNumbers))
(setq f (MATH:C+R f (car RealNumbers)))
(setq f (MATH:C*C f Z))
)
(setq f (MATH:C+R f (cadr RealNumbers)))
)
( (car RealNumbers) (list (car RealNumbers) 0))
)
)
;;;=============================================================
;;;虚系数的复数多项式运算
;;;输入:Z --复数, ImaginaryNumbers--复系数列表
;;;输出:根据复系数多项式复数方程值,用复数表示。
;;;=============================================================
(defun MATH:CImaginary_Polynomial (Z ImaginaryNumbers / f)
(if (setq f (car ImaginaryNumbers))
(progn
(repeat (1- (length ImaginaryNumbers))
(setq ImaginaryNumbers (cdr ImaginaryNumbers))
(setq f (MATH:C*C f Z))
(setq f (MATH:C+C f (car ImaginaryNumbers)))
)
f
)
)
)
测试代码:
;;;=============================================================
;;;以下程序为测试用:
;;;Test for "Math:Quartic_Equation"
;;;If the difference between the Coefficients is big,it will get
;;;some float point error.
;;;=============================================================
(defun C:t4 (/ a b c d e S z z1 z2 z3 z4)
(initget 1)
(setq a (getreal "\nCoefficient a:"))
(initget 1)
(setq b (getreal "\nCoefficient b:"))
(initget 1)
(setq c (getreal "\nCoefficient c:"))
(initget 1)
(setq d (getreal "\nCoefficient d:"))
(initget 1)
(setq e (getreal "\nCoefficient e:"))
;(MISC:test 1000 '((Math:Quartic_Equation a b c d e)))
(setq S (Math:Quartic_Equation a b c d e)) ;get the solutions
(foreach z S
(princ "\n")
(princ (mapcar '(lambda (x) (rtos x 2 20)) z)) ;print them.
)
(foreach z S ;Check every solution
(setq f (MATH:CReal_Polynomial z (list a b c d e)))
(if (not (equal f '(0 0) 1e-6)) ;if f(z) /= 0.0,maybe it's caused by floating point operation;
(alert
(strcat
"Maybe it's a Wrong solution: f("
(vl-princ-to-string z)
") = "
(VL-PRINC-TO-STRING f)
)
)
)
)
(princ)
)
(defun C:t3 (/ a b c d e S z z1 z2 z3 z4)
(initget 1)
(setq a (getreal "\nCoefficient a:"))
(initget 1)
(setq b (getreal "\nCoefficient b:"))
(initget 1)
(setq c (getreal "\nCoefficient c:"))
(initget 1)
(setq d (getreal "\nCoefficient d:"))
(UTI:BENCH
20000
(list
(list 'Math:Cubic_Equation a b c d)
(list 'Math:Cubic_Equation1 a b c d)
)
)
(foreach n '(Math:Cubic_Equation Math:Cubic_Equation1)
(setq S (apply n (list a b c d))) ;get the solutions
(foreach z S
(princ "\n")
(princ (mapcar '(lambda (x) (rtos x 2 20)) z)) ;print them.
)
(foreach z S ;Check every solution
(setq f (MATH:CReal_Polynomial z (list a b c d)))
(if (not (equal f '(0 0) 1e-6)) ;if f(z) /= 0.0,maybe it's caused by floating point operation;
(alert
(strcat
"Use "
(VL-PRINC-TO-STRING n)
",Maybe it's a Wrong solution: f("
(vl-princ-to-string z)
") = "
(VL-PRINC-TO-STRING f)
)
)
)
)
)
(princ)
)
(princ "\nThe command is : test.")
(vl-load-com)
(princ)
在下面的这个帖子中讨论了椭圆的曲率和转弯半径
http://bbs.mjtd.com/thread-62980-1-1.html
现在我把这个主题深化一下,讨论一下曲线的两个函数:
vlax-curve-getSecondDeriv
vlax-curve-getFirstDeriv
这两个函数是什么意思呢?
我们考察AutoCAD里面的曲线类,主要是圆,椭圆,弧和样条曲线,多段线由这几种组合而成。
椭圆和样条曲线实际上都是由参数形成,因此,对于这类曲线,它们每点的坐标可以由参数方程表达:
譬如椭圆 x=a*cos(t); y=b*sin(t);
样条曲线也有方程,假设样条曲线的参数方程为: X= f(t);
Y=g(t);
因此可以对参数方程求导,得到每一点的切线矢量,曲线上每一点对应于一个参数 t0 ,
这个切线矢量的 的X值就是 f(t)在t0处的一阶导数,Y值就是g(t)在t0处的一阶导数,
即( f'(t0), g'(t0),0)
因而我们就理解了vlax-curve-getFirstDeriv 函数返回值的意义,
对于vlax-curve-getSecondDeriv的意义类似,只不过这次换成了二阶导数。
即( f''(t0), g''(t0),0)
那么如何求样条曲线或者椭圆的每一点的曲率及其半径呢?
根据参数方程的求曲率公式,若曲线由下面参数方程给出:
\[\left\{\begin{array}{l}x=\varphi(t) \\ y=\psi(t) \end{array}\right.\]
那么曲率如何计算?提示:
\[K=\frac{\left|\varphi^{\prime}(t) \psi^{\prime \prime}(t)-\varphi^{\prime \prime}(t) \psi^{\prime}(t)\right|}{\left[\varphi^{\prime 2}(t)+\psi^{\prime 2}(t)\right]^{3 / 2}}\]
下面是求值代码:
;;;=============================================================
;;;功能: 获取椭圆上一点处的曲率和离心半径(适用于单次利用此法)
;;;参数: 椭圆实体和椭圆上的一点
;;;返回: 此处离心圆圆心、离心半径及其曲率(离心率)
;;;说明: 如果要在CAD中几何作图,可以参考此贴:
;;; http://bbs.mjtd.com/thread-62980-1-1.html
;;;=============================================================
(defun ELL:GetCurvature (en pt / obj a b px x y v1 v2 rad cen)
(setq obj (vlax-ename->vla-object en))
(setq a (vla-get-MajorRadius obj)) ;椭圆的半长轴
(setq b (vla-get-MinorRadius obj)) ;椭圆的半短轴
(setq pt (vlax-curve-getclosestpointto en pt)) ;保证此点在椭圆上
(setq px (vlax-curve-getParamAtPoint en pt)) ;此点的椭圆参数
(setq v1 (vlax-curve-getFirstDeriv en px)) ;此点的一阶矢量
(setq v2 (list (- (cadr v1)) (car v1) (caddr v1))) ;此点的切线矢量
(setq x (* a (sin px)))
(setq y (* b (cos px)))
(setq rad (/ (expt (+ (* x x) (* y y)) 1.5) (* a b))) ;得到转弯半径
(setq cen (polar pt (angle '(0 0 0) v2) rad)) ;圆心
(list cen rad) ;圆心及半径
)
;;;=============================================================
;;; 利用参数方程的求椭圆的曲率及其离心半径(适用于多次利用此法)
;;; 功能: 获取椭圆上一点处的转弯半径和离心圆圆心
;;; 参数: 椭圆实体和曲线上的一点
;;; 返回: 此处离心圆圆心、离心半径
;;;=============================================================
(defun ELL:GetCurvature1 (en a b pt / px v1 v2 x y rad cen)
(setq pt (vlax-curve-getclosestpointto en pt)) ;保证此点在椭圆上
(setq px (vlax-curve-getParamAtPoint en pt)) ;此点的椭圆参数
(setq v1 (vlax-curve-getFirstDeriv en px)) ;此点的一阶矢量
(setq v2 (list (- (cadr v1)) (car v1) (caddr v1))) ;此点的切线矢量
(setq x (* a (sin px)))
(setq y (* b (cos px)))
(setq rad (/ (expt (+ (* x x) (* y y)) 1.5) (* a b))) ;得到转弯半径
(setq cen (polar pt (angle '(0 0 0) v2) rad)) ;圆心
(list cen rad) ;圆心及半径
)
;;;=============================================================
;;; 一般平面曲线参数方程的曲率离心公式
;;; 功能: 获取曲线上一点处的离心半径和离心圆圆心
;;; 参数: 曲线实体和曲线上的一点
;;; 返回: 此处离心圆圆心、离心半径
;;;=============================================================
(defun CUR:GetCurvature (en pt / ob px v1 v2 v3 x1 y1 x2 y2 cen rad d1 d2)
(setq ob (vlax-ename->vla-object en))
(setq pt (vlax-curve-getclosestpointto en pt)) ;保证此点在曲线上
(setq px (vlax-curve-getParamAtPoint en pt)) ;此点的曲线参数
(setq v1 (vlax-curve-getFirstDeriv en px)) ;此点的一阶矢量
(setq v2 (vlax-curve-getSecondDeriv en px)) ;此点的二阶矢量
(setq v3 (list (- (cadr v1)) (car v1) (caddr v1))) ;此点的切线矢量
(setq x1 (car v1)) ;一阶导数的 X值
(setq y1 (cadr v1)) ;一阶导数的 Y值
(setq x2 (car v2)) ;二阶导数的 X值
(setq y2 (cadr v2)) ;二阶导数的 Y值
(setq d1 (expt (+ (* y1 y1) (* x1 x1)) 1.5))
(setq d2 (- (* x1 y2) (* x2 y1))) ;转弯内外的判定
(if (/= d2 0) ;如果不为直线段
(progn
(setq rad (/ d1 d2))
(if (vlax-method-applicable-p ob 'GetBulge) ;如果为多段线(含圆弧)
(if (< (vla-GetBulge ob (fix px)) 0) ;如果此段凸度小于0
(setq rad (- rad))
)
)
(list (polar pt (angle '(0 0 0) v3) rad) (abs rad)) ;圆心及半径
)
)
)
一些测试代码:
;;;=============================================================
;;;测试程序1: 获取曲线一点处的曲率和离心半径
;;;=============================================================
(defun c:tt1 (/ ent pnt ret)
(setq ent (car (entsel "\n请选取曲线:")))
(while (setq pnt (getpoint "\n点取一点"))
(setq pnt (trans pnt 1 0))
(setq pnt (vlax-curve-getclosestpointto ent pnt))
(setq ret (CUR:GetCurvature ent pnt))
(princ "\n离心半径是:")
(princ ret)
(and ret (apply 'Ent:Make_Circle ret))
)
(princ)
)
;;;=============================================================
;;;测试程序2: 获取椭圆一点处的曲率和离心半径,并比较两种方法
;;;=============================================================
(defun c:tt2 ( / ent obj a b pnt par 1st 2st r1 r2 r3)
(if (and (setq ent (car (entsel "\n请选取椭圆:")))
(setq obj (vlax-ename->vla-object ent))
(vlax-property-available-p obj 'MajorRadius) ;这个地方应该加出错处理
(setq a (vla-get-MajorRadius obj))
(setq b (vla-get-MinorRadius obj))
)
(while (setq pnt (getpoint "\n点取一点:"))
(setq pnt (trans pnt 1 0))
(setq pnt (vlax-curve-getclosestpointto ent pnt))
(setq par (vlax-curve-getparamatpoint ent pnt))
(setq 1st (vlax-curve-getFirstDeriv ent par))
(setq 2st (vlax-curve-getSecondDeriv ent par))
(setq r1 (distance '(0 0) 2st)) ;这个secondDeriv并不意味着半径
(setq r2 (cadr (ELL:GetCurvature1 ent a b pnt)))
(setq r3 (cadr (CUR:getcurvature ent pnt)))
(princ "\nRadius 1 is:")
(princ r1)
(princ "\nRadius 2 is:")
(princ r2)
(princ "\nRadius 3 is:")
(princ r3)
(UTI:Bench
10000
(list
(list 'ELL:GetCurvature1 ent a b pnt)
(list 'CUR:getcurvature ent pnt)
)
)
)
)
)
;;;=============================================================
;;;测试程序3: 动态演示曲线的离心半径
;;;=============================================================
(defun c:tt3 (/ CIR ENT LIN PNT PT0 RET 1ST 2ST CEN PAR PT2 VEC)
(setq ent (car (entsel "\n请选取曲线:")))
(setq cir (Ent:Make_Circle '(0 0 0) 1))
(setq lin (ent:make_Line '(0 0 0) '(0 0 0)))
(setq vec (ent:make_Line '(0 0 0) '(0 0 0)))
(setq lin (vlax-ename->vla-object lin))
(setq vec (vlax-ename->vla-object vec))
(setq cir (vlax-ename->vla-object cir))
(vla-put-color lin 1)
(vla-put-color cir 3)
(vla-put-color vec 6)
(princ "\n按空格或者回车退出!")
(while (= (car (setq pt0 (grread 5 0))) 5)
(setq pnt (trans (cadr pt0) 1 0))
(setq pnt (vlax-curve-getclosestpointto ent pnt))
(setq par (vlax-curve-getparamatpoint ent pnt))
(setq 1st (vlax-curve-getFirstDeriv ent par))
(setq 2st (vlax-curve-getSecondDeriv ent par))
(if (setq ret (CUR:GetCurvature ent pnt))
(progn
(setq pt2 (mapcar '+ pnt 2st))
(setq cen (car ret))
(vlax-put Cir 'Center cen )
(vlax-put cir 'Radius (cadr ret))
(vlax-put lin 'StartPoint cen)
(vlax-put lin 'EndPoint pnt)
(vlax-put vec 'StartPoint pt2)
(vlax-put vec 'EndPoint pnt)
)
)
)
(vla-erase cir)
(vla-erase lin)
(vla-erase vec)
(princ)
)
;;;=============================================================
;;;测试程序4: 由曲线的离心半径描绘其轨迹
;;;=============================================================
(defun c:tt4 (/ ent lst par px1 px2 pxn pnt Inf cen)
(if (setq ent (car (entsel "\n请选取曲线:")))
(progn
(setq px1 (vlax-curve-getStartParam ent))
(setq px2 (vlax-curve-getendparam ent))
(setq pxn (* (- px2 px1) 0.02))
(setq par px1)
(while (<= par px2)
(setq pnt (vlax-curve-getpointatparam ent par))
(setq Inf (CUR:GetCurvature ent pnt))
(if (setq cen (car Inf))
(setq lst (cons cen lst))
)
(setq par (+ par pxn))
)
(setq lst (reverse lst))
(Ent:Make_LWPoly lst 1)
)
)
)
相关联代码:
;;;************************************************************;
;;;实体创建部分
;;;************************************************************;
;;;-------------------------------------------------------------
;;;创建一个点
;;;输入: 一个三维或者二维的点
;;;输出: 点实体的图元名
;;;-------------------------------------------------------------
(defun Ent:Make_Point (p)
(entmakex (list '(0 . "POINT") (cons 10 p)))
)
;;;-------------------------------------------------------------
;;;创建一个带颜色的点(此函数为测试或者其他用途)
;;;输入: 一个三维或者二维的点表和一个颜色号
;;;输出: 点实体的图元名
;;;-------------------------------------------------------------
(defun Ent:MakePoint-1 (p c)
(entmakex (list '(0 . "POINT") (cons 10 p) (cons 62 c)))
)
;;;-------------------------------------------------------------
;;;创建一条直线段
;;;输入: 两个三维或者二维的点
;;;输出: 线段实体的图元名
;;;-------------------------------------------------------------
(defun Ent:Make_Line (p q)
(entmakeX (list '(0 . "LINE") (cons 10 p) (cons 11 q)))
)
;;;-------------------------------------------------------------
;;;创建一个由三条直线组成的三角形
;;;输入: 三个三维或者二维的点
;;;输出: 由三条直线组成的三角形
;;;-------------------------------------------------------------
(defun Ent:Make_Triangle (p1 p2 p3)
(mapcar 'Ent:Make_Line (list p1 p2 p3) (list p2 p3 p1))
)
;;;-------------------------------------------------------------
;;;创建一个三维多段线
;;;输入: 三维的点集
;;;输出: 三维多段线实体
;;;-------------------------------------------------------------
(defun Ent:Make_Poly (pts Closed / e)
(if Closed
(setq e (Entmake (list '(0 . "POLYLINE") '(70 . 9))))
(setq e (Entmake (list '(0 . "POLYLINE") '(70 . 8))))
)
(foreach p pts
(entmake (list '(0 . "VERTEX") '(70 . 32) (cons 10 p)))
)
(entmake '((0 . "SEQEND")))
(entlast)
)
;;;-------------------------------------------------------------
;;;创建轻多段线
;;;输入: 二维的点集
;;;输出: 轻多段线实体名
;;;-------------------------------------------------------------
(defun Ent:Make_LWPoly (pts closed /)
(entmakeX
(append
'((0 . "LWPOLYLINE")
(100 . "AcDbEntity")
(100 . "AcDbPolyline")
)
(list (cons 90 (length pts))) ;顶点个数
(mapcar (function (lambda (x) (cons 10 x))) pts) ;多段线顶点
(list (cons 70 closed)) ;闭合的
)
)
)
;;;-------------------------------------------------------------
;;; 创建圆实体
;;; 输入: 圆心C和半径R
;;; 输出: 圆的实体名
;;;-------------------------------------------------------------
(defun Ent:Make_Circle (C R)
(entmakex (list '(0 . "CIRCLE") (cons 10 C) (cons 40 R)))
)
程序1效果演示:

程序2效果演示:

表达式一般来说有三种:前缀表达式、中缀表达式、后缀表达式,其中后缀表达式又叫做逆波兰表达式。中缀表达式是最符合人们思维方式的一种表达式,顾名思义,就是操作符在操作数的中间。而前缀表达式和后缀表达式中操作符分别在操作数的前面和操作数的后面。在写表达式,我们一般用中缀表达式,譬如 1+2*3-4/5。并且按照操作符的优先级进行计算。
然而LISP语言是一种前缀表达式,为了把表达式转为LISP函数或者求值,需要进行翻译,添加大量的括号和修改函数的顺序。
这个程序的目的就是使得这一工作变简单。
当然,CAD里面本身也有几种种方式能完成这个,但它们的优缺点容我后面讨论。
程序借鉴了飞诗的一些代码,在此深表感谢。
程序的核心代码如下:
;;;=============================================================
;;; 函数目的: 字符表达式转为函数,主要用于多次调用时提升速度
;;; 输入: expr--字符表达式,sFunc--函数名,sArg--参数符号列表
;;; 输出: 定义函数,并返回其名
;;; 例子: (CAL:Expr2Func "sin(x)+20*y" 'test '(x y))
;;; 结果: 定义了一个名为test的函数,参数符号为x y
;;; 注意: 除法区分整数和浮点数,譬如"2/3"和"2/3.0"结果不同;
;;; 可用自定义函数,前提是首先要加载;
;;; 可用科学计算法,但应满足LISP中的语法。建议用括号;
;;; 表达式应满足语法要求,小数和乘号不能按省略写法。
;;;=============================================================
(defun CAL:Expr2Func (expr sFunc sArgs / lst)
(setq lst (CAL:Separate expr)) ;先按照括号空格和运算符分离字符
(setq lst (CAL:Operators lst '((^ . expt)) ())) ;乘方(幂)最优先
(setq lst (CAL:Operators lst '((* . *) (/ . /) (% . rem)) ()));其次乘除和求模运算
(setq lst (CAL:Operators lst '((+ . +) (- . -)) ())) ;最后处理加减法运算
(defun-q-list-set sFunc (cons sArgs lst)) ;表达成函数
sFunc
)
;;;=============================================================
;;; 函数目的: 字符表达式求值
;;; 输入: expr--字符表达式
;;; 输出: 计算表达式的结果
;;; 例子: (CAL:Expr2Value "sin(1)+20*2")
;;; 结果: 40.8415
;;;=============================================================
(defun CAL:Expr2Value (expr / lst)
(setq lst (CAL:Separate expr)) ;先按照括号空格和运算符分离字符
(setq lst (CAL:Operators lst '((^ . expt)) ())) ;乘方(幂)最优先
(setq lst (CAL:Operators lst '((* . *) (/ . /) (% . rem)) ()));其次乘除和求模运算
(setq lst (CAL:Operators lst '((+ . +) (- . -)) ())) ;最后处理加减法运算
(eval (car lst)) ;求值
)
;;;=============================================================
;;; 函数目的: 先分离出函数和+-*/%^运算符,其余均视作变量或数值,
;;; 并简单检查括号匹配。
;;; 输入: expr--字符表达式
;;; 输出: 函数(包括运算符)和变量及数值的列表
;;;=============================================================
(defun CAL:Separate (expr / CHAR FUNS LASTCHAR LST Temp LBRACKET RBRACKET next)
(setq expr (vl-string-translate "{[]}\t\n," "(()) " expr)) ;替换{[]}\t\n,字符
(setq expr (strcase expr t)) ;全部转为小写
(setq funs '("+" "-" "*" "/" "^" "%" )) ;按照基本运算符分割字符
(setq Temp "")
(setq lst "(")
(setq Lbracket 0) ;左括号计数器
(setq Rbracket 0) ;右括号计数器
(while (/= expr "")
(setq char (substr expr 1 1)) ;字符串的第一个字符
(setq next (substr expr 2 1)) ;字符串的第二个字符
(if (or (= char "(")
(= char ")") ;括号一定是分隔符
(and (= char " ") (/= next "(") (/= next " ")) ;如果不是连续的空格符
(and (member char funs) ;根据运算符进行分割
(not (CAL:isScientific temp lastchar char)) ;忽略科学计数法
)
)
(progn
(if (CAL:IsFunction (Read temp)) ;如果为普通函数
(setq lst (strcat lst "(" Temp " " ) ;则把括号移至函数符号前
Lbracket (1+ Lbracket) ;左括号计数器加1
)
(progn
(and (= char "(") (setq Lbracket (1+ Lbracket))) ;左括号计数器加1
(and (= char ")") (setq Rbracket (1+ Rbracket))) ;右括号计数器加1
(setq lst (strcat lst Temp " " char " "))
)
)
(setq Temp "") ;如果是函数或者括号空格之类,则在此处重新开始
)
(setq Temp (strcat Temp char)) ;否则连取这个字符
)
(setq expr (substr expr 2)) ;字符串剩下的字符
(setq lastchar char)
)
(if (/= Lbracket Rbracket) ;如果括号不平衡
(alert "括号不匹配(Mismatched Brackets)!") ;警告信息
(read (strcat lst Temp ")")) ;否则转为表
)
)
;;;=============================================================
;;; 函数目的: 分析+-*/%^运算符,并组合到表中
;;; 输入: lst-已分割的表,funs-待分析的运算符,Recursive-是否递归
;;; 输出: 函数(包括运算符)和变量及数值的列表
;;;=============================================================
(defun CAL:Operators (lst funs Recursive / fun L n)
(foreach a lst
(if (listp a)
(setq a (CAL:Operators a funs T)) ;如果元素为表,则递归进去
)
(if (= (type a) 'INT)
(setq a (float a))
)
(if (setq fun (cdr (assoc (car L) funs))) ;前一个符号为+-*/%^运算符
(if (or (null (setq n (cadr L))) ;前前一个符号为空
(and (VL-SYMBOLP n) (CAL:IsFunction n)) ;或者是函数符号
)
(setq L (cons (list fun a) (cdr L))) ;无须交换位置
(setq L (cons (list fun n a) (cddr L))) ;交换运算符和操作数位置
)
(setq L (cons a L)) ;其他的不做改变
)
)
(setq n (car L))
(if (and Recursive (not (cadr L)) (or (listp n) (numberp n))) ;如果是递归的,而且只有一个元素,且这个元素为表或者数字
n ;那么就只取这个元素,以防止多余括号出现
(reverse L) ;cons运算后的反转表列
)
)
;;;=============================================================
;;; 函数目的: 判断一个符号是不是普通函数(内部函数或自定义函数)
;;;=============================================================
(defun CAL:IsFunction (n / s)
(setq s (type (eval n)))
(and (or (= s 'SUBR) (= s 'USUBR)) (not (member n '(+ - * /))))
)
;;;=============================================================
;;; 函数目的: 检测一个字符串是否是科学计数法(是否有更好方法?)
;;;=============================================================
(defun CAL:isScientific (temp lastchar char)
(and (= lastchar "e") (numberp (read (strcat temp char "0"))))
)
;;;=============================================================
;;; 函数目的: 检查函数表达式转函数的结果
;;; 输入: lst,用cal:expr2func求得的表
;;; 输出: 如果表达式里有非参数且未赋值的变量符号则返回nil, 否则T
;;; 例子: (CAL:CheckFunc (CAL:Expr2func "sin(a)+20*2" 'fx '(x)))
;;; 结果: nil
;;;=============================================================
(defun CAL:CheckFunc (lst / isOK CAL:TempSym Args)
(setq IsOK T)
(setq Args (car lst))
(while (setq lst (cdr lst))
(setq CAL:TempSym (car lst)) ;对表中的每个元素
(if (listp CAL:TempSym) ;如果这个元素为表
(if CAL:TempSym
(setq IsOk (CAL:CheckFunc (cons Args CAL:TempSym))) ;且不为空则递归进去
(setq IsOk nil) ;否则检测结果为假
)
(if (and (vl-symbolp CAL:TempSym) ;如果是一个符号
(not (member CAL:TempSym Args)) ;且不为参数表中的符号
(not (vl-symbol-value CAL:TempSym)) ;且未赋值
)
(setq IsOk nil) ;则检测结果为假
)
)
(if (null IsOK)
(setq lst nil)
)
)
IsOK
)
;;;=============================================================
;;;以下函数为自定义的一些简单的数学函数
;;;=============================================================
(defun r2d (x) (* 57.2957795130823208768 x)) ;弧度转度
(defun d2r (x) (* 0.01745329251994329577 x)) ;度转弧度
(defun int (x) (atoi (rtos x 2 0))) ;四舍五入取整函数
(defun ceil (x) (1+ (fix x))) ;天花板函数
(defun ln (x) (log x)) ;以e为底的对数函数
(defun log10 (x) (* (log x) 0.43429448190325182765)) ;以10为底的对数函数
(defun exp10 (x) (expt 10 x)) ;以10为底的指数函数
(defun pow (x y) (expt x y)) ;指数函数
(defun tan (x) (/ (sin x) (cos x))) ;正切函数
(defun cot (x) (/ (cos x) (sin x))) ;余切函数
(defun sec (x) (/ 1 (cos x))) ;正割函数
(defun csc (x) (/ 1 (sin x))) ;余割函数
(defun asin (x) (atan x (sqrt (- 1 (* x x))))) ;反正弦函数
(defun acos (x) (atan (sqrt (- 1 (* x x))) x)) ;反余弦函数
(defun sinh (x) (* 0.5 (- (exp x) (exp (- x))))) ;双曲正弦函数
(defun cosh (x) (* 0.5 (+ (exp x) (exp (- x))))) ;双曲余弦函数
(defun tanh (x) (- 1 (/ 2 (1+ (exp (+ x x)))))) ;双曲正切函数
(defun coth (x) (/ 1 (tanh x))) ;双曲余切函数
(defun sech (x) (/ 1 (cosh x))) ;双曲正割函数
(defun csch (x) (/ 1 (sinh x))) ;双曲余割函数
(defun asinh (x) (log (+ x (sqrt (1+ (* x x)))))) ;反双曲正弦函数=log(x+sqrt(x*x+1))
(defun acosh (x) (log (+ x (sqrt (1- (* x x)))))) ;反双曲余弦函数=log(x+sqrt(x*x-1))
(defun atanh (x) (log (sqrt (/ (+ 1 x)(- 1 x))))) ;反双曲正切函数=log(sqrt((1+x)/(1-x)))
(defun revSign (x) (- x)) ;反号函数
(defun reciprocal (x) (/ 1.0 x)) ;倒数
(defun sqr (x) (* x x)) ;平方函数
(defun cube (x) (* x x x)) ;立方函数
(defun cuberoot (x) ;立方根函数
(if (minusp x)
(- (expt (- x) 0.333333333333333333333))
(expt x 0.333333333333333333333)
)
)
(defun round (x / y) ;四舍五入函数
(setq y (fix x))
(if (< (abs (- x y)) 0.5)
y
(if (< x 0)
(1- y)
(1+ y)
)
)
)
以下是一些测试:
;;; 例子:
;;; (CAL:Separate "(sin(-x)-cos(-x+(1+8*(2/7))+2^4-5))*0.5-0.5e-20+20*cos(x)+20")
;;; 结果: ((SIN - X) - (COS - X + (1 + 8 * (2 / 7)) + 2 ^ 4 - 5))
;;; (CAL:Expr2Func "(sin(+x)-cos(-x+(1+8*(2/7))+(2^4)-5))*0.5-0.5e-20+20*cos(x)+20" 'test '(x))
;;; 结果: 定义了一个名为test的函数,参数符号为x
;;; (CAL:Expr2Value "(sin(+0.5)-cos(-pi+(1+8*(2/7))+(2^4)-5))*0.5-0.5e-20+20*cos(pi/2)+20")
;;; 结果: 20.6616
;;; 以下是关于这个程序的其他方法:
;;; 方法一:用cal函数计算
;;; 如:(cal "1+4+5*2+(5+5)/2+((6+6)/2+(5+5)/2)")
;;; 优点:CAD内置函数。
;;; 缺点:这个函数要求先要加载cal函数.并且三角函数会自动把变量或者数值理解为角度。
;;; 方法二:wcs脚本语言法,无痕提出的一种方法
;;; (setq wcs (vla-GetInterfaceObject (vlax-get-acad-object) "ScriptControl"))
;;; (vlax-put-property wcs "language" "vbs")
;;; (vla-eval wcs "1+4+5*2+(5+5)/2+((6+6)/2+(5+5)/2)") ;返回 ->31.0
;;; 优点:能按照vb的语法直接计算。
;;; 缺点:难以定义表达式为函数,不能利用自定义函数,在64位CAD上此法行不通,因为不能创建脚本对象。
;;; 下面例子为在CAD中绘制函数图像
(defun c:test1(/ expr a b d x y e pts)
(setq expr (getstring "\n请输入表达式:"))
(initget 1)
(setq a (getreal "\n上届:"))
(initget 1)
(setq b (getreal "\n下届:"))
(if (CAL:EXPR2FUNC expr 'test '(x))
(progn
(setq d (/ (- b a) 1000.0))
(setq x a)
(setq pts nil)
(repeat 1000
(setq x (+ x d))
(setq y (test x))
(setq pts (cons (list x y 0) pts))
)
(setq pts (reverse pts))
(setq e (Entmake (list '(0 . "POLYLINE") '(70 . 8))))
(foreach p pts
(entmake (list '(0 . "VERTEX") '(70 . 32) (cons 10 p)))
)
(entmake '((0 . "SEQEND")))
(entlast)
)
)
)
;;; 在CAD中绘制参数曲线
;;; x=a*(2*cos(t)-cos(2*t))
;;; y=a*(2*sin(t)-sin(2*t))
(defun c:test2 (/ expr1 expr2 a b d k x y pts e)
(setq expr1 "3*(2*cos(k)-cos(2*k))")
(setq expr2 "3*(2*sin(k)-sin(2*k))")
(setq a 0)
(setq b (+ pi pi))
(CAL:EXPR2FUNC expr1 'fx '(k))
(CAL:EXPR2FUNC expr2 'fy '(k))
(setq d (/ (- b a) 360))
(setq k a)
(setq pts nil)
(repeat 360
(setq k (+ k d))
(setq x (fx k))
(setq y (fy k))
(setq pts (cons (list x y 0) pts))
)
(setq pts (reverse pts))
(setq e (Entmake (list '(0 . "POLYLINE") '(70 . 9))))
(foreach p pts
(entmake (list '(0 . "VERTEX") '(70 . 32) (cons 10 p)))
)
(entmake '((0 . "SEQEND")))
(entlast)
)
;;; 定义为函数后,明显速度快多了
(defun c:test3(/ str1 str2 x)
(setq str1 "(sin(+x)-cos(-x+(1+8*(2/7.0))+(2^4)-5))*0.5-0.5e-20+20*cos(x)+20")
(setq str2 "(sin(r2d(x))-cos(r2d(-x+(1+8*(2/7.0))+(2^4)-5)))*0.5-0.5e-20+20*cos(r2d(x))+20")
(CAL:Expr2Func str1 'f1 '(x))
(setq x 12)
(uti:bench 1000
(list
(cons 'f1 '(12))
(cons 'CAL:Expr2Value (list str1))
(cons 'cal (list str2))
)
)
)