Блог

(TLDR: функция linsolve SymPy не может решить систему линейных уравнений, полученную в результате применения метода конечных разностей к BVP ODE, когда вы передаете уравнения в виде чистого списка Python, но можете сделать это, поместив указанный список уравнений в матричную функцию SymPy. Это может быть ошибка, которую необходимо исправить, особенно учитывая, что в документации SymPy в примере, который они вам приводят, они передают список в качестве аргумента linsolve. )

У меня есть обыкновенное дифференциальное уравнение краевой задачи, которое я намерен решить с помощью метода конечных разностей. Моя ОДА, в представлении сочувствия, такова x*y(x).diff(x,2) y(x).diff(x) 500 = 0 :с y(1)=600 И. y(3.5)=25 Весь код выглядит следующим образом:

 import sympy as sp
from numpy import *
from matplotlib.pyplot import *
y = sp.Function('y')
ti = time.time()
steps = 10
ys = [ y(i) for i in range(steps 1) ]
ys[0], ys[-1] = 600, 25
xs = linspace(1, 3.5, steps 1)
dx = xs[1] - xs[0]
eqns = [ xs[i]*(ys[i-1] - 2*ys[i]   ys[i 1])/dx**2   
        ( ys[i 1] - ys[i-1] )/2/dx   500
        for i in range(1, steps) ]
ys[1:-1] = sp.linsolve(eqns, ys[1:-1]).args[0]
scatter(xs, ys)
tf = time.time()
print(f'Time elapsed: {tf-ti}')
 

Для десяти шагов это работает просто отлично. Однако, если я сразу же поднимусь выше, например, на 11 шагов, SymPy больше не сможет решить систему линейных уравнений. Попытка построить график результатов выбрасывает a TypeError: can't convert expression to float . Изучение списка значений ys показывает, что одна из переменных, в частности предпоследняя, не решалась для linsolve sympy, но вместо этого решения для других переменных решаются в терминах этой нерешенной переменной. Например, если у меня есть 50 шагов, предпоследняя переменная будет y(49) и это показано на решениях других неизвестных и не решенных, например 2.02061855670103*y(49) - 26.1340206185567 . Напротив, в другой ОДЕ BVP, которую я решил y(x).diff(x,2) y(x).diff(x) y(x) - 1 с y(0)=1.5 помощью and y(3)=2.5 , нет никаких проблем, хочу ли я 10, 50 или 200 шагов. Это прекрасно решает все переменные, но это, похоже, своеобразное исключение, поскольку я столкнулся с вышеупомянутой проблемой со многими другими ОДАМИ.

Непоследовательность Симпатии здесь была довольно неприятной. Единственное утешение в том, что до того, как я столкнулся с этой проблемой, я уже несколько раз решал ее с помощью разного количества шагов, которые я хотел. Я заключил eqns переменную в матричную функцию sympy, как в ys[1:-1] = sp.linsolve(sp.Matrix(eqns), ys[1:-1]).args[0] просто потому, что она лучше отображалась в терминале. Но для решения в файле сценария я подумал, что обертывать его внутри sp.Matrix не нужно, и я, естественно, удалил его, чтобы упростить ситуацию.

При форматировании вопроса для SO (или где-либо еще) вежливо предоставлять полный пример кода, не пропуская импорт и т. Д. Также вам следует свести проблему к минимуму и удалить все ненужные детали. Имея это в виду, лучший способ продемонстрировать проблему-это на самом деле выяснить, каковы аргументы linsolve , и представить их напрямую, например:

 from sympy import *

y = Function('y')

eqns = [
    -47.52*y(1)   25.96*y(2)   13436.0,
    25.96*y(1) - 56.32*y(2)   30.36*y(3)   500,
    30.36*y(2) - 65.12*y(3)   34.76*y(4)   500,
    34.76*y(3) - 73.92*y(4)   39.16*y(5)   500,
    39.16*y(4) - 82.72*y(5)   43.56*y(6)   500,
    43.56*y(5) - 91.52*y(6)   47.96*y(7)   500,
    47.96*y(6) - 100.32*y(7)   52.36*y(8)   500,
    52.36*y(7) - 109.12*y(8)   56.76*y(9)   500,
    56.76*y(8) - 117.92*y(9)   61.16*y(10)   500,
    61.16*y(9) - 126.72*y(10)   2139.0
]

syms = [y(1), y(2), y(3), y(4), y(5), y(6), y(7), y(8), y(9), y(10)]

print(linsolve(eqns, syms))
 

Здесь вы надеялись получить простое численное решение для каждой из неизвестных, но вместо этого был получен результат (из SymPy 1.8)::

 FiniteSet((5.88050359812056*y(10) - 5.77315239260531, 10.7643116711359*y(10) - 528.13328974178, 14.9403214726998*y(10) - 991.258097567359, 9.85496358613721e 15*y(10) - 1.00932650309452e 18, 7.35110502818395*y(10) - 312.312287998229, 5.84605452313345*y(10) - 217.293922525318, 4.47908204606922*y(10) - 141.418192750506, 3.22698120573309*y(10) - 81.4678489766327, 2.07194244604317*y(10) - 34.9738391105298, 1.0*y(10)))
 

linsolve Функция вернет решение, включающее одно или несколько неизвестных, если в системе нет уникального решения. Обратите также внимание, что есть некоторые большие числа, подобные 9.85496358613721e 15 которым, предполагают, что здесь могут быть численные проблемы.

На самом деле это ошибка в SymPy, и она уже исправлена в главной ветке: https://github.com/sympy/sympy/pull/21527

Если вы установите SymPy из git, то вместо этого вы можете найти следующее в качестве вывода:

 FiniteSet((596.496767861074, 574.326903264955, 538.901024315178, 493.575084012669, 440.573815245681, 381.447789421181, 317.320815173574, 249.033388036155, 177.23053946471, 102.418085492911))
 

Также обратите внимание, что, как правило, лучше избегать использования поплавков в SymPy, поскольку это библиотека, предназначенная для точных символьных вычислений. Решение системы уравнений с плавающей запятой, подобной этой, может быть выполнено гораздо эффективнее с помощью NumPy или какой-либо другой библиотеки с плавающей запятой фиксированной точности, которая может использовать процедуры BLAS/LAPACK. Чтобы использовать рациональную арифметику с симпатией здесь, вам просто нужно изменить свою linspace строку на

 xs = [sp.Rational(x) for x in linspace(1, 3.5, steps 1)]
 

который затем будет отлично работать с SymPy 1.8 и на самом деле быстрее (по крайней мере, если у вас установлен gmpy2).