Author Topic: solving cubic equation problem...  (Read 3118 times)

0 Members and 1 Guest are viewing this topic.

ribarm

  • Gator
  • Posts: 3225
  • Marko Ribar, architect
solving cubic equation problem...
« on: December 05, 2019, 07:07:23 AM »
Hi all, I am having problem in solving cubic equation ax^3+bx^2+cx+d=0... No matter how deep I looked into Cardano's method, I just can't get correct results in AutoLisp although it looks everything good... Please if someone can see where is the problem let me know...

Code - Auto/Visual Lisp: [Select]
  1. (defun eigenvalues ( mat / eigenvalues-qubiccoefficients a b c d p q t1 t2 t3 x1 x2 x3 ) ; mat = 3x3 matrix
  2.  
  3.   (defun eigenvalues-qubiccoefficients ( mat / a11 a12 a13 a21 a22 a23 a31 a32 a33 a b c d ) ; mat = 3x3 matrix
  4.     (setq a11 (nth 0 (nth 0 mat)))
  5.     (setq a12 (nth 1 (nth 0 mat)))
  6.     (setq a13 (nth 2 (nth 0 mat)))
  7.     (setq a21 (nth 0 (nth 1 mat)))
  8.     (setq a22 (nth 1 (nth 1 mat)))
  9.     (setq a23 (nth 2 (nth 1 mat)))
  10.     (setq a31 (nth 0 (nth 2 mat)))
  11.     (setq a32 (nth 1 (nth 2 mat)))
  12.     (setq a33 (nth 2 (nth 2 mat)))
  13.     ;|
  14.     (setq det (- (+ (* (- a11 x) (- a22 x) (- a33 x)) (* a12 a23 a31) (* a13 a21 a32)) (+ (* a31 (- a22 x) a13) (* a32 a23 (- a11 x)) (* (- a33 x) a21 a12))))
  15.     (setq det (- (+ (* (+ (* a11 a22) (* (- a11) x) (* (- a22) x) (* x x)) (- a33 x)) (* a12 a23 a31) (* a13 a21 a32)) (+ (- (* a31 a22 a13) (* a31 a13 x)) (- (* a32 a23 a11) (* a32 a23 x)) (- (* a33 a21 a12) (* a21 a12 x)))))
  16.     (setq det (- (+ (- (+ (* a11 a22 a33) (* (- a11) a33 x) (* (- a22) a33 x) (* a33 x x)) (+ (* a11 a22 x) (* (- a11) x x) (* (- a22) x x) (* x x x))) (* a12 a23 a31) (* a13 a21 a32)) (- (+ (* a31 a22 a13) (* a32 a23 a11) (* a33 a21 a12)) (+ (* a31 a13 x) (* a32 a23 x) (* a21 a12 x)))))
  17.     (setq det (- (+ (* a11 a22 a33) (* a12 a23 a31) (* a13 a21 a32) (* x (+ (* (- a11) a33) (* (- a22) a33) (* (- a11) a22))) (* x x (+ a33 a11 a22)) (- (* x x x))) (- (+ (* a31 a22 a13) (* a32 a23 a11) (* a33 a21 a12)) (* x (+ (* a31 a13) (* a32 a23) (* a21 a12))))))
  18.     (setq det (+ (- (+ (* a11 a22 a33) (* a12 a23 a31) (* a13 a21 a32)) (+ (* a31 a22 a13) (* a32 a23 a11) (* a33 a21 a12))) (* x (+ (* (- a11) a33) (* (- a22) a33) (* (- a11) a22) (* a31 a13) (* a32 a23) (* a21 a12))) (* x x (+ a33 a11 a22)) (* x x x (- 1))))
  19.     |;
  20.     (setq a -1.0)
  21.     (setq b (+ a11 a22 a33))
  22.     (setq c (- (+ (* a31 a13) (* a32 a23) (* a21 a12)) (+ (* a11 a33) (* a22 a33) (* a11 a22))))
  23.     (setq d (- (+ (* a11 a22 a33) (* a12 a23 a31) (* a13 a21 a32)) (+ (* a31 a22 a13) (* a32 a23 a11) (* a33 a21 a12))))
  24.     (list a b c d)
  25.   )
  26.  
  27.   (mapcar 'set '(a b c d) (eigenvalues-qubiccoefficients mat))
  28.   (setq p (+ (* 9.0 a c) (* -3.0 (expt b 2.0))))
  29.   (setq q (+ (* 2.0 (expt b 3.0)) (* 27.0 (expt a 2.0) d) (* -9.0 a b c)))
  30.   (if (>= (+ (* 729.0 (expt q 2.0)) (* 108.0 (expt p 3.0))) 0.0)
  31.     (setq t1 (+ (if (>= (/ (+ (* -27.0 q) (sqrt (+ (* 729.0 (expt q 2.0)) (* 108.0 (expt p 3.0))))) 54.0) 0.0) (expt (/ (+ (* -27.0 q) (sqrt (+ (* 729.0 (expt q 2.0)) (* 108.0 (expt p 3.0))))) 54.0) (/ 1.0 3.0)) (- (expt (abs (/ (+ (* -27.0 q) (sqrt (+ (* 729.0 (expt q 2.0)) (* 108.0 (expt p 3.0))))) 54.0)) (/ 1.0 3.0)))) (if (>= (/ (+ (* -27.0 q) (sqrt (+ (* 729.0 (expt q 2.0)) (* 108.0 (expt p 3.0))))) 54.0) 0.0) (expt (/ (+ (* -27.0 q) (sqrt (+ (* 729.0 (expt q 2.0)) (* 108.0 (expt p 3.0))))) 54.0) (/ 1.0 3.0)) (- (expt (abs (/ (+ (* -27.0 q) (sqrt (+ (* 729.0 (expt q 2.0)) (* 108.0 (expt p 3.0))))) 54.0)) (/ 1.0 3.0))))))
  32.   )
  33.   (if (>= (+ (* 729.0 (expt q 2.0)) (* 108.0 (expt p 3.0))) 0.0)
  34.     (setq t2 (+ (if (>= (/ (+ (* -27.0 q) (sqrt (+ (* 729.0 (expt q 2.0)) (* 108.0 (expt p 3.0))))) 54.0) 0.0) (expt (/ (+ (* -27.0 q) (sqrt (+ (* 729.0 (expt q 2.0)) (* 108.0 (expt p 3.0))))) 54.0) (/ 1.0 3.0)) (- (expt (abs (/ (+ (* -27.0 q) (sqrt (+ (* 729.0 (expt q 2.0)) (* 108.0 (expt p 3.0))))) 54.0)) (/ 1.0 3.0)))) (if (>= (/ (- (* -27.0 q) (sqrt (+ (* 729.0 (expt q 2.0)) (* 108.0 (expt p 3.0))))) 54.0) 0.0) (expt (/ (- (* -27.0 q) (sqrt (+ (* 729.0 (expt q 2.0)) (* 108.0 (expt p 3.0))))) 54.0) (/ 1.0 3.0)) (- (expt (abs (/ (- (* -27.0 q) (sqrt (+ (* 729.0 (expt q 2.0)) (* 108.0 (expt p 3.0))))) 54.0)) (/ 1.0 3.0))))))
  35.   )
  36.   (if (>= (+ (* 729.0 (expt q 2.0)) (* 108.0 (expt p 3.0))) 0.0)
  37.     (setq t3 (+ (if (>= (/ (- (* -27.0 q) (sqrt (+ (* 729.0 (expt q 2.0)) (* 108.0 (expt p 3.0))))) 54.0) 0.0) (expt (/ (- (* -27.0 q) (sqrt (+ (* 729.0 (expt q 2.0)) (* 108.0 (expt p 3.0))))) 54.0) (/ 1.0 3.0)) (- (expt (abs (/ (- (* -27.0 q) (sqrt (+ (* 729.0 (expt q 2.0)) (* 108.0 (expt p 3.0))))) 54.0)) (/ 1.0 3.0)))) (if (>= (/ (- (* -27.0 q) (sqrt (+ (* 729.0 (expt q 2.0)) (* 108.0 (expt p 3.0))))) 54.0) 0.0) (expt (/ (- (* -27.0 q) (sqrt (+ (* 729.0 (expt q 2.0)) (* 108.0 (expt p 3.0))))) 54.0) (/ 1.0 3.0)) (- (expt (abs (/ (- (* -27.0 q) (sqrt (+ (* 729.0 (expt q 2.0)) (* 108.0 (expt p 3.0))))) 54.0)) (/ 1.0 3.0))))))
  38.   )
  39.   (if t1
  40.     (setq x1 (/ (- t1 b) (* 3.0 a)))
  41.   )
  42.   (if t2
  43.     (setq x2 (/ (- t2 b) (* 3.0 a)))
  44.   )
  45.   (if t3
  46.     (setq x3 (/ (- t3 b) (* 3.0 a)))
  47.   )
  48.   (princ (list x1 x2 x3))
  49.   (princ)
  50. )
  51.  
  52. (setq a1 '((2.0 1.0 0.0) (1.0 3.0 1.0) (0.0 1.0 2.0)))
  53. (setq a2 '((7.0 3.0 -2.0) (3.0 4.0 1.0) (-2.0 -1.0 3.0)))
  54. (setq a3 '((4.0 1.0 0.0) (1.0 0.0 -1.0) (1.0 1.0 -4.0)))
  55. (setq a4 '((1.0 1.0 0.5) (1.0 1.0 0.25) (0.5 0.25 2.0)))
  56. (setq a5 '((2.0 3.0 8.0) (1.0 0.0 5.0) (-1.0 5.0 0.0)))
  57. (setq a6 '((-3.0 1.0 -1.0) (-7.0 5.0 -1.0) (-6.0 6.0 -2.0)))
  58. (setq a7 '((1.0 0.0 0.0) (0.0 2.0 0.0) (0.0 0.0 3.0)))
  59. ;|
  60. (eigenvalues a1)
  61. (eigenvalues a2)
  62. (eigenvalues a3)
  63. (eigenvalues a4)
  64. (eigenvalues a5)
  65. (eigenvalues a6)
  66. (eigenvalues a7)
  67. |;
  68.  

Eigenvalue coefficients of cubic equation are correct... The problem is when I try to solve eigenvalue - cubic equation...
Cardano's method is in attached *.html
Marko Ribar, d.i.a. (graduated engineer of architecture)

:)

M.R. on Youtube

dgpuertas

  • Newt
  • Posts: 80
Re: solving cubic equation problem...
« Reply #1 on: December 05, 2019, 09:46:24 AM »

¿by numerical methods?
Bisector, newton...?

bisector method is very simple to find one solution (initial value near solution value)

ribarm

  • Gator
  • Posts: 3225
  • Marko Ribar, architect
Re: solving cubic equation problem...
« Reply #2 on: December 05, 2019, 10:23:32 AM »

¿by numerical methods?
Bisector, newton...?

bisector method is very simple to find one solution (initial value near solution value)

Like I did something wrong...
Quadratic equation has solution based on numerical method...

I thought it should work as well as simple algebra math... I just don't see where I made mistake, everything looks good to me...
Marko Ribar, d.i.a. (graduated engineer of architecture)

:)

M.R. on Youtube

ribarm

  • Gator
  • Posts: 3225
  • Marko Ribar, architect
Re: solving cubic equation problem...
« Reply #3 on: December 09, 2019, 02:59:29 AM »
I've found some workaround in attachment of this post :
http://www.theswamp.org/index.php?topic=43453.msg486924#msg486924

But then again I dived a little deeper and I decided to share my mod. of chlh_jd's [E] sub function...
Now the new problem arises :
We know how to get eigenvalues from square matrix, but how to get coefficients a*x^n+b*x^(n-1)+c*x^(n-2)+...+zz*x+zza=0
from square matrix - I did it for only 3x3 matrix...
And one more question even more important - we know equation - polynomial : a*x^n+b*x^(n-1)+c*x^(n-2)+...+zz*x+zza=0, but how then to construct square matrix that chlh_jd used in his examples...

I'll post my mod in above linked topic... I hope you'll find it useful as for me QR-decomposition don't work as expected, and Jakobi and Power methods are also somewhat wrong on my PC... I decided to follow first method of chlh_jd and changed his (Eigenvalue) sub...
Marko Ribar, d.i.a. (graduated engineer of architecture)

:)

M.R. on Youtube

myloveflyer

  • Newt
  • Posts: 152
Re: solving cubic equation problem...
« Reply #4 on: December 09, 2019, 09:28:45 PM »
Code: [Select]
;;;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!
The above code is written by highflybir!
hope this helps you!
Never give up !

myloveflyer

  • Newt
  • Posts: 152
Re: solving cubic equation problem...
« Reply #5 on: December 09, 2019, 09:30:54 PM »
Code: [Select]
(defun MATH:CubeRoot (x)
  (if (< x 0)
    (- (expt (- x) 0.33333333333333333333333))
    (expt x 0.33333333333333333333333)
  )
)
There is also a sub-function
Never give up !

ribarm

  • Gator
  • Posts: 3225
  • Marko Ribar, architect
Re: solving cubic equation problem...
« Reply #6 on: December 10, 2019, 07:23:54 AM »
Actually, I've found more comprehensive version of Highflyingbird's lisp - COMPLEX.LSP... I'll attach it here, but thanks anyway for your input...
Marko Ribar, d.i.a. (graduated engineer of architecture)

:)

M.R. on Youtube

ribarm

  • Gator
  • Posts: 3225
  • Marko Ribar, architect
Re: solving cubic equation problem...
« Reply #7 on: December 12, 2019, 01:31:20 PM »
For those that are in the mood of solving mathematical problems...

Quote
(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))
« Last Edit: December 12, 2019, 05:49:50 PM by ribarm »
Marko Ribar, d.i.a. (graduated engineer of architecture)

:)

M.R. on Youtube

ribarm

  • Gator
  • Posts: 3225
  • Marko Ribar, architect
Re: solving cubic equation problem...
« Reply #8 on: December 13, 2019, 10:47:47 AM »
I think I solved this one...

Code - Auto/Visual Lisp: [Select]
  1. (defun [E]-[C] ( lst / str+ foo k r rr )
  2.  
  3. ;; Command: ([e]-[c] '(2.0 2.0 2.0))
  4. ;; (1.0 -6.0 12.0 -8.0)
  5. ;; Command: ([e]-[c] '(1.0 3.0 5.0))
  6. ;; (1.0 -9.0 23.0 -15.0)
  7. ;; Command: ([e]-[c] '(0.0660769 4.18334 7.50996 18.2577 94.9829))
  8. ;; (1.0 -125.0 3098.0 -24040.0 56056.9 -3600.0)
  9.  
  10.   (defun str+ ( source / prefix rchar )
  11.     (cond
  12.       ( (= source "") "A")
  13.       ( (= (setq prefix (substr source 1 (1- (strlen source)))
  14.                   rchar (substr source (strlen source))
  15.             )
  16.             "Z"
  17.         )
  18.         (strcat (str+ prefix) "A")
  19.       )
  20.       ( (strcat prefix (chr (1+ (ascii rchar)))))
  21.     )
  22.   )
  23.  
  24.   (defun foo ( l n / x xx s f zzzzz k a b )
  25.     (setq s "(defun f ( l / zzzz ) ")
  26.     (setq k 0)
  27.     (repeat n (setq k (1+ k)) (if (= k 1) (setq s (strcat s "(foreach " (progn (if (null x) (setq x (str+ "")) (setq x (str+ x))) x) " l ")) (setq s (strcat s "(foreach " (progn (if (null x) (setq x (str+ "")) (setq x (str+ x))) x) (if (= k 2) (strcat (setq a (strcat " (vl-remove " (progn (if (null xx) (setq xx (str+ "")) (setq xx (str+ xx))) xx))) " l " (setq b ")")) (strcat (setq a (strcat a " (vl-remove " (progn (if (null xx) (setq xx (str+ "")) (setq xx (str+ xx))) xx))) " l " (setq b (strcat b ")"))))))))
  28.     (setq s (strcat s "(setq zzzz (list "))
  29.     (setq xx nil)
  30.     (repeat n (setq s (strcat s (progn (if (null xx) (setq xx (str+ "")) (setq xx (str+ xx))) xx) " ")))
  31.     (setq s (strcat s "))"))
  32.     (setq s (strcat s "(if (not (vl-some '(lambda ( zzz ) (vl-every '(lambda ( yyy ) (vl-position yyy zzz)) zzzz)) zzzzz)) (setq zzzzz (cons zzzz zzzzz)))"))
  33.     (repeat n (setq s (strcat s ") ")))
  34.     (setq s (strcat s ")"))
  35.     (eval (read s))
  36.     (f l)
  37.     zzzzz
  38.   )
  39.  
  40.   (setq k 0)
  41.   (setq lst (mapcar '(lambda ( x ) (cons (setq k (1+ k)) x)) lst))
  42.   (setq k 0)
  43.   (setq rr (cons 1.0 rr))
  44.   (repeat (length lst)
  45.     (setq r (foo (mapcar 'car lst) (setq k (1+ k))))
  46.     (setq r ((if (= (rem k 2) 0) + -) (apply '+ (mapcar '(lambda ( x ) (if (= (length x) 1) (cdr (assoc (car x) lst)) (apply '* (mapcar '(lambda ( y ) (cdr (assoc y lst))) x)))) r))))
  47.     (setq rr (cons r rr))
  48.   )
  49.   (reverse rr)
  50. )
  51.  
« Last Edit: December 18, 2019, 08:30:59 AM by ribarm »
Marko Ribar, d.i.a. (graduated engineer of architecture)

:)

M.R. on Youtube

ribarm

  • Gator
  • Posts: 3225
  • Marko Ribar, architect
Re: solving cubic equation problem...
« Reply #9 on: December 18, 2019, 07:41:23 AM »
I've found Lee's sub for combinations without repetitions, so this my last code could be just simple and length of initial list can now converge to infinity :

Code - Auto/Visual Lisp: [Select]
  1. (defun [E]-[C] ( lst / nCr k r rr )
  2.  
  3. ;; Command: ([e]-[c] '(2.0 2.0 2.0))
  4. ;; (1.0 -6.0 12.0 -8.0)
  5. ;; Command: ([e]-[c] '(1.0 3.0 5.0))
  6. ;; (1.0 -9.0 23.0 -15.0)
  7. ;; Command: ([e]-[c] '(0.0660769 4.18334 7.50996 18.2577 94.9829))
  8. ;; (1.0 -125.0 3098.0 -24040.0 56056.9 -3600.0)
  9.  
  10.   (defun nCr ( l r )
  11.   ;https://www.cadtutor.net/forum/topic/59422-combinations-without-repetition-3-and-4-element/
  12.      (cond
  13.          (   (< r 2)
  14.              (mapcar 'list l)
  15.          )
  16.          (   l
  17.              (append
  18.                  (mapcar '(lambda ( x ) (cons (car l) x)) (nCr (cdr l) (1- r)))
  19.                  (nCr (cdr l) r)
  20.              )
  21.          )
  22.      )
  23.   )
  24.  
  25.   (setq k 0)
  26.   (setq lst (mapcar '(lambda ( x ) (cons (setq k (1+ k)) x)) lst))
  27.   (setq k 0)
  28.   (setq rr (cons 1.0 rr))
  29.   (repeat (length lst)
  30.     (setq r (nCr (mapcar 'car lst) (setq k (1+ k))))
  31.     (setq r ((if (= (rem k 2) 0) + -) (apply '+ (mapcar '(lambda ( x ) (if (= (length x) 1) (cdr (assoc (car x) lst)) (apply '* (mapcar '(lambda ( y ) (cdr (assoc y lst))) x)))) r))))
  32.     (setq rr (cons r rr))
  33.   )
  34.   (reverse rr)
  35. )
  36.  

Regards, and thanks Lee...
« Last Edit: December 18, 2019, 08:24:55 AM by ribarm »
Marko Ribar, d.i.a. (graduated engineer of architecture)

:)

M.R. on Youtube

ribarm

  • Gator
  • Posts: 3225
  • Marko Ribar, architect
Re: solving cubic equation problem...
« Reply #10 on: December 19, 2019, 05:41:19 AM »
For complex eigenvalues, coefficients of polynomia are also complex...

Code - Auto/Visual Lisp: [Select]
  1. (defun [E]-[C] ( lst / nCr cxc k r c rr )
  2.  
  3. ;; Command: ([e]-[c] '((2.0 0) (2.0 0) (2.0 0)))
  4. ;; ((1.0 0.0) (-6.0 0.0) (12.0 0.0) (-8.0 0.0))
  5. ;; Command: ([e]-[c] '((1.0 0) (3.0 0) (5.0 0)))
  6. ;; ((1.0 0.0) (-9.0 0.0) (23.0 0.0) (-15.0 0.0))
  7. ;; Command: ([e]-[c] '((0.0660769 0) (4.18334 0) (7.50996 0) (18.2577 0) (94.9829 0)))
  8. ;; ((1.0 0.0) (-125.0 0.0) (3098.0 0.0) (-24040.0 0.0) (56056.9 0.0) (-3600.0 0.0))
  9. ;;
  10. ;; Command: ([e]-[c] '((1.0 0.5) (3.0 -1.5) (5.0 2.0)))
  11. ;; ((1.0 0.0) (-9.0 -1.0) (25.75 3.0) (-18.75 -7.5))
  12.  
  13. ;; 0 = (1.0 + 0.0i)(1.0 + 0.5i)^3+(-9.0 - 1.0i)(1.0 + 0.5i)^2+(25.75 + 3.0i)(1.0 + 0.5i)-18.75-7.5i
  14. ;; 0 = (1.0 + 0.0i)(1.0^3+3*0.5i+3*(0.5i)^2+(0.5i)^3)+(-9.0 - 1.0i)(1.0 + 1.0i - 0.25)+25.75+12.875i+3i-1.5-18.75-7.5i
  15. ;; 0 = 1+1.5i-0.75-0.125i-6.75-9i-0.75i+1+25.75+12.875i+3i-1.5-18.75-7.5i
  16. ;; 0 = 0 + 0i = 0
  17.  
  18.   (defun nCr ( l r )
  19.   ;https://www.cadtutor.net/forum/topic/59422-combinations-without-repetition-3-and-4-element/
  20.      (cond
  21.          (   (< r 2)
  22.              (mapcar 'list l)
  23.          )
  24.          (   l
  25.              (append
  26.                  (mapcar '(lambda ( x ) (cons (car l) x)) (nCr (cdr l) (1- r)))
  27.                  (nCr (cdr l) r)
  28.              )
  29.          )
  30.      )
  31.   )
  32.  
  33.   (defun cxc ( c1 c2 / a1 a2 b1 b2 r i )
  34.     (setq a1 (car c1) b1 (cadr c1))
  35.     (setq a2 (car c2) b2 (cadr c2))
  36.     (setq r (- (* a1 a2) (* b1 b2)))
  37.     (setq i (+ (* a1 b2) (* a2 b1)))
  38.     (list r i)
  39.   )
  40.  
  41.   (setq k 0)
  42.   (setq lst (mapcar '(lambda ( x ) (list (setq k (1+ k)) x)) lst))
  43.   (setq k 0)
  44.   (setq rr (cons (list 1.0 0.0) rr))
  45.   (repeat (length lst)
  46.     (setq r (nCr (mapcar 'car lst) (setq k (1+ k))))
  47.     (setq r (if (= (rem k 2) 1) (mapcar '(lambda ( x ) (* x -1)) (apply 'mapcar (cons '+ (mapcar '(lambda ( x ) (if (= (length x) 1) (mapcar 'float (cadr (assoc (car x) lst))) (progn (setq c (mapcar 'float (cadr (assoc (car x) lst)))) (foreach cc (mapcar '(lambda ( y ) (cadr (assoc y lst))) (cdr x)) (setq c (cxc c cc))) c))) r)))) (apply 'mapcar (cons '+ (mapcar '(lambda ( x ) (if (= (length x) 1) (cadr (assoc (car x) lst)) (progn (setq c (cadr (assoc (car x) lst))) (foreach cc (mapcar '(lambda ( y ) (cadr (assoc y lst))) (cdr x)) (setq c (cxc c cc))) c))) r)))))
  48.     (setq rr (cons r rr))
  49.   )
  50.   (reverse rr)
  51. )
  52.  
Marko Ribar, d.i.a. (graduated engineer of architecture)

:)

M.R. on Youtube