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

对一元三次或者四次方程,是有数学公式求精确解的,可以不用迭代法。参考了维基的上的方法,现在我贴出一元二次、三次或者四次方程的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源码:

  1. ;;;=============================================================
  2. ;;;一元二次方程的解                                     	
  3. ;;;f(x) = a*x^2+b*x+c = 0                               	
  4. ;;;Input: the coefficients a, b, c are real numbers     	
  5. ;;;Output: when a /= 0,one or Two solutions,all of them like 	
  6. ;;;        this: (x y),means: x + y * i,if it's a real number,  
  7. ;;;        then y = 0.Otherwise ,  return a real number or nil  
  8. ;;;Ref: http://en.wikipedia.org/wiki/Quadratic_equation 	
  9. ;;;=============================================================
  10. (defun Math:Quadratic_Equation (a b c / d e g)
  11.   (if (zerop a)
  12.     (if (not (zerop b))
  13.       (list (list (/ (- c) (float b)) 0))
  14.     )
  15.     (progn
  16.       (setq a (float a))
  17.       (if (not (equal a 1 1e-14))
  18.         (setq b (/ b a)
  19.               c (/ c a)
  20.         )
  21.       )
  22.       (setq d (- (* b b) (* 4 c)))
  23.       (setq e (* b -0.5))
  24.       (cond
  25.         ( (equal d 0 1e-14)
  26.           (list (list e 0) (list e 0))
  27.         )
  28.         ( (> d 0)
  29.           (setq g (* (sqrt d) -0.5))
  30.           (list (list (- e g) 0) (list (+ e g) 0))
  31.         )
  32.         ( (< d 0)
  33.           (setq g (* (sqrt (- d)) -0.5))
  34.           (list (list e (- g)) (list e g))
  35.         )
  36.       )
  37.     )
  38.   )
  39. )

解一元三次方程的LISP源码:

  1. ;;;=============================================================
  2. ;;;一元三次方程的解                                     	
  3. ;;;f(x) = a*x^3+b*x^2+c*x+d = 0                         	
  4. ;;;Input: the coefficients a, b, c, d are real numbers  	
  5. ;;;Output: when a /= 0,Three solutions,all of them like this: 	
  6. ;;;        (x y),means: x+y*i,if it's a real number,then y = 0; 
  7. ;;;        otherwise goto "Math:Quadratic_Equation"     	
  8. ;;;Ref: http://en.wikipedia.org/wiki/Cubic_function     	
  9. ;;;=============================================================
  10. (defun Math:Cubic_Equation (a b c d / e f g h u w x y)
  11.   (cond
  12.     ( (zerop a)
  13.       (Math:Quadratic_Equation b c d)
  14.     )
  15.     ( (zerop d)
  16.       (cons '(0 0) (Math:Quadratic_Equation a b c))
  17.     )
  18.     (t
  19.       (setq b (/ b a 3.0))
  20.       (setq c (/ c a 6.0))
  21.       (setq d (/ d a 2.0))
  22.  
  23.       (setq e (- (* b (- (+ c c c) (* b b))) d))        	;Alpha
  24.       (setq f (- (* b b) c c))                          	;Beta
  25.       (setq g (- (* e e) (* f f f)))                    	;Delta,The nature of the roots
  26.       (cond
  27.         ( (equal g 0 1e-14)
  28.           (setq u (MATH:CubeRoot e))
  29.           (list (list (- (+ u u) b) 0.0)
  30.                 (list (- (+ b u)) 0.0)
  31.                 (list (- (+ b u)) 0.0)
  32.           )
  33.         )
  34.         ( (> g 0)
  35.           (setq h (sqrt g))
  36.           (setq u (MATH:CubeRoot (+ e h)))
  37.           (setq w (MATH:CubeRoot (- e h)))
  38.           (setq x (- (+ b (* (+ u w) 0.5))))
  39.           (setq y (* (sqrt 3) (- u w) 0.5))
  40.           (list (list (+ u w (- b)) 0)
  41.                 (list x y)
  42.                 (list x (- y))
  43.           )
  44.         )
  45.         ( (< g 0)
  46.           (setq x (/ e f (sqrt f)))
  47.           (setq y (sqrt (abs (- 1 (* x x)))))
  48.  
  49.           (setq u (/ (atan y x) 3))
  50.           (setq w (/ (+ pi pi) 3))
  51.           (setq h (* 2 (sqrt f)))
  52.           (list (list (- (* (cos u) h) b) 0)
  53.                 (list (- (* (cos (+ u w)) h) b) 0)
  54.                 (list (- (* (cos (- u w)) h) b) 0)
  55.           )
  56.         )
  57.       )
  58.     )
  59.   )
  60. )

解一元四次方程的解:

  1. ;;;=============================================================
  2. ;;;一元四次方程的解                                     	
  3. ;;;f(x) = a*x^4+b*x^3+c*x^2+d*x+e= 0                    	
  4. ;;;Input: the coefficients a,b,c,d,e are real numbers.  	
  5. ;;;Output: when a /= 0,Three solutions,all of them like this: 	
  6. ;;;        (x y),means: x+y*i,if it's a real number,then y = 0, 
  7. ;;;        otherwise goto "Math:Quadratic_Equation".            
  8. ;;;Ref: http://en.wikipedia.org/wiki/Quartic_function   	
  9. ;;;=============================================================
  10. (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)
  11.   (setq eps 1e-8)
  12.   (cond
  13.     ( (equal a 0 eps)
  14.       (Math:Cubic_Equation b c d e)
  15.     )
  16.     ( (equal e 0 eps)
  17.       (cons '(0 0) (Math:Cubic_Equation a b c d))
  18.     )
  19.     ( (and (equal b 0 eps) (equal d 0 eps))
  20.       (foreach x (Math:Quadratic_Equation a c e)
  21.         (foreach y (Math:CSqrt x)
  22.           (setq s (cons y s))
  23.         )
  24.       )
  25.     )  
  26.     ( (setq a (float a))
  27.       (setq b (/ b a))
  28.       (setq c (/ c a))
  29.       (setq d (/ d a))
  30.       (setq e (/ e a))
  31.       (setq b2 (* b b))
  32.       (setq b3 (* b2 b))
  33.       (setq b4 (* b3 b))
  34.  
  35.       (setq f (+ (* b2 -0.375) c))                      	;alpha
  36.       (setq g (+ (* b3 0.125) (* b c -0.5) d))          	;beta
  37.       (setq h (+ (* -0.01171875 b4) (* 0.0625 c b2) (* -0.25 b d) e))
  38.  
  39.       (if (equal g 0 eps)
  40.         (progn 
  41.           (setq x (* b -0.25))
  42.           (mapcar
  43.             (function (lambda (e) (list (+ (car e) x) (cadr e))))
  44.             (Math:Quartic_Equation 1 0 f 0 h)           	;Recursion 
  45.           )
  46.         )
  47.         (progn 
  48.           (setq p (- (* f f 2) h))
  49.           (setq q (- (* f f f 0.5) (* f h 0.5) (* g g 0.125)))
  50.           (setq r (Math:Cubic_Equation 1 (* 2.5 f) p q))
  51.           (while (not (equal (cadar r) 0 eps))
  52.             (setq r (cdr r))
  53.           )
  54.           (setq r (caar r))
  55.           (setq w (sqrt (abs (+ f r r))))
  56.           (foreach i (list + -)
  57.             (foreach j (list + -)
  58.               (setq x (i (* b -0.25) (* w 0.5)))
  59.               (setq y (- (+ f f f r r (i (/ g 0.5 w)))))
  60.               (if (< y 0)
  61.                 (setq z (list x (j (* (sqrt (- y)) 0.5))))
  62.                 (setq z (list (+ x (j (* (sqrt y) 0.5))) 0)) 
  63.               )
  64.               (setq S (cons z S))
  65.             )
  66.           )
  67.         )
  68.       )
  69.     )
  70.   )
  71. )

复数相关代码:

  1. ;;;=============================================================
  2. ;;;立方根(为解三次和四次方程编写,因expt函数第一个参数不能为负) 
  3. ;;;输入:x 实数                                                 
  4. ;;;输出:x 的立方根                                             
  5. ;;;=============================================================
  6. (defun MATH:CubeRoot (x)
  7.   (if (< x 0)
  8.     (- (expt (- x) 0.33333333333333333333333))
  9.     (expt x 0.33333333333333333333333)
  10.   )
  11. )
  12.  
  13. ;;;=============================================================
  14. ;;;复数加复数                                                   
  15. ;;;输入:Z1--复数,Z2--复数                                     
  16. ;;;输出:复数跟实数相加的结果                                   
  17. ;;;=============================================================
  18. (defun Math:C+C (z1 z2)
  19.   (mapcar '+ z1 z2)
  20. )
  21.  
  22. ;;;=============================================================
  23. ;;;复数减复数                                                   
  24. ;;;输入:Z1--复数,Z2--复数                                     
  25. ;;;输出:复数跟实数相减的结果                                   
  26. ;;;=============================================================
  27. (defun math:C-C (z1 z2)
  28.   (mapcar '- Z1 Z2)
  29. )
  30.  
  31. ;;;=============================================================
  32. ;;;复数乘复数                                                   
  33. ;;;输入:Z1--复数,Z2--复数                             	
  34. ;;;输出:复数跟实数相乘的结果                           	
  35. ;;;=============================================================
  36. (defun Math:C*C (Z1 Z2)
  37.   (list
  38.     (- (* (car Z1) (car z2)) (* (cadr z1) (cadr z2)))
  39.     (+ (* (car Z1) (cadr Z2)) (* (cadr z1) (car Z2)))
  40.   )
  41. )
  42.  
  43. ;;;=============================================================
  44. ;;;复数除以复数                                         	
  45. ;;;输入:Z1--复数,Z2--复数                             	
  46. ;;;输出:复数跟实数相除的结果                           	
  47. ;;;=============================================================
  48. (defun Math:C/C (Z1 Z2 / a b c d N)
  49.   (setq a (car  z1))
  50.   (setq b (cadr z1))
  51.   (setq c (car  z2))
  52.   (setq d (cadr z2))
  53.   (setq N (float (+ (* c c) (* d d))))
  54.   (if (not (zerop N))
  55.     (list
  56.       (/ (+ (* a c) (* b d)) N)
  57.       (/ (- (* b c) (* a d)) N)
  58.     )
  59.   )
  60. )
  61.  
  62. ;;;=============================================================
  63. ;;;复数加实数                                           	
  64. ;;;输入:Z --复数,R --实数                             	
  65. ;;;输出:复数跟实数相加的结果                           	
  66. ;;;=============================================================
  67. (defun Math:C+R (Z R)
  68.   (list (+ (car z) R) (cadr z))
  69. )
  70.  
  71. ;;;=============================================================
  72. ;;;复数减实数                                           	
  73. ;;;输入:Z --复数,R --实数                             	
  74. ;;;输出:复数跟实数相减的结果                           	
  75. ;;;=============================================================
  76. (defun Math:C-R (Z R)
  77.   (list (- (car z) R) (cadr z))
  78. )
  79.  
  80. ;;;=============================================================
  81. ;;;复数乘实数                                           	
  82. ;;;输入:Z --复数,R --实数                             	
  83. ;;;输出:复数跟实数相乘的结果                           	
  84. ;;;=============================================================
  85. (defun Math:C*R (Z R)
  86.   (list (* (car z) R) (* (cadr z) R))
  87. )
  88.  
  89. ;;;=============================================================
  90. ;;;复数除以实数                                         	
  91. ;;;输入:Z --复数,R --实数                             	
  92. ;;;输出:复数跟实数相除的结果                           	
  93. ;;;=============================================================
  94. (defun Math:C/R (Z R)
  95.   (if (not (zerop R))
  96.     (list (/ (car z) R) (/ (cadr z) R))
  97.   )
  98. )
  99.  
  100. ;;;=============================================================
  101. ;;;复数的模                                             	
  102. ;;;输入:Z --复数                                       	
  103. ;;;输出:复数的模,即复数代表的矢量长度                 	
  104. ;;;=============================================================
  105. (defun MATH:CNormal (Z)
  106.   (distance '(0 0) Z)
  107. )
  108.  
  109. ;;;=============================================================
  110. ;;;复数的平方根                                         	
  111. ;;;输入:Z --复数                                       	
  112. ;;;输出:复数的平方根,以复数形式表达,有两个解         	
  113. ;;;=============================================================
  114. (defun Math:CSqrt (Z / x y a r)
  115.   (setq x (car z))
  116.   (setq y (cadr z))
  117.   (if (equal y 0 1e-14)
  118.     (if (> x 0)
  119.       (list (list (setq x (sqrt x)) 0) (list (- x) 0))
  120.       (list (list 0 (setq y (sqrt (- x)))) (list 0 (- y)))
  121.     )
  122.     (progn
  123.       (setq a (* (atan y x) 0.5))
  124.       (setq r (sqrt (distance '(0 0) Z)))
  125.       (setq x (* r (cos a)))
  126.       (setq y (* r (sin a)))
  127.       (list (list x y) (list (- x) (- y)))
  128.     )
  129.   )
  130. )
  131.  
  132. ;;;=============================================================
  133. ;;;复数的方根                                           	
  134. ;;;输入:Z --复数, n --正整数                           	
  135. ;;;输出:复数的n次方根,以复数形式表达,有n个解。       	
  136. ;;;=============================================================
  137. (defun MATH:CRoot (Z n / r a b i s 2Pi)
  138.   (if (and (= (type n) 'INT) (> n 0))
  139.     (progn
  140.       (setq r (expt (distance '(0 0) z) (/ 1. n)))
  141.       (setq a (atan (cadr z) (car z)))
  142.       (setq i 0)
  143.       (setq 2Pi (+ pi pi))
  144.       (repeat n
  145.         (setq b (/ (+ a (* i 2Pi)) n))
  146.         (setq s (cons (list (* r (cos b)) (* r (sin b))) s))
  147.         (setq i (1+ i))
  148.       )
  149.       (reverse s)
  150.     )
  151.   )
  152. )
  153.  
  154. ;;;=============================================================
  155. ;;;复数的实幂指数                                       	
  156. ;;;输入:Z --复数, x -- 实数                            	
  157. ;;;输出:复数Z 的x次幂,以复数形式表达。                	
  158. ;;;=============================================================
  159. (defun MATH:CRPower (Z x / a r)
  160.   (setq a (atan (cadr Z) (car Z)))
  161.   (setq r (expt (distance '(0 0) Z) x))
  162.   (list (* r (cos (* a x))) (* r (sin (* a x))))
  163. )
  164.  
  165. ;;;=============================================================
  166. ;;;复数的复幂指数                                       	
  167. ;;;输入:Z1 --复数, Z2 -- 复数                          	
  168. ;;;输出:复数Z1的Z2次幂,以复数形式表达。               	
  169. ;;;=============================================================
  170. (defun MATH:CZPower (Z1 Z2 / a r x y k)
  171.   (if (> (setq r (distance '(0 0) Z1)) 1e-20)
  172.     (progn 
  173.       (setq a (atan (cadr Z1) (car Z1)))
  174.       (setq x (car z2))
  175.       (setq y (cadr z2))
  176.       (setq k (log r))
  177.       (setq r (exp (- (* k x) (* y a))))
  178.       (setq a (+ (* k y) (* x a)))
  179.       (list (* r (cos a)) (* r (sin a)))
  180.     )
  181.   )
  182. )
  183.  
  184. ;;;=============================================================
  185. ;;;复数的自然对数                                        	
  186. ;;;输入:Z --复数                                       	
  187. ;;;输出:复数Z 的自然对数,以复数形式表达。             	
  188. ;;;=============================================================
  189. (defun MATH:CExp (Z / r)
  190.   (if (> (setq r (distance '(0 0) Z)) 1e-20)
  191.     (list (log r) (atan (cadr Z) (car Z)))
  192.   )
  193. )
  194.  
  195. ;;;=============================================================
  196. ;;;复数的余弦                                           	
  197. ;;;输入:Z --复数                                          	
  198. ;;;输出:复数Z 的余弦,以复数形式表达。                 	
  199. ;;;=============================================================
  200. (defun MATH:CCOS (Z / x y c s u v)
  201.   (setq x (car z))
  202.   (setq y (cadr z))
  203.   (if (equal y 0 1e-20)
  204.     (list (cos x) 0)
  205.     (progn
  206.       (setq c (* 0.5 (cos x)))
  207.       (setq s (* 0.5 (sin x)))
  208.       (setq u (exp y))
  209.       (setq v (exp (- y)))
  210.       (list (* c (+ v u)) (* s (- v u)))
  211.     )
  212.   )
  213. )
  214.  
  215. ;;;=============================================================
  216. ;;;复数的正弦                                          		
  217. ;;;输入:Z --复数                                       	
  218. ;;;输出:复数Z 的正弦,以复数形式表达。                 	
  219. ;;;=============================================================
  220. (defun MATH:CSIN (Z / x y c s u v)
  221.   (setq x (car z))
  222.   (setq y (cadr z))
  223.   (if (equal y 0 1e-20)
  224.     (list (sin x) 0)
  225.     (progn
  226.       (setq s (* 0.5 (sin x)))
  227.       (setq c (* 0.5 (cos x)))
  228.       (setq u (exp (- y)))
  229.       (setq v (exp y))
  230.       (list (* s (+ v u)) (* c (- v u)))
  231.     )
  232.   )
  233. )
  234.  
  235. ;;;=============================================================
  236. ;;;复数的正切                                           	
  237. ;;;输入:Z --复数                                       	
  238. ;;;输出:复数Z 的正切,以复数形式表达。                 	
  239. ;;;=============================================================
  240. (defun MATH:CTAN (Z / x y c s u v)
  241.   (MATH:C/C (MATH:CSIN Z) (MATH:CCOS Z))
  242. )
  243.  
  244. ;;;=============================================================
  245. ;;;实系数的复数多项式运算                               	
  246. ;;;输入:Z --复数, RealNumbers --实系数列表             	
  247. ;;;输出:根据实系数多项式复数方程值,用复数表示。       	
  248. ;;;=============================================================
  249. (defun MATH:CReal_Polynomial (Z RealNumbers / f)
  250.   (cond
  251.     ( (cdr RealNumbers) 
  252.       (setq f (MATH:C*R Z (car RealNumbers)))
  253.       (repeat (- (length RealNumbers) 2)
  254.         (setq RealNumbers (cdr RealNumbers))
  255.         (setq f (MATH:C+R f (car RealNumbers)))
  256.         (setq f (MATH:C*C f Z))
  257.       )
  258.       (setq f (MATH:C+R f (cadr RealNumbers)))
  259.     )
  260.     ( (car RealNumbers) (list (car RealNumbers) 0))
  261.   )
  262. )
  263.  
  264. ;;;=============================================================
  265. ;;;虚系数的复数多项式运算                               	
  266. ;;;输入:Z --复数, ImaginaryNumbers--复系数列表         	
  267. ;;;输出:根据复系数多项式复数方程值,用复数表示。       	
  268. ;;;=============================================================
  269. (defun MATH:CImaginary_Polynomial (Z ImaginaryNumbers / f)
  270.   (if (setq f (car ImaginaryNumbers))
  271.     (progn
  272.       (repeat (1- (length ImaginaryNumbers))
  273.         (setq ImaginaryNumbers (cdr ImaginaryNumbers))
  274.         (setq f (MATH:C*C f Z)) 
  275.         (setq f (MATH:C+C f (car ImaginaryNumbers)))
  276.       )
  277.       f
  278.     )
  279.   )
  280. )

测试代码:

  1. ;;;=============================================================
  2. ;;;以下程序为测试用:                                    	
  3. ;;;Test for "Math:Quartic_Equation"                     	
  4. ;;;If the difference between the Coefficients is big,it will get
  5. ;;;some float point error.                   			
  6. ;;;=============================================================
  7. (defun C:t4 (/ a b c d e S z z1 z2 z3 z4)
  8.   (initget 1)
  9.   (setq a (getreal "\nCoefficient a:"))
  10.   (initget 1)
  11.   (setq b (getreal "\nCoefficient b:"))
  12.   (initget 1)
  13.   (setq c (getreal "\nCoefficient c:"))
  14.   (initget 1)
  15.   (setq d (getreal "\nCoefficient d:"))
  16.   (initget 1)
  17.   (setq e (getreal "\nCoefficient e:"))
  18.   ;(MISC:test 1000 '((Math:Quartic_Equation a b c d e)))
  19.   (setq S (Math:Quartic_Equation a b c d e))            	;get the solutions
  20.   (foreach z S
  21.     (princ "\n")
  22.     (princ (mapcar '(lambda (x) (rtos x 2 20)) z))      	;print them.
  23.   )
  24.   (foreach z S                                          	;Check every solution
  25.     (setq f (MATH:CReal_Polynomial z (list a b c d e)))
  26.     (if (not (equal f '(0 0) 1e-6))                     	;if f(z) /= 0.0,maybe it's caused by floating point operation;
  27.       (alert
  28.         (strcat
  29.           "Maybe it's a Wrong solution: f("
  30.           (vl-princ-to-string z)
  31.           ") = "
  32.           (VL-PRINC-TO-STRING f)
  33.         )
  34.       )
  35.     )
  36.   )
  37.   (princ)
  38. )
  39.  
  40. (defun C:t3 (/ a b c d e S z z1 z2 z3 z4)
  41.   (initget 1)
  42.   (setq a (getreal "\nCoefficient a:"))
  43.   (initget 1)
  44.   (setq b (getreal "\nCoefficient b:"))
  45.   (initget 1)
  46.   (setq c (getreal "\nCoefficient c:"))
  47.   (initget 1)
  48.   (setq d (getreal "\nCoefficient d:"))
  49.   (UTI:BENCH
  50.     20000 
  51.     (list
  52.       (list 'Math:Cubic_Equation a b c d)
  53.       (list 'Math:Cubic_Equation1 a b c d)
  54.     )
  55.   )
  56.   (foreach n '(Math:Cubic_Equation Math:Cubic_Equation1)
  57.     (setq S (apply n (list a b c d)))            		;get the solutions
  58.     (foreach z S
  59.       (princ "\n")
  60.       (princ (mapcar '(lambda (x) (rtos x 2 20)) z))      	;print them.
  61.     )
  62.     (foreach z S                                          	;Check every solution
  63.       (setq f (MATH:CReal_Polynomial z (list a b c d)))
  64.       (if (not (equal f '(0 0) 1e-6))                     	;if f(z) /= 0.0,maybe it's caused by floating point operation;
  65.         (alert
  66.           (strcat
  67. 	    "Use "
  68. 	    (VL-PRINC-TO-STRING n)
  69.             ",Maybe it's a Wrong solution: f("
  70.             (vl-princ-to-string z)
  71.             ") = "
  72.             (VL-PRINC-TO-STRING f)
  73.           )
  74.         )
  75.       )
  76.     )
  77.   )
  78.   (princ)
  79. )
  80. (princ "\nThe command is : test.")
  81. (vl-load-com)
  82. (princ)

发表评论

您的电子邮箱地址不会被公开。 必填项已用*标注