Как решить слау в питоне

Как решить систему уравнений в Python (3 примера)

Чтобы решить систему уравнений в Python, мы можем использовать функции из библиотеки NumPy .

В следующих примерах показано, как использовать NumPy для решения нескольких различных систем уравнений в Python.

Пример 1. Решение системы уравнений с двумя переменными

Предположим, у нас есть следующая система уравнений, и мы хотели бы найти значения x и y:

В следующем коде показано, как использовать NumPy для поиска значений x и y:

import numpy as np #define left-hand side of equation left_side = np.array([[5, 4], [2, 6]]) #define right-hand side of equation right_side = np.array([35, 36]) #solve for x and y np.linalg.inv (left_side). dot (right_side) array([3., 5.])

Это говорит нам о том, что значение x равно 3 , а значение y равно 5 .

Пример 2. Решение системы уравнений с тремя переменными

Предположим, у нас есть следующая система уравнений, и мы хотели бы найти значения x, y и z:

В следующем коде показано, как использовать NumPy для определения значений x, y и z:

import numpy as np #define left-hand side of equation left_side = np.array([[4, 2, 1], [3, 5, -2], [2, 2, 4]]) #define right-hand side of equation right_side = np.array([34, 41, 30]) #solve for x, y, and z np.linalg.inv (left_side). dot (right_side) array([5., 6., 2.])

Это говорит нам о том, что значение x равно 5 , значение y равно 6 , а значение z равно 2 .

Пример 3. Решение системы уравнений с четырьмя переменными

Предположим, у нас есть следующая система уравнений, и мы хотели бы найти значения w, x, y и z:

6ш + 2х + 2у + 1з = 37

2ш + 1х + 1у + 0з = 14

3ш + 2х + 2у + 4з = 28

2ш + 0х + 5у + 5з = 28

В следующем коде показано, как использовать NumPy для определения значений w, x, y и z:

import numpy as np #define left-hand side of equation left_side = np.array([[6, 2, 2, 1], [2, 1, 1, 0], [3, 2, 2, 4], [2, 0, 5, 5]]) #define right-hand side of equation right_side = np.array([37, 14, 28, 28]) #solve for w, x, y, and z np.linalg.inv (left_side). dot (right_side) array([4., 3., 3., 1.])

Это говорит нам о том, что значение w равно 4 , x равно 3 , y равно 3 и z равно 1 .

Дополнительные ресурсы

В следующих руководствах объясняется, как решить систему уравнений с помощью другого статистического программного обеспечения:

Numpy: решение системы линейных уравнений методом Гаусса

Пытаюсь в Python создать алгоритм расчета СЛАУ методом Гаусса. Метод заключается в следующем. Составляем матрицу коэффициентов, включая свободный член. Затем приводим матрицу к треугольному виду. Для этого сначала по первому столбцу (с индексом 0) каждый элемент делим на диагональный (a0,0) (в примере — 3,8), вычисляя индекс m, а после пересчитываем строку 2 по формуле: каждый ее элемент (без элемента свободного члена из последнего столбца) минус произведение элемента над ним (из нулевой строки) и индекса m для второй строки. Отдельно отработаем со столбцом свободного члена (здесь алгоритм неважен). Следом аналогичные действия надо проделать для третьей строки элементов (но учитывая, что на первой итерации элементы второй строки преобразованы вышеописанным алгоритмом, а коэффициент m будет считаться по второму столбцу: соответственно делим все его элементы на диагональный элемент из 2-й строки a1,1) (в примере 1,3). Вопрос: я рассчитал вектор-столбец m: m = ([1,000, 1,684, 0,632]) И теперь надо отработать со второй строкой матрицы. И вот здесь сложность с индексацией. Во-первых, не могу перебрать значения m, тип которых float. Во-вторых, неверно задаю индексацию элементов второй строки (по сути — после нулевой это первая строка)

import numpy as np matrix = np.array([[3.8, 6.7, -1.2, 5.2], [6.4, 1.3, -2.7, 3.8], [2.4, -4.5, 3.5, -0.6]]) def gaussFunc(matrix): # расчет len1 (3) и len2 (4) - здесь не приводится # код расчета m по нулевому столбцу: for row in range(len1): for column in range(len2-3): m = matrix[row][column] / matrix[column][column] elem = row-1 # значения столбцов по нулевой строке кладем в # переменную elem for i in range(len(m)-1): # идем циклом по диапазону трех значений m минус #последнее третье: ошибка по len для float while row < (len1-1): # пока строка первая или вторая (в len2 их всего # 3): while column < (len2-1): # пока колонка первая, вторая или третья # (минус столбец свободного # члена): # пересчитанные коэффициенты второй (первой в numpy) строки: # текущий элемент - m по данной строке*верхний элемент в данном # столбце (из строки 0): a = matrix[row][column] - m[i]*matrix[elem][column]

Отслеживать
15.4k 1 1 золотой знак 8 8 серебряных знаков 32 32 бронзовых знака
задан 17 апр 2021 в 11:08
Alex_Kazantsev Alex_Kazantsev
629 1 1 золотой знак 8 8 серебряных знаков 19 19 бронзовых знаков

2 ответа 2

Сортировка: Сброс на вариант по умолчанию

В конце приведена ссылка на jupyter notebook с более-менее полным решателем СЛАУ. Плюс трюк, как считать на pyton почти так же быстро, как на Фортране 🙂

Если не обращать внимание на возможное деление на ноль, то привидение к треугольному виду можно записать очень просто:

def gaussFunc(matrix): # функция меняет матрицу через побочные эффекты # если вам нужно сохранить прежнюю матрицу, скопируйте её np.copy for nrow, row in enumerate(matrix): # nrow равен номеру строки # row содержит саму строку матрицы divider = row[nrow] # диагональный элемент # делим на диагональный элемент. row /= divider # теперь надо вычесть приведённую строку из всех нижележащих строчек for lower_row in matrix[nrow+1:]: factor = lower_row[nrow] # элемент строки в колонке nrow lower_row -= factor*row # вычитаем, чтобы получить ноль в колонке nrow # все строки матрицы изменились, в принципе, можно и не возвращать return matrix

Результат для вашего примера:

array([[ 1. , 1.76315789, -0.31578947, 1.36842105], [-0. , 1. , 0.06800211, 0.49657354], [ 0. , 0. , 1. , 0.09309401]])

В чём засада. В ходе вычислений может оказаться ноль на диагонали.

matrix = np.array([[1, 0, 0, 1], [0, 0, 1, 2], [0, 1, 0, 3]])

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

Функция make_identity берёт матрицу, полученную методом Гаусса, и доводит её до единичной. Для этого строки перебираются снизу вверх и вычитаются из вышележащих строк, чтобы обнулить соответствующие колонки.

def make_identity(matrix): # перебор строк в обратном порядке for nrow in range(len(matrix)-1,0,-1): row = matrix[nrow] for upper_row in matrix[:nrow]: factor = upper_row[nrow] upper_row -= factor*row return matrix

Результат для вашего примера: make_identity(gaussFunc(np.copy(matrix)))

array([[ 1. , 0. , 0. , 0.53344344], [-0. , 1. , 0. , 0.49024295], [ 0. , 0. , 1. , 0.09309401]])

Вырезаем последний столбец, получим строку корней: roots = make_identity(gaussFunc(np.copy(matrix)))[:,3]

array([0.53344344, 0.49024295, 0.09309401])

Умножаем найденные корни на исходную матрицу и сравниваем с последним столбцом исходной задачи: np.matmul(matrix[. 3], roots.T) - matrix[:,3]

Результат вычисления array([ 0.00000000e+00, -4.44089210e-16, -2.22044605e-16])

Следовательно, корни найдены правильно. С чем и поздравляю.

UPDATE

Метод Гаусса с выбором главного элемента. Плюс добавлена обработка нуля на диагонали.

def gaussPivotFunc(matrix): for nrow in range(len(matrix)): # nrow равен номеру строки # np.argmax возвращает номер строки с максимальным элементом в уменьшенной матрице # которая начинается со строки nrow. Поэтому нужно прибавить nrow к результату pivot = nrow + np.argmax(abs(matrix[nrow:, nrow])) if pivot != nrow: # swap # matrix[nrow], matrix[pivot] = matrix[pivot], matrix[nrow] - не работает. # нужно переставлять строки именно так, как написано ниже matrix[[nrow, pivot]] = matrix[[pivot, nrow]] row = matrix[nrow] divider = row[nrow] # диагональный элемент if abs(divider) < 1e-10: # почти нуль на диагонали. Продолжать не имеет смысла, результат счёта неустойчив raise ValueError(f"Матрица несовместна. Максимальный элемент в столбце : ") # делим на диагональный элемент. row /= divider # теперь надо вычесть приведённую строку из всех нижележащих строчек for lower_row in matrix[nrow+1:]: factor = lower_row[nrow] # элемент строки в колонке nrow lower_row -= factor*row # вычитаем, чтобы получить ноль в колонке nrow # приводим к диагональному виду make_identity(matrix) return matrix

В этой функции главный фокус в том, как переставлять строки. Простая формула matrix[nrow], matrix[pivot] = matrix[pivot], matrix[nrow] не работает. При таком присваивании справа стоят указатели на строку, а слева - адреса, куда нужно скопировать значения. То есть при первом присваивании в строку nrow будет скопирована строка pivot , а в строку pivot - содержимое строки nrow -- но оно уже переписано из строки pivot ! То есть фактически строка pivot скопируется в саму себя. И вместо перестановки двух строк будет копия одной строки.

matrix[[nrow, pivot]] = matrix[[pivot, nrow]] - работает. И c явным копированием тоже работает: matrix[nrow], matrix[pivot] = matrix[pivot], np.copy(matrix[nrow])

UPDATE 2

Помимо собственно решателя дано сравнение с Си-шным решателем numpy.linalg.solve и трюк как ускорить скорость счёта решателя на пайтоне в 20 раз, что почти так же быстро, как чисто Си-шный решатель.

Строго говоря, решатель в numpy даже не Си-шный, а фортрановский LAPACK gesv

Как я могу решить уравнения в Python?

Скажем, у меня есть уравнение: 2x + 6 = 12 С алгеброй видно, что x = 3 . Как я могу сделать программу на Python, которая может решить для x ? Я новичок в программировании, и я смотрел на eval() и exec() , но я не могу понять, как заставить их делать то, что я хочу. Я не хочу использовать внешние библиотеки (например, SAGE), я хочу сделать это в простом Python.

user1221937 08 май 2012, в 16:09
Поделиться

Я хочу автомобиль, который делает 0 до 60 за 4,5 секунды и получает 45 миль на галлон. Может быть, вы могли бы рассмотреть возможность удаления требования для простого Python и использовать внешние библиотеки

Mike Pennington 08 май 2012, в 13:51

Если вы хотите решить какое-либо уравнение, вам придется создать свою собственную библиотеку. Также 4,5 секунды недостаточно для этого примера: D

jamylak 08 май 2012, в 13:52
Будут ли проблемы всегда выглядеть так: solve y = mx + c for x ?
Li-aung Yip 08 май 2012, в 13:53

@MikePennington: если он хочет разрешить только очень ограниченный набор уравнений, это вполне выполнимо. В абстрактном случае, конечно, вам придется создать свой собственный символический математический движок в духе maxima или Mathematica , но я не думаю, что это цель.

Li-aung Yip 08 май 2012, в 13:54

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

Lostsoul 08 май 2012, в 13:55

@Li-aungYip Li-aungYip, если он хочет решить только ограниченный набор уравнений, он должен указать это в вопросе

Mike Pennington 08 май 2012, в 13:56
это домашнее задание? или почему вы не хотите использовать внешние библиотеки?
Niek de Klein 08 май 2012, в 14:02
Показать ещё 5 комментариев
Поделиться:
calculator

5 ответов

Как насчет SymPy? Их solver выглядит так, как вам нужно. Посмотрите их исходный код, если вы хотите сами построить библиотеку.

eumiro 08 май 2012, в 15:16
Поделиться

Забавно видеть ответ, подобный этому, через несколько минут после всех глупых комментариев по этому вопросу: D

naught101 30 янв. 2015, в 02:30

Существует два способа подхода к этой проблеме: численно и символически.

Чтобы решить его численно, вы должны сначала закодировать его как функцию "runnable" - вставьте значение в значение, получите значение. Например,

def my_function(x): return 2*x + 6

Весьма возможно проанализировать строку для автоматического создания такой функции; скажем, разберите 2x + 6 в список, [6, 2] (где индекс списка соответствует мощности x - так 6 * x ^ 0 + 2 * x ^ 1). Тогда:

def makePoly(arr): def fn(x): return sum(c*x**p for p,c in enumerate(arr)) return fn my_func = makePoly([6, 2]) my_func(3) # returns 12

Затем вам понадобится другая функция, которая повторно подключает значение x в вашу функцию, анализирует разницу между результатом и тем, что она хочет найти, и изменяет ее значение x (надеюсь) минимизирует разницу.

def dx(fn, x, delta=0.001): return (fn(x+delta) - fn(x))/delta def solve(fn, value, x=0.5, maxtries=1000, maxerr=0.00001): for tries in xrange(maxtries): err = fn(x) - value if abs(err) < maxerr: return x slope = dx(fn, x) x -= err/slope raise ValueError('no solution found')

Здесь есть много потенциальных проблем - найти хорошее начальное значение x, считая, что функция фактически имеет решение (т.е. нет никаких однозначных ответов на x ^ 2 + 2 = 0), попав в пределы точность вычислений и т.д. Но в этом случае функция минимизации ошибок подходит и получается хороший результат:

solve(my_func, 16) # returns (x =) 5.000000000000496

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

Символьный решатель, такой как Mathematica или Maple, представляет собой экспертную систему с множеством встроенных правил ( "знаний" ) об алгебре, исчислении и т.д.; он "знает", что производная от sin is cos, что производная от kx ^ p равна kpx ^ (p-1) и т.д. Когда вы даете ему уравнение, он пытается найти путь, набор правил-приложений, откуда он (уравнение), туда, где вы хотите быть (простейшая форма уравнения, которое, мы надеемся, является решением).

Ваше примерное уравнение довольно простое; символическое решение может выглядеть так:

=> LHS([6, 2]) RHS([16]) # rule: pull all coefficients into LHS LHS, RHS = [lh-rh for lh,rh in izip_longest(LHS, RHS, 0)], [0] => LHS([-10,2]) RHS([0]) # rule: solve first-degree poly if RHS==[0] and len(LHS)==2: LHS, RHS = [0,1], [-LHS[0]/LHS[1]] => LHS([0,1]) RHS([5])

и есть ваше решение: x = 5.

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

Решаем систему линейных алгебраических уравнений с Python-пакетом scipy.linalg (не путать с numpy.linalg)

Аппарат линейной алгебры применяют в самых разных областях — в линейном программировании, эконометрике, в естественных науках. Отдельно отмечу, что этот раздел математики востребован в машинном обучении. Если, например, вам нужно поработать с матрицами и векторами, то, вполне возможно, на каком-то шаге вам придётся решать систему линейных алгебраических уравнений (СЛАУ).

СЛАУ — мощный инструмент моделирования процессов. Если некая хорошо изученная модель на основе СЛАУ годится для формализации задачи, то с высокой вероятностью её можно будет решить. А вот насколько быстро — это зависит от вашего инструментария и вычислительных мощностей.

Я расскажу про один из таких инструментов — Python-пакет scipy.linalg из библиотеки SciPy, который позволяет быстро и эффективно решать многие задачи с использованием аппарата линейной алгебры.

В этом туториале вы узнаете:

как установить scipy.linalg и подготовить среду выполнения кода;
как работать с векторами и матрицами с помощью NumPy;
почему scipy.linalg лучше, чем numpy.linalg;
как формализовать задачи с использованием систем линейных алгебраических уравнений;
как решать СЛАУ с помощью scipy.linalg (на реальном примере).

Когда речь идёт о математике, изложение материала должно быть последовательным — таким, чтобы одно следовало из другого. Эта статья не исключение: сначала будет много подготовительной информации и только потом мы перейдём непосредственно к делу.

Если готовы к этому — приглашаю под кат. Хотя, честно говоря, некоторые разделы можно пропускать — например, основы работы с векторами и матрицами в NumPy (если вы хорошо знакомы с ним).

Установка scipy.linalg

SciPy — это библиотека Python с открытым исходным кодом для научных вычислений: решение СЛАУ, оптимизация, интеграция, интерполяция и обработка сигналов. Помимо linalg, она содержит несколько других пакетов — например, NumPy, Matplotlib, SymPy, IPython и pandas.

Среди прочего, scipy.linalg содержит функции для с работы с матрицами — вычисление определителя, инверсия, вычисление собственных значений и векторов, а также сингулярное разложение.

Чтобы использовать scipy.linalg, вам необходимо установить и настроить библиотеку SciPy. Это можно сделать с помощью дистрибутива Anaconda, а также системы управления пакетами и инсталлятора Conda.

Для начала подготовьте среду выполнения для библиотеки:

$ conda create --name linalg $ conda activate linalg

Установите необходимые пакеты:

(linalg) $ conda install scipy jupyter

Эта команда может работать долго. Не пугайтесь!

Я предлагаю использовать Jupyter Notebook для запуска кода в интерактивной среде. Это не обязательно, но лично мне он облегчает работу.

Перед открытием Jupyter Notebook вам необходимо зарегистрировать экземпляр conda linalg, чтобы использовать его в качестве ядра при создании ноутбука. Для этого выполните следующую команду:

(linalg) $ python -m ipykernel install --user --name linalg

Теперь можно открыть Jupyter Notebook:

$ jupyter notebook

Когда он загрузится в вашем браузере, создайте новый notebook, нажав New → linalg:

Чтобы убедиться, что установка библиотеки SciPy прошла успешно, введите в ноутбуке:

>>> In [1]: import scipy

NumPy для работы с векторами и матрицами (куда же без него)

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

NumPy — это наиболее популярный пакет для работы с матрицами и векторами в Python. Часто его применяют в сочетании с scipy.linalg. Чтобы начать работу с матрицами и векторами, нужно импортировать пакет NumPy:

>>> In [2]: import numpy as np

Для представления матриц и векторов NumPy использует специальный тип, называемый ndarray. Чтобы создать объект ndarray, вы можете использовать array ().

Например, вам нужно создать следующую матрицу:

Создадим матрицу как набор вложенных списков (векторов-строк):

>>> In [3]: A = np.array([[1, 2], [3, 4], [5, 6]]) . A Out[3]: array([[1, 2], [3, 4], [5, 6]])

Заметьте, что приведённый выше вывод (Outp[3]) достаточно наглядно показывает получившуюся матрицу.

И ещё: все элементы матрицы должны и будут иметь один тип. Это можно проверить с помощью dtype.

>>> In [4]: A.dtype Out[4]: dtype('int64')

Здесь элементы являются целыми числами, поэтому их общий тип по умолчанию — int64. Если бы среди них было хотя бы одно число с плавающей точкой, все элементы получили бы тип float64:

>>> In [5]: A = np.array([[1.0, 2], [3, 4], [5, 6]]) . A Out[5]: array([[1., 2.], [3., 4.], [5., 6.]]) In [6]: A.dtype Out[6]: dtype('float64')

Чтобы вывести на экран размерность матрицы, можно использовать метод shape:

>>> In [7]: A.shape Out[7]: (3, 2)

Как и ожидалось, размерность матрицы A 3×2, то есть A имеет три строки и два столбца.

При работе с матрицами часто приходится использовать операцию транспонирования, которая столбцы превращает в строки и наоборот. Чтобы транспонировать вектор или матрицу (представленную объектом типа ndarray), вы можете использовать .transpose () или .T. Например:

>>> In [8]: A.T Out[8]: array([[1., 3., 5.], [2., 4., 6.]])

Чтобы создать вектор, также можно использовать.array(), передав туда список значений в качестве аргумента:

>>> In [9]: v = np.array([1, 2, 3]) . v Out[9]: array([1, 2, 3])

По аналогии с матрицами, используем .shape(), чтобы вывести на экран размерность вектора:

>>> In [10]: v.shape Out[10]: (3,)

Заметьте, что она выглядит как (3,), а не как (3, 1) или (1, 3). Разработчики NumPy решили сделать отображение размерности векторов так же, как в MATLAB.

Чтобы получить на выходе размерность (1, 3), нужно было бы создать вот такой массив:

>>> In [11]: v = np.array([[1, 2, 3]]) . v.shape Out[11]: (1, 3)

Для (3, 1) — вот такой:

>>> In [12]: v = np.array([[1], [2], [3]]) . v.shape Out[12]: (3, 1)

Как видите, они не идентичны.

Часто возникает задача из вектора-строки сделать вектор-столбец. Как вариант, можно сначала создать вектор-строку, а потом использовать .reshape() для его преобразования в столбец:

>>> In [13]: v = np.array([1, 2, 3]).reshape(3, 1) . v.shape Out[13]: (3, 1)

В приведённом выше примере мы использовали .reshape() для получения вектора-столбца с размерностью (3, 1) из вектора с размерностью(3,). Стоит отметить, что .reshape() делает преобразование с учётом того, что количество элементов (100% заполненных мест) в массиве с новой размерностью должно быть равно количеству элементов в исходном массиве.

В практических задачах часто требуется создавать матрицы, полностью состоящие из нулей, единиц или случайных чисел. Для этого NumPy предлагает несколько удобных функций, о которых я расскажу в следующем разделе.

Заполнение массивов данными

NumPy позволяет быстро создавать и заполнять массивы. Например, чтобы создать массив, заполненный нулями, можно использовать .zeros():

>>> In [15]: A = np.zeros((3, 2)) . A Out[15]: array([[0., 0.], [0., 0.], [0., 0.]])

В качестве аргумента .zeros() нужно указать размерность массива, упакованную в кортеж (перечислить значения через запятую и обернуть это в круглые скобки). Элементы созданного массива получат тип float64.

Точно так же, для создания массивов из единиц можно использовать .ones ():

>>> In [16]: A = np.ones((2, 3)) . A Out[16]: array([[1., 1., 1.], [1., 1., 1.]])

Элементы созданного массива также получат тип float64.

Создать массив, заполненный случайными числами, поможет .random.rand():

>>> In [17]: A = np.random.rand(3, 2) . A Out[17]: array([[0.8206045 , 0.54470809], [0.9490381 , 0.05677859], [0.71148476, 0.4709059 ]])

Говоря точнее, метод .random.rand() возвращает массив с псевдослучайными значениями (от 0 до 1) из множества, сгенерированного по закону равномерного распределения. Обратите внимание, что в отличие от .zeros() и .ones(), .random.rand () на вход принимает не кортеж, а просто два значения через запятую.

Чтобы получить массив с псевдослучайными значениями, взятыми из множества, сгенерированного по закону нормального распределения с нулевым средним и единичной дисперсией, вы можете использовать .random.randn():

>>> In [18]: A = np.random.randn(3, 2) . A Out[18]: array([[-1.20019512, -1.78337814], [-0.22135221, -0.38805899], [ 0.17620202, -2.05176764]])

Почему scipy.linalg лучше, чем numpy.linalg

NumPy имеет встроенный модуль numpy.linalg для решения некоторых задач, связанных с аппаратом линейной алгебры. Обычно scipy.linalg рекомендуют использовать по следующим причинам:

В следующем разделе мы применим scipy.linalg для работы с системами линейных алгебраических уравнений. Наконец-то практика!

Формализация и решение задач с scipy.linalg

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

В этом разделе вы узнаете, как использовать scipy.linalg.solve() для решения СЛАУ. Но прежде чем приступить к работе с кодом, займёмся формализацией задачи и далее рассмотрим простой пример.

Система линейных алгебраических уравнений — это набор из m уравнений, n переменных и вектора свободных членов. Прилагательное «линейных» означает, что все переменные имеют первую степень. Для простоты рассмотрим СЛАУ, где m и n равны 3:

Есть ещё одно требование к «линейности»: коэффициенты K₁ … K₉ и вектор b₁ … b₃ должны быть константами (в математическом смысле этого слова).

В реальных задачах СЛАУ обычно содержат большое количество переменных, что делает невозможным решение систем вручную. К счастью, такие инструменты, как scipy.linalg, могут выполнить эту тяжелую работу.

Задача 1

Мы сначала разберёмся с основами scipy.linalg.solve() на простом примере, а в следующем разделе возьмём задачу посложнее.

Перейдём к матричной записи нашей системы и введём соответствующие обозначения — A, x и b:

Заметьте: левая часть в исходной записи системы — это обычное произведение матрицы A на вектор x.

Всё, теперь можно переходить к программированию.

Пишем код, используя scipy.linalg.solve()

Входными данными для scipy.linalg.solve() будут матрица A и вектор b. Их нужно представить в виде двух массивов: A — массив 2х2 и b — массив 2х1. В этом нам как раз поможет NumPy. Таким образом, мы можем решить систему так:

>>> 1In [1]: import numpy as np 2 . from scipy.linalg import solve 3 4In [2]: A = np.array( 5 . [ 6 . [3, 2], 7 . [2, -1], 8 . ] 9 . ) 10 11In [3]: b = np.array([12, 1]).reshape((2, 1)) 12 13In [4]: x = solve(A, b) 14 . x 15 Out[4]: 16 array([[2.], 17 [3.]])

Разберём приведённый выше код:

Строки 1 и 2: импортируем NumPy и функцию solve() из scipy.linalg.
Строки с 4 по 9: создаём матрицу коэффициентов как двумерный массив с именем A.
Строка 11: создаём вектор-строку b как массив с именем b. Чтобы сделать его вектор-столбцом, вызываем .reshape ((2, 1)).
Строки 13 и 14: вызываем .solve() для нашей системы, результат сохраняем в х и выводим его. Обратите внимание, что .solve() возвращает вектор из чисел с плавающей запятой, даже если все элементы исходных массивов являются целыми числами.

Далее возьмём более сложный пример из реальной практики.

Задача 2: составление плана питания

Это одна из типовых задач, встречающихся на практике: найти пропорции компонентов для получения определенной смеси.

Ниже мы сформируем план питания, смешивая разные продукты, чтобы получить сбалансированную диету.

Нам даны нормы содержания витаминов в пище:

170 единиц витамина А;
180 единиц витамина B;
140 единиц витамина С;
180 единиц витамина D;
350 единиц витамина Е.

Продукт	Витамин A	Витамин B	Витамин C	Витамин D	Витамин E
1	1	10	1	2	2
2	9	1	0	1	1
3	2	2	5	1	2
4	1	1	1	2	13
5	1	1	1	9	2

Необходимо скомбинировать продукты питания так, чтобы их концентрация была оптимальной и соответствовала нормам содержания витаминов в пище.

Обозначим оптимальную концентрацию (количество единиц) для продукта 1 как x1, для продукта 2 — как x2 и так далее. Так как мы будем смешивать продукты, то для каждого витамина (столбца таблицы) можно просто просуммировать значения, по всем продуктам. Учитывая, что сбалансированная диета должна включать 170 единиц витамина А, то, используя данные из столбца «Витамин А», составим уравнение:

Аналогичные уравнения можно составить и для витаминов B, C, D, E, объединив всё в систему:

Запишем полученную СЛАУ в матричной форме:

Теперь для решения системы можно использовать scipy.linalg.solve():

>>> In [1]: import numpy as np . from scipy.linalg import solve In [2]: A = np.array( . [ . [1, 9, 2, 1, 1], . [10, 1, 2, 1, 1], . [1, 0, 5, 1, 1], . [2, 1, 1, 2, 9], . [2, 1, 2, 13, 2], . ] . ) In [3]: b = np.array([170, 180, 140, 180, 350]).reshape((5, 1)) In [4]: x = solve(A, b) . x Out[4]: array([[10.], [10.], [20.], [20.], [10.]])

Мы получили результат. Давайте его расшифруем. Сбалансированная диета должна включать:

10 единиц продукта 1;
10 единиц продукта 2;
20 единиц продукта 3;
20 единиц продукта 4;
10 единиц продукта 5.

Облачные VPS-серверы с быстрыми NVMе-дисками и посуточной оплатой. Загрузка своего ISO.

python
линейные уравнения
математика
линейная алгебра