Решение кубического полинома с помощью numpy ⇐ Python

[Anonymous]

Я планирую спроектировать фундамент приподнятого резервуара для воды, состоящий из плота, форма которого показана на изображении ниже. Для этого мне нужно определить диаметр плота Dr, решив следующее неравенство:

Код: Выделить всё

𝑁𝑇∙𝐷𝑟 ≥ 8∙𝑀

(где NT и M — входные данные в моем коде)


Основываясь на моих ручных вычислениях, после учета всех необходимых параметров для оценки NT, окончательное уравнение, которое нужно решить, представляет собой кубический многочлен со следующими выражение:

Код: Выделить всё

5.694·Dr³ - 0.942·Dr² +2540.505·Dr - 73612.480 = 0

Я могу легко решить это кубическое полиномиальное уравнение, применив следующий код, который возвращает реальное решение Dr =17,360, соответствующее моим расчетам вручную.

Код: Выделить всё

import sys
from scipy.optimize import fsolve

def equation(Dr):
return 5.694 * Dr**3 - 0.942 * Dr**2 + 2540.505 * Dr - 73612.480

# initial guess for the solver
initial_guess = 0.0

root = fsolve(equation, initial_guess)
OUT = root[0]

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

• Кубическое полиноминальное уравнение: 6,676·Dr³ - 0,942·Dr² +2540,505·Dr - 73612,480 = 0
• Реальное решение: Dr = 16,744 м

Код: Выделить всё

import sys
import os
sys.path.append(r'C:\Users\nono\AppData\Local\Programs\Python\Python38\Lib\site-packages')
import numpy as np
import math

def vertical_load(H, a, b, c, di, de, gb, gs, P0, Dr):
# calculating the entire vertical force at the base
# H is embedment depth
# de: shaft outer diameter
# di: shaft inner diameter
# gb is the concrete density = 2.5T/m3
# gs is the soil density = 1.9T/m3
# P0 vertical load from the superstructure
# Dr raft diameter to find
h = H / 4

d1 = de + 2*a
d2 = d1 + 2*b
d3 = Dr - 2*c

# raft area calculation

s1 = de**2 - di**2
s2 = Dr**2
s3 = (d3**2 - d2**2) / 2
s4 = d2**2
s5 = (d2**2 - d1**2) / 2
s6 = d1**2

# raft weight calculation
Pr = gb * h * math.pi / 4 * (s1 + s2 + s3 + s4 + s5 + s6)

# soil area calculation

t1 = Dr**2 - d3**2
t2 = (d3**2 - d2**2) / 2
t3 = Dr**2 - d2**2
t4 = (d2**2 - d1**2) / 2
t5 = Dr**2 - de**2

# soil weight
Pt = gs * h * math.pi / 4 * (t1 + t2 + t3 + t4 + t5)

# Total vertical load
NT = P0 + Pr + Pt
return NT

def Equi_inequality(Dr, H, a, b, c, di, de, gb, gs, P0, M):
NT = vertical_load(H, a, b, c, di, de, gb, gs, P0, Dr)
return NT * Dr - 8 * M

def polynom(H, a, b, c, di, de, gb, gs, P0, M):
Dr_vals = np.array([0.0, 1.0, 2.0, 3.0])
F_vals = np.array([
Equi_inequality(Dr, H, a, b, c, di, de, gb, gs, P0, M)
for Dr in Dr_vals
])

A = np.array([
[0, 0, 0, 1],
[1, 1, 1, 1],
[8, 4, 2, 1],
[27, 9, 3, 1]
])

coeffs = np.linalg.solve(A, F_vals)
return coeffs

def solve_Dr(H, a, b, c, di, de, gb, gs, P0, M):
a3, a2, a1, a0 = polynom(
H, a, b, c, di, de, gb, gs, P0, M
)

roots = np.roots([a3, a2, a1, a0])

real_roots = [
r.real for r in roots
if abs(r.imag) < 1e-8 and r.real > 0
]

if not real_roots:
raise ValueError("no solution found")

return min(real_roots), (a3, a2, a1, a0)

# Inputs
H  = 4.00
a  = 1.05
b  = 2.00
c  = 1.00
di = 6.80
de = 7.60
gb = 2.5
gs = 1.9
P0 = 2492.52
M  = 9201.56

Dr, coeffs = solve_Dr(H, a, b, c, di, de, gb, gs, P0, M)

print(f"{coeffs[0]:.3f}·Dr³ {coeffs[1]:+.3f}·Dr² "
f"{coeffs[2]:+.3f}·Dr {coeffs[3]:+.3f}")

print(f"Dr = {Dr:.3f} m")

OUT = Dr, coeffs = solve_Dr(H, a, b, c, di, de, gb, gs, P0, M)

Как исправить этот код, чтобы вернуть правильное решение?
Примечание. Я использую механизм Cpython3 в Dynamo для Revit для выполнения кода.
Изменить:
Согласно моим ручным расчетам и для того, чтобы показать, как я пришел к окончательному выражению полиномиального уравнения, Хотелось бы выделить следующие ключевые моменты:

• Код: Выделить всё

  𝑁𝑇∙𝐷𝑟 ≥ 8∙𝑀

  :Неравенство равновесия, которое необходимо проверить, где

  Код: Выделить всё

  M

  — это входные данные, а NT — это общая вертикальная нагрузка, определяемая как NT = P0 + Pr + Pt, где:

• Код: Выделить всё

  P0

  — входное значение в коде

• Код: Выделить всё

  Pr

  : вес плота рассчитывается в соответствии с геометрией, показанной на изображении.

• Код: Выделить всё

  Pt

  : вес почвы, рассчитанный в соответствии с геометрией, показанной на изображении.

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



Подробнее здесь: https://stackoverflow.com/questions/798 ... sing-numpy