[eugen81]

Доброго времени суток!
Уже которую неделю ломаю голову, а не могу решить на лиспе казалось бы простую задачу.
Собственно задача:
Имеется матрица X, над которой надо произвести разложение Холецкого и получить матрицу L.
Программа в нынешнем виде может считать каждую строку матрицы L и надо бы их объединить в одно целое, но проблема состоит в следующем.
если нумеровать строки матрицы с нуля, то:
сначала вычисляем строку L(0)
затем L(1), причем для ее расчета используется L(0)
далее L(2) и для ее расчета используются уже обе строки L(0,1)
и так далее в том же духе.

Функции вычисления отдельной строки выглядят следующим образом:

Код:

[Выделить все]

(defun mrowL (j X L)
  (list (mrowL_ 0 j X L))
  );_defun

(defun mrowL_ (i j X L)
  (cond
   ((= i (length X)) nil)
   (t
    (cons
     (evLij i j X L)
     (mrowL_ (1+ i) j X L)
     );_cons
    );t
   );_cond
  );_defun

Здесь:
X - упомянутая выше исходная матрица Х,
L - часть матрицы L, используемая для вычисления элемента, а значит и строки.
Функция evLij вычисляет каждый элемент матрицы L,
i, j - позиция элемента в матрице (номер столбца и строки соответственно).

Итог:
1 строка: (setq r0 (mrowL 0 Х nil))
2 строка: (setq r1 (mrowL 1 X (append nil r0)))
3 строка: (setq r2 (mrowL 2 X (append r0 r1)))
4 строка: (setq r3 (mrowL 3 X (append r0 r1 r2)))
и так далее

Для конкретного числа строк программа составляется очень легко, но хотелось бы иметь универсальную программу

P.S. забыл оговориться, что матрицы имеют вид (list (list x11 x21 x31...) (list x12 x22 x32...)...)

[kp+]

Рекурсия нужна. Поищи по этому слову.
И на caduser.ru была тема "Уроки создания рекурсивных функций"

[eugen81]

Вот я и пытаюсь сделать рекурсию, у меня почти все функции обработки матриц сделаны через рекурсию. но именно эта часть задачи не получается.
Кстати, если интересно - выкладываю лисп, который в готовом виде должен находить функцию из ряда координат точек, какие довольно часто встречаются в справочниках по гидравлике и всяких СНиПах.
Жду критики кода и предложений по его улучшению

[Елпанов Евгений]

Сообщение от eugen81 Вот я и пытаюсь сделать рекурсию, у меня почти все функции обработки матриц сделаны через рекурсию. но именно эта часть задачи не получается.

Напиши подробно, что и как должна делать программа с примерами матрицы и пошаговых вычислений. Я помогу тебе собрать универсальную программу.
Извини, но самостоятельно разбираться с необходимым тебе алгоритмом, мне некогда...

[eugen81]

Евгений, спасибо за отзыв!
Постараюсь описать, что требуется от программы.
Наверное стоит это делать с момента, когда исходная квадратная матрица А уже готова и имеет вид списка

Код:

[Выделить все]

((а11 а21 а31 a41)
 (а12 а22 а32 a42)
 (а13 а23 а33 a43)
 (a14 a24 a34 a44))

,
причем размерность ее может любой.

Чтобы получить матрицу L используются формулы во вложении.
Первая формула определяет элементы главной диагонали L, т.е. L11, L22, L33 и т.д.
Вторая формула определяет все остальные элементы, причем элементы матрицы, находящиеся ниже главной диагонали равны нулю.
После расчета матрица должна принять вид

Код:

[Выделить все]

((L11 L21 L31 L41)
 (0   L22 L32 L42)
 (0   0   L33 L43)
 (0   0   0   L44)

Исходя из 1й формулы элемент L11 равен квадратному корню а11.
По второй формуле находим элементы первой строки:
L21 = 1/L11 * a21
L31 = 1/L11 * a31
L41 = 1/L11 * a41

Далее 2я строка:
L22 = sqrt(a22 - L21^2)
L32 = 1/L22 * (a32 - L31*L21)
L42 = 1/L22 * (a42 - L41*L21)

третья строка:
L33 = sqrt(a33 - L31^2 - L32^2)
L43 = 1/L43 * (a43 - L41*L31 - L42*L32)

и четвертая строка:
L44 = sqrt(a44 - L41^2 - L42^2 - L43^2)

Теперь что сделал я)

Элемент Lii на главной диагонали находится как квадратный корень из aii минус сумму квадратов элементов, расположенных над Lii.
Эту сумму квадратов вычисляет функция Li:

Код:

[Выделить все]

(defun Li (i mtrl)
"Возвращает сумму квадратов элементов Lij над элементом Lii (i - позиция колонки в матрице начиная с нуля)
 Функция используется для вычичления очередного элемента Lii"
  (Li_ i i mtrl)
  );_defun

  (defun Li_ (i j mtrl)
    (cond
     ((= j 0) 0)
     (t
      (+
       (expt (nth i (nth (1- j) mtrl)) 2);_expt Элемент над Lii возводим в квадрат
       (Li_ i (1- j) mtrl)
       );+
      );t
     );_cond
    );_defun

Вычисление элементов Lii выполняет одноименная функция из моего лиспа (только нумерация строк и столбцов здесь начинается с нуля):

Код:

[Выделить все]

(defun Lii (i mtrx mtrl)
"Вычисляет элемент Lii (i - позиция в матрице начиная с нуля)"
  (cond
   ((= i 0) (expt (nth i (nth i mtrx)) 0.5)) ; L00 равен квадратному корню из а00
   (t (expt (-
	     (nth i (nth i mtrx))
	     (Li i mtrl)
	     );-
	    0.5);_expt
      );t
   );_cond
  );_defun

Вычисление элементов Lij:

Код:

[Выделить все]

(defun Lij (i j mtrx mtrl)
"Вычисляет элемент Lij (i<>j - позиция в матрице начиная с нуля)"
  (cond
   ((= j 0) (/                        ; для 0й строки Li0 = ai0/L00
	     (nth i (nth j mtrx))
	     (Lii j mtrx nil)
	     );/
    );=
   (t (/
       (-
	(nth i (nth j mtrx))
	(sumLik i j j mtrl)
	);-
       (Lii j mtrx mtrl)
       );/
      );_t
   );_cond
  );_defun

Здесь sumLik - функция вычисляющая сумму произведений элементов Lik*Ljk во 2й формуле:

Код:

[Выделить все]

  (defun sumLik (i j k mtrl)
    (cond
     ((= k 0) 0)
     (t (+
	 (*
	  (nth i (nth (1- k) mtrl))
	  (nth j (nth (1- k) mtrl))
	  );*
	 (sumLik i j (1- k) mtrl)
	 );+
	);t
     );_cond
    );_defun

Затем все это я объединил в одну функцию evLij, которая определяет позицию искомого элемента и выполняет соответствующую функцию:

Код:

[Выделить все]

(defun evLij (i j mtrx mtrl)
"Вычисляет элемент матрицы Lij методом Холецкого (i,j - позиция в матрице начиная с нуля)"
  (cond
   ((= i j)
    (Lii i mtrx mtrl)
    );1
   ((> i j)
    (Lij i j mtrx mtrl)
    );2
   (t 0)
   );_cond
  );_defun

И следующая функция собирает элементы в одну строку:

Код:

[Выделить все]

(defun mrowL (j mtrx mtrl)
"Создает строку матрицы L"
  (list (mrowL_ 0 j mtrx mtrl))
  );_defun

  (defun mrowL_ (i j mtrx mtrl)
    (cond
     ((= i (length mtrx)) nil)
     (t
      (cons
       (evLij i j mtrx mtrl)
       (mrowL_ (1+ i) j mtrx mtrl)
       );_cons
      );t
     );_cond
    );_defun

З.Ы. Наконец-то дошло как применять лямбда-функцию, скоро перепишу код, чтобы убрать вспомогательные функции типа mrowL_

[Елпанов Евгений]

Жаль, ты не дал пример исходной матрицы и необходимого результата.
результат не проверял, вычисления выполнял на своей тестовой матрице...

Код:

[Выделить все]

(defun ai (a b) (sqrt (apply (function -) (cons a (mapcar (function *) b b)))))
(defun aii (a b c)
 (cons
  a
  (mapcar
   (function (lambda (b c d) (/(apply (function -) (cons b (mapcar (function *) c d)))a)))
   b
   (cdr c)
   (cddr c)
  )
 )
)
(defun cholesky (a l / ll)
 (cond
  ((not l)
   (cholesky (mapcar (function (lambda (a) (cons 0. a))) a)
        (list (cons 0. (mapcar (function (lambda (a) 0.)) (car a))))
   )
  )
  (a
   (cons (setq ll (aii (ai (cadar a) (mapcar (function cadr) l))
                       (cddar a)
                       (apply (function mapcar) (cons 'list l))
                  )
         )
         (mapcar (function (lambda (a) (cons 0. a)))
                 (cholesky (mapcar (function cdr) (cdr a)) (cons ll (mapcar (function cdr) l)))
         )
   )
  )
  (nil)
 )
)

Пример вызова:

Код:

[Выделить все]

(cholesky '((81. -45. 45.) (-45. 50. -15.) (45. -15. 38.)) nil)

возвращает:

Код:

[Выделить все]

'((9. -5. 5.) (0. 5. 2.) (0. 0. 3.))

проверка:
Правая часть матрицы:
b = '(531. -460. 193.)
Ly = b
y = '(59. -33. -12.)
Lty = x
Lt = (apply
(function mapcar)
(cons (function list) (cholesky '((81. -45. 45.) (-45. 50. -15.) (45. -15. 38.)) nil))
)

x = '(6. -5. -4.)

[eugen81]

Спасибо за помощь, все получилось! теперь осталось разобраться)

[CB]

Женя!

Жаль, ты не дал пример исходной матрицы и необходимого результата.

Исходную матрицу можно легко составить и самому:
1. Результирующая матрица

Код:

[Выделить все]

 
(setq lst
       '((2.0 1.0 3.0 4.0 1.0)
         (0.0 3.0 2.0 3.0 2.0)
         (0.0 0.0 1.0 2.0 5.0)
         (0.0 0.0 0.0 5.0 3.0)
         (0.0 0.0 0.0 0.0 2.0)
        )
) ;_ end of setq

2. Транспонированная матрица

Код:

[Выделить все]

(defun transpon-matrix (matrix)
  (apply 'mapcar (cons 'list matrix))
)
(setq lst_t (transpon-matrix lst))
=>
((2.0 0.0 0.0 0.0 0.0)
   (1.0 3.0 0.0 0.0 0.0)
   (3.0 2.0 1.0 0.0 0.0)
   (4.0 3.0 2.0 5.0 0.0)
   (1.0 2.0 5.0 3.0 2.0)
 )

3. Перемножаем lst_t на lst

Код:

[Выделить все]

(defun multiplication-matrix (a b)
  (mapcar '(lambda (x)
             (mapcar '(lambda (y)
                        (apply '+ (mapcar '* x y))
                      ) ;_ end of lambda
                     (apply 'mapcar (cons 'list b))
             ) ;_ end of mapcar
           ) ;_ end of lambda
          a
  ) ;_ end of mapcar
) ;_ end of defun
(setq matr (multiplication-matrix lst_t lst))
=>
((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)
)

Это и есть исходная матрица...
Если подставить эту матрицу в твой код, то получается

Код:

[Выделить все]

(cholesky matr nil)
=>
$ 
; error: function undefined for argument: -15.8889
_$

Накидал свой вариант:

Код:

[Выделить все]

(defun cholesky (lst / str temp f test)
  (defun test (lst)
    (if lst
      (cons (f (car lst))
            (mapcar '(lambda (x) (cons 0. x))
                    (test (mapcar 'cdr (cdr lst)))
            ) ;_ end of mapcar
      ) ;_ end of cons
      lst
    ) ;_ end of if
  ) ;_ end of defun
  (defun f (c / b)
    (if str
      (setq temp (mapcar 'cdr (cons str temp)))
    ) ;_  if
    (setq
      b (sqrt
          (- (car c)
             (apply '+
                    (mapcar '(lambda (x) (* x x))
                            (mapcar 'car temp)
                    ) ;_ end of mapcar
             ) ;_ end of apply
          ) ;_ end of -
        ) ;_ end of sqrt
    ) ;_ end of setq
    (setq str
           (cons
             b
             (if str
               (mapcar
                 '(lambda (x y) (/ (- x y) b))
                 (cdr c)
                 (apply
                   'mapcar
                   (cons
                     '+
                     (mapcar
                       '(lambda (x)
                          (mapcar '(lambda (y) (* (car x) y))
                                  (cdr x)
                          ) ;_  mapcar
                        ) ;_  lambda
                       temp
                     ) ;_  mapcar
                   ) ;_  cons
                 ) ;_  apply
               ) ;_  mapcar
               (mapcar '(lambda (x) (/ x b)) (cdr c))
             ) ;_  if
           ) ;_ end of cons
    ) ;_  setq
  ) ;_  defun
  (test lst)
) ;_ end of defun

Код:

[Выделить все]

(cholesky matr)
=>
((2.0 1.0 3.0 4.0 1.0) (0.0 3.0 2.0 3.0 2.0) (0.0 0.0 1.0 2.0 5.0) (0.0 0.0 0.0 5.0 3.0) (0.0 0.0 0.0 0.0 2.0))

Естественно в этот код нужно добавить как минимум две проверки - на симметричность исходной матрицы и на положительность элементов главной диагонали разложенной матрицы...
Если первая проверка решается легко, со второй лично у меня пока заминка, т.к. такая проверка производится непосредственно в процессе вычислений и не получается выйти из рекурсии...

[Елпанов Евгений]

CB, спасибо, поправил код...

Код:

[Выделить все]

(defun ai (a b) (sqrt (apply (function -) (cons a (mapcar (function *) b b)))))
(defun aii (a b c)
 (cons a
       (mapcar
        (function
         (lambda (b d) (/ (apply (function -) (cons b (mapcar (function *) (car c) d))) a))
        )
        b
        (cdr c)
       )
 )
)
(defun cholesky (a l / ll)
 (cond
  ((not l) (cholesky a (list (cons 0. (mapcar (function (lambda (a) 0.)) (car a))))))
  (a
   (cons
    (setq ll (aii (ai (caar a) (mapcar (function car) l))
                  (cdar a)
                  (apply (function mapcar) (cons 'list l))
             )
    )
    (mapcar (function (lambda (a) (cons 0. a)))
            (cholesky (mapcar (function cdr) (cdr a)) (mapcar (function cdr) (cons ll l)))
    )
   )
  )
  (nil)
 )
)

[eugen81]

Ух-ты) не ожидал такой помощи! Всем большое спасибо!

[ETCartman]

FEM - Лисп, лицензия BSD
http://savannah.nongnu.org/projects/femlisp

[eugen81]

Спасибо, посмотрим
Но все же и самому неплохо уметь такое делать

[ETCartman]

BSD это самая вольная лицензия. "Бери, используй, перерабатывай, продавай (можешь стать миллионером )