¿by numerical methods?
Bisector, newton...?
bisector method is very simple to find one solution (initial value near solution value)
;;;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))
)
)
)
)
)
;;;----------------------------------------------------; ;
;;;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))
)
)
)
)
)
)
)
HI,Ribarm!(defun MATH:CubeRoot (x)
(if (< x 0)
(- (expt (- x) 0.33333333333333333333333))
(expt x 0.33333333333333333333333)
)
)
There is also a sub-function
(defun decomposemat ( mat ) ; mat - nxn matrix
;; mat = ((1.0 0.0 0.0) (0.0 2.0 0.0) (0.0 0.0 3.0))
;; eigenvalues = 1.0; 2.0; 3.0
;; [E]det = 0; [E]det = (1-x)(2-x)(3-x) = 0
;; (2-x-2x+x^2)*(3-x) = (2-3x+x^2)*(3-x) = 6-2x-9x+3x^2+3x^2-x^3 = -x^3+6x^2-11x+6 = 0
;; (x-1)(x-2)(x-3) = 0 ;; (x^2-2x-x+2)(x-3) = 0 ;; x^3-3x^2-3x^2+9x+2x-6 = 0 ;; x^3-6x^2+11x-6 = 0 ;; -1*[E]det = 0
;; a=1; b=-6; c=11; d=-6 ;; (1.0 -6.0 11.0 -6.0)
;;
;; mat = ((2.0 3.0 8.0) (1.0 0.0 5.0) (-1.0 5.0 0.0))
;; eigenvalues = -4.19258; 1.19258; 5.0
;; (x+4.19258)(x-1.19258)(x-5.0) = 0 ;; (x^2-1.19258x+4.19258x-5.0)(x-5.0) = 0 ;; (x^2+3x-5)(x-5) = 0 ;; x^3-5x^2+3x^2-15x-5x+25 = 0 ;; x^3-2x^2-20x+25 = 0
;; a=1; b=-2; c=-20; d=25 ;; (1.0 -2.0 -20.0 25.0)
;;
;; mat = ((1.0 2.0) (2.0 1.0))
;; eigenvalues = -1.0; 3.0
;; (x+1)(x-3) = 0 ;; x^2-3x+x-3 = 0 ;; x^2-2x-3 = 0
;; a=1; b=-2; c=-3 ;; (1.0 -2.0 -3.0)
;;
;; mat = ((1.0 3.0 2.0 0.5) (-5.0 7.0 -8.0 6.0) (9.0 12.0 -10.0 11.0) (13.0 -15.0 14.0 -16.0))
;; eigenvalues = -18.4868; -5.79211
;;
;; mat - 4x4 sometimes has also complex numbers eigenvalues
;;
;; mat = ((0.5 0.25 1.5 -2.7 -4.2 3.0 2.4) (1.5 1.25 2.5 -5.7 -1.2 1.0 3.4) (2.5 2.25 0.5 -1.7 -3.2 2.0 4.4) (1.5 3.25 3.5 -4.7 -2.2 4.0 0.4) (3.5 4.25 4.5 -0.7 -0.2 0.0 1.4) (4.5 2.25 1.5 -3.7 -2.2 1.0 3.4) (5.5 5.25 5.5 -0.7 -5.2 5.0 5.4))
;; eigenvalues = -2.14872; 3.80431; 9.59801
;;
;; mat - 7x7 sometimes has also complex numbers eigenvalues
;;
;;
;; mat = ((4.0 2.0 6.0 8.0 2.0) (2.0 10.0 9.0 13.0 7.0) (6.0 9.0 14.0 20.0 12.0) (8.0 13.0 20.0 54.0 35.0) (2.0 7.0 12.0 35.0 43.0))
;; eigenvalues = 0.0660769; 4.18334; 7.50996; 18.2577; 94.9829
;;
;; this mat - 5x5 has only real eigenvalues
;; (x-0.0660769)(x-4.18334)(x-7.50996)(x-18.2577)(x-94.9829) = 0
;; (x^2-4.2494169x+0.276422138)(x^2-25.76766x+137.1145967)(x-94.9829) = 0
;; (x^4-25.76766x^3+137.1145967x^2-4.2494169x^3+109.4975299x^2-582.6570845x+0.276422138x^2-7.122751668x+37.90150997)(x-94.9829) = 0
;; (x^4-30.0170769x^3+246.8885487x^2-589.7798362x+37.90150997)(x-94.9829) = 0
;; x^5-30.0170769x^4+246.8885487x^3-589.7798362x^2+37.90150997x-94.9829x^4+2851.109013x^3-23450.19033x^2+56018.9992x-3599.995331 = 0
;; x^5-125x^4+3098x^3-24040x^2+56056.9x-3600 = 0
;; a=1; b=-125; c=3098; d=-24040; e=56056.9; f=-3600 ;; (1.0 -125.0 3098.0 -24040.0 56056.9 -3600.0)
;;
;; (car rtn) = 1
;; (cadr rtn) = (- (apply '+ eigenvalues))
;; (nth 2 rtn) = (+ (* [E1] [E2]) (* [E1] [E3]) (* [E2] [E3])) - for 3x3;
;; (+ (* [E1] [E2]) (* [E1] [E3]) (* [E1] [E4]) (* [E1] [E5]) (* [E2] [E3]) (* [E2] [E4]) (* [E2] [E5]) (* [E3] [E4]) (* [E3] [E5]) (* [E4] [E5])) - for 5x5
;; (nth 3 rtn) = (- (+ (* [E1] [E2] [E3]) (* [E1] [E2] [E4]) (* [E1] [E2] [E5]) (* [E1] [E3] [E4]) (* [E1] [E3] [E5]) (* [E1] [E4] [E5]) (* [E2] [E3] [E4]) (* [E2] [E3] [E5]) (* [E2] [E4] [E5]) (* [E3] [E4] [E5]))) - for 5x5
;; (nth 4 rtn) = (+ (* [E1] [E2] [E3] [E4]) (* [E1] [E2] [E3] [E5]) (* [E1] [E2] [E4] [E5]) (* [E1] [E3] [E4] [E5]) (* [E2] [E3] [E4] [E5])) - for 5x5
;; ... = ?? - for nxn
;; (last rtn) = (* (* (- (rem (1+ (length mat)) 2) 0.5) 2.0) (apply '* eigenvalues))