Author Topic: Solve Multiple Equation (Read 5971 times)

chlh_jd · « **on:** August 31, 2012, 02:26:55 AM »

We know a cubic or quartic equation can be solved by root formula , how to solve higher power equation by vlisp ?
Any suggestions are welcome .

irneb · « **Reply #1 on:** August 31, 2012, 03:20:29 AM »

You mean something like this: http://rosettacode.org/wiki/Roots_of_a_quadratic_function

The basic theory is using Lagrange resolvents. Note the hint for higher order functions at the end of that section, basically turning this into recursion.

chlh_jd · « **Reply #2 on:** August 31, 2012, 06:10:33 AM »

Quote from: irneb on August 31, 2012, 03:20:29 AM

You mean something like this: http://rosettacode.org/wiki/Roots_of_a_quadratic_function
The basic theory is using Lagrange resolvents. Note the hint for higher order functions at the end of that section, basically turning this into recursion.

Thank you very much .

chlh_jd · « **Reply #3 on:** September 01, 2012, 01:13:46 AM »

solve cubic and quartic equation ,
Need your help

for test or optimization .
For cubic equation

Code - Auto/Visual Lisp: [Select]

(defun ->eq^2  (a b c / d)
  ;; Solve quadratic equation : ax^2+bx+c = 0
  ;; Vieta theorem
  ;; ax^2+bx+c=0
  (setq a (float a)
        b (float b)
        c (float c))
  (if (/= a 0)
    (cond ((= (setq d (- (* b b) (* 4 a c))) 0)
           (list (/ b -2. a)))
          ((> d 0)
           (list (/ (+ (- b) (sqrt d)) 2. a)
                 (/ (- (- b) (sqrt d)) 2. a))))
    (if (/= b 0)
      (list (/ (- c) b)))))
;;
;;(->eq^3 8 -18908 14648440 -3730814500)-> (987.5 655. 721.)
(defun ->eq^3  (a b c d / aa bb cc dd y1 y2 x1 x2 x3 m)
  ;; solve cubic equation : ax^3+bx^2+cx+d = 0
  ;; &#30427;&#37329;&#20844;&#24335;&#27714;&#35299;&#27861; Though Fan ShengJin Method . (Fan ShengJin a Chinese mathematician)
  ;; code by GSLS(SS) 2012-9-1
  ;; Here only return Real solutions .
  ;;
  (setq a (float a)
        b (float b)
        c (float c)
        d (float d))
  (if (/= a 0)
    (progn
      (setq aa (- (* b b) (* 3 a c))
            bb (- (* b c) (* 9 a d))
            cc (- (* c c) (* 3 b d))
            dd (- (* bb bb) (* 4 aa cc)))
      (cond ((= aa bb 0)
             (setq x1 (/ b -3 a))
             (if (/= b 0)
               (setq x2 (/ c b -1)))
             (if (/= c 0)
               (setq x3 (/ (* 3 d) c))
               )
             (vl-remove nil (list x1 x2 x3))
             )
            ((> dd 0)
             (setq Y1 (+ (* aa b) (* 1.5 a (+ (- bb) (sqrt dd))))
                   Y2 (+ (* aa b) (* 1.5 a (- (- bb) (sqrt dd)))))
             (setq X1 (/ (+ b (expt3 Y1) (expt3 Y2)) -3. a))
             (list X1))
            ((and (= dd 0) (/= aa 0))
             (setq X1 (- (/ bb aa) (/ b a))
                   X2 (/ bb aa -2))
             (list X1 X2))
            ((and (< dd 0) (> aa 0.))
             (setq m
                    (acos
                      (/ (- (* 2. aa b) (* 3. a bb))
                         2.
                         (sqrt(expt aa 3)))));_m &#8712; [0 , pi]
             (setq x1 (/ (+ (- b) (* -2. (sqrt aa) (cos (/ m 3.))))
                         3.
                         a)                
                   x2 (/ (+ (- b) (* (sqrt aa) (+ (cos (/ m 3.)) (* (sqrt 3) (sin (/ m 3.)))))
                            ) 
                         3.
                         a)
                   x3 (/ (+ (- b) (* (sqrt aa) (- (cos (/ m 3.)) (* (sqrt 3) (sin (/ m 3.)))))
                            ) 
                         3.
                         a)
                   )
             (list x1 x2 x3))
            )
      )
    (->eq^2 b c d)
    )
  )
;; eq (expt num (/ 1. 3.))
(defun expt3 (num)
    (if (>= num 0)
      (expt num (/ 1. 3.))
      (/ (* num (expt (abs num) (/ 1. 3.))) (abs num))))
;;arccos(num)
;;
(defun acos (a / b)
  (if (and (numberp a)
           (<= (abs a) 1.0)
      )
    (cond ((= a 0.0)
           (* pi 0.5)
          )       
          (t
           (setq b (atan (/ (sqrt (- 1 (* a a)))
                    a
                 )
           ))
           (if (>= b 0)
             b
             (+ b pi))
          )
    )
  )
)
 

For quartic equation .

Code - Auto/Visual Lisp: [Select]

;;(->eq^4 1 118 4727 76538 426816)-> (-54.0 -32.0 -19.0 -13.0)
(defun ->eq^4    (a b c d e / yl sm m n xl res)
  ;; solve quartic equation
  ;; Here only return Real solutions .
  ;; code by GSLS(SS) 2012-9-1
  (if (and (/= a 0)
           (setq b  (/ (float b) a)
                 c  (/ (float c) a)
                 d  (/ (float d) a)
                 e  (/ (float e) a)
                 a  (/ (float a) a)
                 yl (->eq^3
                      8.
                      (* -4. c)
                      (- (* 2. b d) (* 8. e))
                      (- (* e (- (* 4. c) (* b b))) (* d d)))
                 ))
    (foreach y  yl
      (setq sm (+ (* 8 y) (* b b) (* -4 c)))
      (if (> sm 0)
        (progn
          (setq m (sqrt sm)
                n (- (* b y) d))
          (setq xl
                 (append (->eq^2 2. (+ b m) (* 2. (+ y (/ n m))))
                         (->eq^2 2. (- b m) (* 2. (- y (/ n m))))))
          (foreach x  xl
            (if (not (eqmember x res 1e-9))
              (setq res (cons x res)))))))
    (if (= a 0)
      (setq res (->eq^3 b c d e)))
    )
  (vl-sort res (function <))
  )
;;
;;;eqmember
(defun eqmember (p l eps)
  (cond ((not l)  nil)
          ((equal p (car l) eps) l)
          (t (eqmember p (cdr l) eps))
          ))
 

kdub:edit code=cadlisp-7

Stefan · « **Reply #4 on:** September 01, 2012, 06:42:52 PM »

Using Newton's Method, here is a general solver.
Good results so far, but still need rigorous testing.

Code - Auto/Visual Lisp: [Select]

;Polynomial equations solver 
;input: a list containing equation coefficients
;ex:  x^4-5x^3+5x^2+5x-6=0
;    (ph:solver '(1 -5 5 5 -6))
;returns unique real solutions
;Stefan M. 01.09.2012 - v. 1.01
 
;;;(defun test (n / l)
;;;  (repeat n (setq l (cons (- (rem (dos_random) 20) 10) l)))
;;;  (foreach x (ph:solver l)
;;;    (print (rtos (val l x) 2 16))
;;;    )
;;;  (princ)
;;;  )
 
(defun ph:solver (l / lst r)
  (while (zerop (car l)) (setq l (cdr l)))
  (setq l (unit (mapcar 'float l)))
  (while (> (length l) 2)
    (setq lst (cons l lst)
          l (unit (drv l))
          )
    )
  (setq r (list (- (/ (cadr l) (car l)))))
  (while lst
    (setq r (solve (car lst) (reverse r))
          lst (cdr lst))
    )
  r
  )
 
(defun solve (l xd / newton d s a b lst r y y1)
  (defun newton (x)
    (if (equal (val l x) 0.0 1e-9) x
      (newton (- x (/ (val l x) (val d x))))
      )
    )
  (or xd (setq xd '(0.0)))
  (setq d (drv l)
        a (car xd)
        y (val l a)
        )
  (cond
    ((minusp (setq s (* (sign (car l)) (sign y))))
     (setq lst (list (+ a (* 0.1 (abs a)))))
    )
    ((equal y 0.0 1e-9)
     (setq lst (list a))
    )
  )
  (while (cdr xd)
    (setq a (car xd) y (val l a) b (cadr xd) y1 (val l b))
    (cond ((minusp (* (sign y) (sign y1)))
           (setq lst (cons (/ (+ a b) 2.0) lst))
          )
          ((equal y1 0.0 1e-9)
           (setq lst (cons b lst))
          )
    )
    (setq xd (cdr xd))
    )
  (setq a (car xd) y (val l a))
  (cond
    ((minusp (- (* (sign (car l)) (expt -1.0 (length l)) (sign y))))
     (setq lst (cons (- a (* 0.1 (abs a))) lst)))
    )
  (mapcar 'newton lst)
  )
 
(defun grd (l / i r)
  (setq i -1)
  (repeat (length l) (setq r (cons (setq i (1+ i)) r)))
)
 
(defun val (l x)
  (if (zerop x)
    (float (last l))
    (apply '+ (mapcar (function (lambda (a b) (* a 1. (expt x b)))) l (grd l)))
  )
)
 
(defun drv (l)
  (mapcar '* l (reverse (cdr (reverse (grd l)))))
)
 
(defun sign (x)
  (if (minusp x) -1.0
    (if (zerop x) 0.0 1.0)
  )
)
 
(defun unit (l)
  (setq d (sqrt (apply '+ (mapcar '* l l))))
  (mapcar '(lambda (x) (/ x d)) l)
  )
 

kdub:edit code=cadlisp-7

edit: code updated by Stefan

chlh_jd · « **Reply #5 on:** September 02, 2012, 03:25:15 AM »

Thanks Kerry for Editing the codes

Quote from: Stefan on September 01, 2012, 06:42:52 PM

Using Newton's Method, here is a general solver.
Good results so far, but still need rigorous testing.
Code - Auto/Visual Lisp: [Select]
;Polynomial equations solver
;input: a list containing equation coefficients
;ex: x^4-5x^3+5x^2+5x-6=0
; (ph:solver '(1 -5 5 5 -6))
;returns unique real solutions
;Stefan M. 01.09.2012
(defun ph:solver (l)
(while (zerop (car l)) (setq l (cdr l)))
......

Code - Auto/Visual Lisp: [Select]
_$ (ph:solver '(8 -18908 14648440 -3730814500))
(655.0 721.0 987.5)
_$ (ph:solver '(1 118 4727 76538 426816))
(-54.0 -32.0 -19.0 -13.0)

Nice Routine , Stefan .

chlh_jd · « **Reply #6 on:** September 02, 2012, 07:05:47 AM »

The Newton mode efficiency is relatively low, Is it can be improved ?

Code: [Select]

(defun f1 (l)
  (->eq^4 (car l) (cadr l) (caddr l) (cadddr l) (last l)))
(setq l (list 1 -5 5 5 -6))
(quickbench (mapcar (function (lambda (f) (list f 'l))) (list f1 ph:solver)))

Test result

Code: [Select]

_$ 
Benchmarking .. done for 8192 iterations. Sorted from fastest.
Statement                                Increment  Time(ms) Normalize  Relative
--------------------------------------------------------------------------------
(#<USUBR @000000002cd081d8 F1> L)             8192      1653      1653     24.14
(#<USUBR @0000000032caacc8 PH:SOLVER> L)       256      1247     39904      1.00
--------------------------------------------------------------------------------
 
_$

Stefan · « **Reply #7 on:** September 02, 2012, 08:32:43 AM »

Quote from: chlh_jd on September 02, 2012, 07:05:47 AM

The Newton mode efficiency is relatively low, Is it can be improved ?

Of course it is slower because it uses an iterative calculation. It is also inaccurate and it fails most of the times

However, the differences in your benchmark test doesn't comes from newton's method but rather from my algorithm.
Newton's method needs an initial guess in order to find a root, so I tried to isolate roots by solving f'(x)=0 (between 2 consecutive roots of first derivate can be at most 1 root of function). So my algorithm solves n-1 equations and this is the cause of it's slowness. Of course, it is useful for higher degree polynomial equations, where you can not find an exact solution.

Updated the code in my previous post... It is not faster (au contraire) but it didn't fail so often

In some cases, that initial guess is not good so newton's method is not convergent. I'll try to solve this, but my initial guess is that I can not make it "full proof".

chlh_jd · « **Reply #8 on:** September 02, 2012, 08:54:01 AM »

Hi , Strfan . Thank you for your wonderful interpretation ! I am looking forward your improvements very much

Is this way more efficient :
Though cal Eigendecomposition of a matrix http://en.wikipedia.org/wiki/Eigendecomposition_of_a_matrix

Code: [Select]

;;High-order Equation .
;;a1x^n+a2x^(n-1)+...+akx^(n+1-k)+...+an = 0
;;Create the matrix 
(( -a1 1 0 ... 0)
 (  -a2 0 1 ... 0)
   ...
 ( - a(n-1) ... 1)
 ( -an  0  ...   0 ))
;;Cal Matrix Eigenvalue

chlh_jd · « **Reply #9 on:** September 02, 2012, 08:25:21 PM »

Hi , Strfan . your new's is faster .
I change some codes , I hope you don't mind .

Code: [Select]

;;Polynomial equations solver 
;;input: a list containing equation coefficients
;;ex:  x^4-5x^3+5x^2+5x-6=0
;;    (ph:solver '(1 -5 5 5 -6))
;;returns unique real solutions
;;Stefan M. 01.09.2012 - v. 1.01

;;;(defun test (n / l)
;;;  (repeat n (setq l (cons (- (rem (dos_random) 20) 10) l)))
;;;  (foreach x (ph:solver l)
;;;    (print (rtos (val l x) 2 16))
;;;    )
;;;  (princ)
;;;  )

(defun ph:solver  (l / lst r0 r)
  (while (zerop (car l)) (setq l (cdr l)))
  ;;add last item = 0 case 
  (while (zerop (last l))
    (setq r0 (cons 0.0 r0)
	  l  (butlast l)))
  ;;
  (setq l (unit (mapcar 'float l)))
  (while (> (length l) 2)
    (setq lst (cons l lst)
	  l   (unit (drv l))
	  )
    )
  (setq r (list (- (/ (cadr l) (car l)))))
  (while lst
    (setq r   (solve (car lst) (reverse r))
	  lst (cdr lst))
    )
  (append r0 r)
  )
(defun solve  (l xd / d a b lst r y y1)
  ;; l ---- Equation coefficients List   
  (or xd (setq xd '(0.0)))
  (setq	d (drv l)
	a (car xd)
	y (val l a)
	)
  (if (equal y 0.0 1e-9)
    (setq lst (list a))
    (if	(xminusp (list (car l) y))
      (setq lst (list (+ a (* 0.1 (abs a)))))))
  (while (cdr xd)
    (setq b  (cadr xd)
	  y1 (val l b))
    (if	(equal y1 0.0 1e-9)
      (setq lst (cons b lst))
      (if (xminusp (list y y1))
	(setq lst (cons (/ (+ a b) 2.0) lst))))
    (setq xd (cdr xd)
	  a  b
	  y  y1))
					;(setq a (car xd) y (val l a))
  (if (xplusp (list (car l) (expt -1.0 (rem (length l) 2)) y))
    (setq lst (cons (- a (* 0.1 (abs a))) lst)))
  (mapcar (function (lambda (x) (newton x 1. 1e-9))) lst)
  )
;;add Acceleration factor  
(defun newton  (x af eps / y)
  ;; x ---- Initial value    
  ;; af ---- Acceleration factor , Such as 0.9
  ;; eps ---- Solution accuracy
  (while (not (equal (setq y (val l x)) 0.0 eps))
    (setq x (- x (* af (/ y (val d x)))))
    ))
 ;|
;;(grd '(1 -5 5 5 -6))
(defun grd  (l / n)
  (setq n (length l))
  (mapcar (function (lambda (x)
		      (setq n (1- n))))
	  l)
  )
|;
(defun val  (l x)
  (if (zerop x)
    (last l)
    ((lambda (n)
       (apply
	 (function +)
	 (mapcar
	   (function (lambda (a) (* a 1. (expt x (setq n (1- n))))))
	   l))
       )
      (length l)
      )
    )
  )

(defun drv  (l)
  ((lambda (n)
     (butlast
       (mapcar (function (lambda (a) (* a (setq n (1- n))))) l)))
    (length l)))

(defun unit  (l)
  (setq d (sqrt (apply '+ (mapcar '* l l))))
  (mapcar '(lambda (x) (/ x d)) l)
  )

(defun butlast	(l)
  (reverse (cdr (reverse l))))

(defun xminusp	(l)
  (and (not (vl-some (function zerop) l))
       (= (rem (length (vl-remove-if-not (function minusp) l)) 2)
	  1)))

(defun xplusp  (l)
  (and (not (vl-some (function zerop) l))
       (zerop
	 (rem (length (vl-remove-if-not (function minusp) l)) 2))))

Test result .

Code: [Select]

_$ 
Benchmarking .. done for 8192 iterations. Sorted from fastest.
Statement                                Increment  Time(ms) Normalize  Relative
--------------------------------------------------------------------------------
(#<USUBR @000000002cc864d0 F1> L)             8192      1716      1716     11.06
(#<USUBR @000000002c6a1930 PH:SOLVER> L)       512      1186     18976      1.00
-------------------------------------------------------------------------------- 
_$

Gasty · « **Reply #10 on:** September 03, 2012, 02:07:19 AM »

Hi,

This is another way using Laguerre's method:

Code: [Select]

;;2012 Gaston Nunez
;;Laguerre's Method
;;;;;Params;;;
;;Px:=Polynomy in list form, example: 4x^3 + 3x^2 + 2x + 1:=(4 3 2 1)
;;m :=Iterations, keep small usually 5-10 is enought
;;guess:= Initial guessed root
;;Use example: (getroot '(4 3 2 1) 5 7.2)



(defun fd(px / Cn n fdd)
  (setq n (- (length px) 1))
  (setq Cn (getCn (+ n 1)))
  (setq fdd (product Cn px))
  (append '(0) (reverse (cdr(reverse fdd))))
)


(defun getCn(n / i)
  (setq i 0)
  (repeat n
    (setq Cn (cons i Cn))
    (setq i (+ i 1))
  )
Cn)

(defun product(a b)
  (mapcar '(lambda(x y)(* x y)) a b)
)

(defun evaluate(px x / mono)
  (setq mono (getmonoeval x (length px)))
  (apply '+ (product px mono))
)

(defun getmonoeval(x n / mono i)
 (setq i 0) 
 (repeat n
   (setq mono (cons (expt x i) mono))
   (setq i (+ i 1))
 )
mono  
)

(defun Gk(px ffd x0)
  (/ (evaluate ffd x0)(evaluate px x0))
)

(defun Hk(px ffd ssd x0 / hhx)
  (setq hhx (/ (evaluate ssd x0) (evaluate px x0)))
  (- (* (Gk px ffd x0)(Gk px ffd x0)) hhx)
)  

(defun Disc(px ffd ssd x0 / n nHk A1)
  (setq n (- (length px) 1))
  (setq nHk (* n (Hk px ffd ssd x0)))
  (setq A1 (- nHk (* (Gk px ffd x0)(Gk px ffd x0))))
  (setq A1 (* (- n 1) A1))
  (sqrt (abs A1))
)

(defun Ak(px ffd ssd x0 n / F1 F2)
  (setq F1 (+ (Gk px ffd x0) (Disc px ffd ssd x0)))
  (setq F2 (- (Gk px ffd x0) (Disc px ffd ssd x0)))
  (if (> (abs F1) (abs F2))
    (/ n F1)
    (/ n F2)
  )
)


(defun GetRoot(px m guess / n)
  (setq n (- (length px) 1))
  (setq ffd (fd px))
  (setq ssd (fd ffd))
  (setq Xk guess)
  (repeat m
    (setq Xk (- Xk (Ak px ffd ssd Xk n)))
  )
 (print Xk) 
)

chlh_jd · « **Reply #11 on:** September 03, 2012, 05:12:58 AM »

Quote from: Gasty on September 03, 2012, 02:07:19 AM

Hi,

This is another way using Laguerre's method:

Code: [Select]
;;2012 Gaston Nunez ;;Laguerre's Method ;;;;;Params;;; ;;Px:=Polynomy in list form, example: 4x^3 + 3x^2 + 2x + 1:=(4 3 2 1) ;;m :=Iterations, keep small usually 5-10 is enought ;;guess:= Initial guessed root ;;Use example: (getroot '(4 3 2 1) 5 7.2) (defun fd(px / Cn n fdd) ...

Thank you for give a new method .
It dosen't work for bellow code ?

Code: [Select]

(grtroot (list 1 -5 5 5 -6) 5 0.0)

News:

Author Topic: Solve Multiple Equation (Read 5971 times)

chlh_jd

Solve Multiple Equation

irneb

Re: Solve Multiple Equation

chlh_jd

Re: Solve Multiple Equation

chlh_jd

Re: Solve Multiple Equation

Stefan

Re: Solve Multiple Equation

chlh_jd

Re: Solve Multiple Equation

chlh_jd

Re: Solve Multiple Equation

Stefan

Re: Solve Multiple Equation

chlh_jd

Re: Solve Multiple Equation

chlh_jd

Re: Solve Multiple Equation

Gasty

Re: Solve Multiple Equation

chlh_jd

Re: Solve Multiple Equation