Lösen eines kubischen Polynoms mit Numpy ⇐ Python

[Anonymous]

Ich beabsichtige, das Fundament eines Hochwassertanks zu entwerfen, bestehend aus einem Floß mit der im Bild unten gezeigten Form. Dazu muss ich den Floßdurchmesser Dr bestimmen, indem ich die folgende Ungleichung löse:

Code: Select all

𝑁𝑇∙𝐷𝑟 ≥ 8∙𝑀

(wobei NT und M Eingaben in meinem Code sind)


Basierend auf meinen Handberechnungen und nach Berücksichtigung aller notwendigen Parameter zur Auswertung von NT ist die endgültige zu lösende Gleichung ein kubisches Polynom mit dem folgenden Ausdruck:

Code: Select all

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

Ich kann diese kubische Polynomgleichung leicht lösen, indem ich den folgenden Code anwende, der die tatsächliche Lösung Dr =17,360 zurückgibt, die meiner Handberechnung entspricht.

Code: Select all

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]

Ich möchte jedoch die Polynomkoeffizienten automatisch bestimmen, anstatt sie manuell zu berechnen. Dies würde es mir ermöglichen, den Code später wiederzuverwenden, um ein weiteres Floß basierend auf anderen Eingaben zu entwerfen, ohne diese Koeffizienten manuell neu berechnen zu müssen. Dafür habe ich den folgenden Code angewendet, der, wie Sie sehen können, falsche Ergebnisse zurückgibt:

• Kubische Polynomgleichung: 6,676·Dr³ - 0,942·Dr² +2540,505·Dr - 73612.480 = 0
• Echte Lösung: Dr = 16,744 m

Code: Select all

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)

Wie kann ich diesen Code korrigieren, um die richtige Lösung zurückzugeben?
Anmerkung: Ich verwende die Cpython3-Engine in Dynamo für Revit, um den Code auszuführen.
Bearbeiten:
Nach meinen Handberechnungen und um zu zeigen, wie ich zum endgültigen Ausdruck der Polynomgleichung gekommen bin, Ich möchte die folgenden Kernpunkte hervorheben:

• Code: Select all

  𝑁𝑇∙𝐷𝑟 ≥ 8∙𝑀

  :Überprüfen Sie die Gleichgewichtsungleichung wo

  Code: Select all

  M

  ist eine Eingabe und NT ist die gesamte vertikale Last, definiert als NT = P0 + Pr + Pt wobei:

• Code: Select all

  P0

  ist ein Eingabewert im Code

• Code: Select all

  Pr

  : Floßgewicht, berechnet gemäß der im Bild gezeigten Geometrie

• Code: Select all

  Pt

  : Bodengewicht, berechnet gemäß der im Bild gezeigten Geometrie.

Durch Berücksichtigung aller Eingabedaten und Anwendung der Gleichgewichtsungleichung erhalten wir die kubische Polynomgleichung, die aus meinen Handberechnungen abgeleitet wurde und in diesem Bild gezeigt wird: