一元二次、三次、四次方程求解和复数的运算

highflybird 周四, 01/21/2021 - 20:59

对一元三次或者四次方程，是有数学公式求精确解的，可以不用迭代法。参考了维基的上的方法，现在我贴出一元二次、三次或者四次方程的LISP求解方法。使得在求解效率可以得到极大提高。
注明： 因为这几个方程的解有可能是复数，所以我对每个解都用表的形式来列出。
如果这个表的第二项为0，那么这个解是实数，否则是复数。
譬如 ：\({x^4+3x^3+7x^2+2x-5 = 0}\)
 (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)

标签

方程 LISP 复数
编程