Блог

Я пытаюсь решить следующее уравнение для r:

 from sympy import pi, S, solve, solveset, nsolve, symbols

(n_go, P_l, T, gamma_w, P_g, r, R_mol) = symbols(
    'n_go, P_l, T, gamma_w, P_g, r, R_mol', real=True)

expr = -P_g   P_l - 3*R_mol*T*n_go/(4*r**3*pi)   2*gamma_w/r
soln = solveset(expr, r, domain=S.Reals)
soln1 = solve(expr, r)
  

soln имеет форму Complement(Intersection(FiniteSet(...))) , с которой я действительно не знаю, что делать.
soln1 представляет собой список из 3 выражений, два из которых являются сложными. На самом деле, если я подставлю значения для символов и вычислю решения для soln1, все они будут сложными:

 vdict = {n_go: 1e-09, P_l: 101325, T: 300, gamma_w: 0.07168596252716256, P_g: 3534.48011713030, R_mol: 8.31451457896800}
for result in soln1:
    print(result.subs(vdict).n())
  

ВОЗВРАТ:

 -9.17942953565355e-5   0.000158143657514283*I
-9.17942953565355e-5 - 0.000158143657514283*I
0.000182122477993494   1.23259516440783e-32*I
  

Интересно, что сначала подстановка значений, а затем использование solveset() или solve() дает реальный результат:

 solveset(expr.subs(vdict), r, domain=S.Reals).n()
  

{0.000182122477993494}

И наоборот, nsolve завершается ошибкой с этим уравнением, если только начальная точка не содержит первые 7 значащих цифр решения (!):

 nsolve(expr.subs(vdict), r,0.000182122 )
  

Ошибка значения: не удалось найти корень в пределах заданного допуска. (9562985778.9619347103> 2.16840434497100886801e-19)

Это не должно быть так сложно, вот сюжет:

Мои вопросы:

1. Почему здесь nsolve так бесполезно?
2. Как я могу использовать решение, возвращенное из solveset , для вычисления любых численных решений?
3. Почему я не могу получить реальное решение, solve если я сначала решаю, а затем подставляю значения?

Ответ от Maelstrom хороший, но я просто хочу добавить несколько моментов.

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

Кроме того, перед заменой символов solve не может узнать, какой из трех корней является реальным. Вот почему вы получаете три решения solve(expr, r) и только одно решение solve(expr.subs(vdict), r) . Третье решение, которое является реальным после подстановки, такое же (игнорируя крошечную мнимую часть), возвращаемое solve после подстановки:

 In [7]: soln1[2].subs(vdict).n()                                                                                                                              
Out[7]: 0.000182122477993494   1.23259516440783e-32⋅ⅈ

In [8]: solve(expr.subs(vdict), r)                                                                                                                            
Out[8]: [0.000182122477993494]
  

Поскольку многочлен плохо обусловлен и имеет большой градиент в корне nsolve , ему трудно найти этот корень. Однако nsolve можно найти корень, если задан достаточно узкий интервал:

 In [9]: nsolve(expr.subs(vdict), r, [0.0001821, 0.0001823])                                                                                                   
Out[9]: 0.000182122477993494
  

Поскольку это, по сути, многочлен, вам лучше всего преобразовать его в многочлен и использовать nroots . Самый быстрый способ сделать это — использовать as_numer_denom , хотя в этом случае это вводит ложный корень в ноль:

 In [26]: Poly(expr.subs(vdict).as_numer_denom()[0], r).nroots()                                                                                               
Out[26]: [0, 0.000182122477993494, -9.17942953565356e-5 - 0.000158143657514284⋅ⅈ, -9.17942953565356e-5   0.000158143657514284⋅ⅈ]
  

Ваше выражение по сути является кубическим уравнением.

Применение subs до или после решения не должно существенно ничего менять.

soln

soln имеет вид Complement(Intersection(FiniteSet(<3 cubic solutions>), Reals), FiniteSet(0)) , т.е. кубическое решение в вещественной области, исключая 0.

Следующее должно дать вам простой FiniteSet ответ, но evalf , похоже, не очень хорошо реализовано для наборов.

 print(soln.subs(vdict).evalf())
  

Надеюсь, что-то будет сделано с этим в ближайшее время.

1

Причина nsolve , по которой это бесполезно, заключается в том, что график почти асимптотически вертикальный. Согласно вашему графику, градиент составляет примерно 1.0e8. Я не думаю nsolve , что это полезно для таких крутых графиков.

Выводя ваше замещенное выражение, мы получаем: уменьшая масштаб, мы получаем:

Это довольно дикая функция, и я подозреваю nsolve , что использует эпсилон, который слишком велик, чтобы быть полезным в этой ситуации. Чтобы исправить это, вы могли бы предоставить более разумные числа, которые при замене ближе к 1. (Рассмотрите возможность предоставления разных единиц измерения. например. вместо метров / год рассмотрим км / час)

2

Трудно сказать вам, как обращаться с выводом solveset в целом, потому что с каждым типом набора нужно обращаться по-разному. Это также не имеет математического смысла, поскольку soln.args[0].args[0].args[0] должно давать первое кубическое решение, но оно забывает, что оно должно быть реальным и отличным от нуля.

Вы можете использовать args или preorder_traversal или вещи для навигации по дереву. Также должно помочь чтение документации по различным наборам. solve и solveset должны использоваться «в интерактивном режиме», потому что существует множество возможных результатов с множеством способов их понимания.

3

Я считаю soln1 , что имеет 3 решения вместо 4, как вы заявляете. В противном случае ваш цикл будет печатать 4 строки вместо 3. Все они технически сложны (как и природа с плавающими числами в Python). Однако третье решение, которое вы предоставляете, имеет очень маленький воображаемый компонент. Чтобы удалить подобные привередливые вещи, есть аргумент, chop который должен помочь:

 for result in soln1:
    print(result.subs(vdict).n(chop=True))
  

Один из результатов 0.000182122477993494 — это то, что выглядит как ваш корень.

Вот ответ на основной вопрос: как эффективно вычислить корни приведенного выше уравнения?

Основываясь на предложении @OscarBenjamin, мы можем работать еще лучше и быстрее, используя Poly и roots вместо nroots . Ниже sympy мгновенно вычисляет корни уравнения для 100 различных значений P_g, сохраняя при этом все остальное постоянным:

 from sympy import pi, Poly, roots, solve, solveset, nsolve, nroots, symbols
(n_go, P_l, T, gamma_w, P_g, r, R_mol) = symbols(
    'n_go, P_l, T, gamma_w, P_g, r, R_mol', real=True)
vdict = {pi:pi.n(), n_go:1e-09, P_l:101325, T:300, gamma_w:0.0717, R_mol: 8.31451457896800}
expr = -P_g   P_l - 3*R_mol*T*n_go/(4*r**3*pi)   2*gamma_w/r
expr_poly = Poly(expr.as_numer_denom()[0], n_go, P_l, T, gamma_w, P_g, r, R_mol, domain='RR[pi]')
result = [roots(expr_poly.subs(vdict).subs(P_g, val)).keys() for val in range(4000,4100)]
  

Все, что остается, это проверить, соответствуют ли решения нашим условиям (положительным, реальным). Спасибо, @OscarBenjamin!

PS: Должен ли я расширить тему выше, чтобы включить nroots и корни?