@helter

Лисп лиспом, но линейную алгебру надо знать!

Десять уравнений ― это маленькие системы. В численных методах счёт уравнений может идти на десятки и сотни тысяч. Но это отдельная наука, я плохо её знаю. Тем не менее, по Крамеру всё равно не хочется решать ничего выше второго порядка. Нерациональный он.

Вообще, у систем может быть разное поведение. Они могут не иметь решения. Если у них есть решения, то либо ровно одно, либо бесконечно много. Во втором случае все решения можно получить, выразив так называемые главные неизвестные через свободные, и давая свободным произвольные значения. Какое API вы предполагаете для решения этих вопросов, трудно сказать. Мне кажется, имеет смысл написать две функции: приведение матрицы к ступенчатому виду и к «приведённому» ступенчатому виду (не припомню, как он по-русски, а по-английски ― reduced row echelon form). В частности, по приведённому ступенчатому виду расширенной матрицы системы невооружённым глазом видно решение системы. Неприведённый ступенчатый вид матрицы можно использовать и для других целей: вычисление ранга, например, или определителя, если при получении этого вида не используется умножение строк на число.

Начнём программу на лиспе так:
Пакет ― это навроде пространства имён. Когда вы открываете лисп, пакетом по умолчанию является CL-USER. Однако когда вы пишете что-то, вы обычно определяете свои пакеты и работаете в них. В данном случае мы определяем пакет GAUSS, импортируем в него пакет CL (это никнейм для CL-USER) и экспортируем несколько символов.

Вообще, чтобы получить доступ к символу mysym из пакета mypackage, можно использовать префиксную запись mypackage::mysym. Пакет может экспортировать некоторые символы. Такие символы называются внешними, и в префиксной записи можно ставить одно двоеточие: mypackage:my-external-sym. Внешние символы можно рассматривать как API пакета. Когда вы видите в коде префикс с двойным двоеточием, это всегда подозрительно: почему пользуются неэспортированным символом? Не деталь ли это реализации, которая может измениться в следующем релизе?

Если один пакет импортирует другой, то все внешние символы импортируемого пакета можно использовать без префиксов. Это простая, но удобная система.

Мы импортируем пакет CL, чтобы можно было без префиксов использовать обычные лисповые функции и т. п. Экспортируем мы имена функций, представляющие ценность для внешнего мира. Это row-echelon и reduced-row-echelon, вычисляющие соответственно ступенчатую и приведённую ступенчатую формы, а также некоторые утилиты: число строк и столбцов матрицы и элементарные операции со столбцами.

Потом мы пишем форму in-package, которая делает пакет GAUSS текущим.

Обычно все формы defpackage собирают в отдельном файле, по традиции называемом package.lisp, а файлы с кодом начинают сразу строчкой in-package. Однако если вас интересует организация библиотек в виде файлов, лучше в отдельной теме обсудить. Пока и так сойдёт.

Теперь переходим к коду.

В лиспе есть многомерные массивы, поэтому, не мудрствуя, будем работать с двумерными массивами.

Размеры многомерного массива можно достать функцией array-dimension. Раз мы работаем с матрицами, удобно завести функции row-dimension и column-dimension, возвращающие число строк и столбцов. Функции маленькие, поэтому заинлайним их (предварительная оптимизация, так сказать).

Теперь ― метод Гаусса.

Мы имеем право применять три вида элементарных преобразований строк: перестановка двух строк, умножение строки на число, и прибавление к строке другой строки, умноженной на число. Нетрудно понять, что соответствующие операции для систем уравнений не влияют на множество решений (если только во втором случае не умножать на 0).

Напишем функции, которые выполняют эти преобразования.

Поскольку создавать на каждый чих новую матрицу накладно, мы будем преобразовывать матрицы in place, то есть изменять их сами. Чтобы не забывать, что функции изменяют свои аргументы (что для функций, вообще говоря, нетипично), мы в докстринг включим слово destructively.

При этом функциям неплохо и возвращать что-нибудь. Договоримся, что они будут возвращать саму матрицу.
Реализация очевидна: счётчик цикла проходит от 0 до числа числа столбцов минус 1 (нумерация массивов в лиспе всегда с нуля) и выполняются соответствующие присваивания.

dotimes ― один из немудрёный операторов цикла в лиспе, который запускает счётчик всегда с нуля и перебирает указанное число значений счётчика. То есть dotimes i n эквивалентно for (i = 0; i < n; ++i). В лиспе есть обычай отовсюду возвращать значения, это очень удобно. По умолчанию dotimes возвращает неинтересное значение nil, но факультативный третий аргумент в шапке является выражением, которое нужно вернуть. Мы указываем вернуть a, то есть саму матрицу.

Вообще, функция возвращает значение последней вычисленной формы. В наших функциях по одной форме ― dotimes; их значение и будет возвращено.

Доступ к элементам массива осуществляется функцией aref: aref a i j эквивалентно a[i, j].

Когда мы меняем местами два элемента, мы можем использовать временную переменную для хранения одной из них. Вместо этого мы воспользуемся оператором параллельного присваивания psetf. Обычный оператор присваивания ― setf.

incf эквивалентно +=. К сожалению, аналога *= нет, поэтому в multiply-row честно пишем выражение полностью.

Каков алгоритм приведения матрицы к ступенчатому виду?

Берём первый нулевой столбец. Если в нём одни нули, переходим к следующему столбцу. Если в нём есть ненулевой элемент, то его можно считать первым нулевым по номеру, чего всегда можно добиться перестановкой строк. Тогда мы зануляем все элементы столбца № 0, стоящие ниже ненулевого, путём прибавления к последующим строкам строки 0, умноженной на подходящие числа. Занулив столбец, двигаемся дальше, зануляя следующие столбцы.

Я не в ударе, почитайте любой учебник линейной алгебры.

Напишем вспомогательные функции.

Сначала ― занулить все элементы под a[i,j]. Здесь нужно делать цикл по номерам строк от i+1 до числа строк. К сожалению, dotimes умеет только с нуля. Поэтому мы воспользуемся мудрёным макросом ― loop с ключевыми словами. Он представляет большие возможности ― например, больше, чем сишный for. На самом деле это специальный язык для циклов ― язык внутри языка, который надо учить отдельно. Он не без недостатков и не всем нравится, но часто удобен. Здесь мы используем его наподобие паскалевского for; однако, below означает, что счётчик принимает только значения, строго меньшие верхнего предела.
Клауза finally выполняется в конце итерации; в нашем случае она выполняет форму (return a), возвращающую из цикла значение a. Так как цикл ― единственная форма в функции, его значение становится значением функции.

Заодно напишем зануление элементов над a[i,j], пригодится для обратного хода метода Гаусса. Его можно было бы сделать с помощью dotimes, но я просто подправил код предыдущей функции. В for нижний предел 0 можно не писать.
Теперь нам надо уметь искать ненулевые элементы в столбце j, начиная со строки i. Если все нулевые, вернём ложь ― nil. Значение nil удобно возвращать в случае «неудачи» ― конечно, если оно не может появиться в «удачном» случае.
unless ― одна из условных клауз цикла loop. В данном случае мы перебираем номера строки, если находим ненулевой (zerop проверяет, нулевой ли у неё аргумент), возвращаем номер строки из цикла; если же мы не выйдем из цикла досрочно, то нулевых элементов не было, и возвращаем nil. Вообще, насчёт loop почитайте «Практический Common Lisp» или ещё что-то; поначалу можете использовать более ортодоксальные средства.

Сестра-близнец пригодится для обратного хода. (Вообще, мне начинает не нравиться дублирование кода.)
Теперь можно написать приведение к ступенчатому виду. Опять не на лиспе, а на «лупе». А что делать, задача простая, продвинутые средства для неё не нужны, достаточно массивов и циклов. Мы не будем использовать умножение строк на числа, чтобы функцию можно было использовать для вычисления определителей.
Реализовано примерно то, что я пытался написать выше. Мы начинаем со строки и столбца 0. Если в текущем столбце нет ненулевых элементов, увеличиваем номер текущего столбца. Если есть ненулевой элемент, но не в текущей строке, меняем строки; убедившись, что в текущих строке и столбце не 0, зануляем столбец под этим элементом и увеличиваем на 1 текущую строку и столбец.

Мы используем следующие возможности цикла loop. Во-первых, with ... = вводит локальные переменные для цикла. Во-вторых, мы совмещаем одновременно делаем цикл типа while, и цикл по номерам строк с ненулевыми элементами. Дело в том, что клауза for может иметь вид типа for i from 1 to 10, а может иметь вид типа for i = выражение, в последнем случае на каждой итерации переменной i присваивается значение выражения. Тело цикла у нас разбито в два do. Сначала ― что делать, если нашли-таки строку с ненулевым элементом (я воспользовался лисповым оператором while, эквивалентным if с одной веткой, а можно было воспользоваться одноимённой клаузой loop); а во втором do ― увеличение текущего столбца, которое нужно делать при любом раскладе. Да, incf по умолчанию увеличивает на 1.

Следующая функция выполняет обратный ход метода Гаусса в предположении, что аргумент уже в ступенчатом виде. Алгоритм простой: в каждой строке ищется ненулевой элемент и зануляются все элементы над ним. Здесь мы совмещаем два счётчика и цикл while (если в строке уже нет ненулевых элементов, нужно заканчивать).
Наконец, следующая функция возвращает приведённую ступенчатую форму. Она, в частности, и решает уравнения.
Пример 567 из Проскурякова:

Это значит, исходная система эквивалентна следующей:

Нетрудно понять, какое решение.

Неопределённая система, № 690:

Значит, система эквивалентна такой:

Кто немного знает линейную алгебру, понимает, что второе и четвёртое неизвестные ― свободные, они могут принимать любые значения, и общее решение имеет ви

Из приведённой формы всё видно невооружённым глазом.

У Проскурякова в ответе в качестве главных взяты x1 и x2. С его ответом можно сравнить, если эти столбцы записать последними в матрице системы.

Несовместная система:

Исходная система эквивалентна такой:

несовместность которой очевидна.

Учите алгебру!

Вернуться к обсуждению:
Метод Крамера или метод Гаусса. Реализация Lisp

2