Как решать сложные уравнения численно в Python? ⇐ Python

[Anonymous]

Изначально у меня было следующее уравнение:
2mgv×sin(\alpha) = CdA×\rho(v^2 + v_{wind}^2 + 2vv_{wind}cos(\phi))^(3/2)


которое я мог бы выразить в виде следующего нелинейного уравнения:
(K × v)^(2/3) = v^2 + v_{wind} + 2vv_{wind}cos(\phi)

Чтобы решить эту проблему, мне нужно использовать численный подход.
Я попробовал написать код для этого на Python, используя fsolve из scipy.optimize. Однако результаты, которые я получил, не слишком многообещающие. Что еще мне следует попробовать? Должен ли я использовать другой подход/пакет или просто нужно улучшить мой код? Я также заметил, что результат сильно зависит от v_initial_guess.
Обратите внимание, что я считаю себя новичком в программировании.
Я тоже попробовал написать код для решения уравнения для v, используя метод Ньютона-Рафсона, но мне это не удалось.
Вот мой код:
import numpy as np
from scipy.optimize import fsolve

m = 80
g = 9.81
alpha = np.radians(10) #incline
CdA = 0.65
rho = 1.225
v_w = 10
phi = np.radians(30) #wind angle with the direction of motion

sin_alpha = np.sin(alpha)
cos_phi = np.cos(phi)

def equation(v):
K = ((m * g * sin_alpha) / ((CdA * rho))**(2/3))
return K * v**(2/3) - v**2 - 2*v*v_w*cos_phi - v_w**2

v_initial_guess = 30

v_solution = fsolve(equation, v_initial_guess, xtol=1e-3)

print("v:", v_solution[0])type here

РЕДАКТИРОВАТЬ:
Вот как сейчас выглядит мой код:
import numpy as np
from scipy.optimize import fsolve
import matplotlib.pyplot as plt

m = 80
g = 9.81
alpha = np.radians(2) # incline
CdA = 0.321
rho = 1.22
v_w = 5
phi = np.radians(180) # wind angle with the direction of motion

sin_alpha = np.sin(alpha)
cos_phi = np.cos(phi)

def lhs(v):
return m * g * v * sin_alpha

def rhs(v):
return 0.5 * CdA * rho * (v**2 + v_w**2 + 2*v*v_w*cos_phi)**(3)

def difference(v):
return lhs(v) - rhs(v)

# fsolve to find the intersection
v_initial_guess = 8
v_intersection = fsolve(difference, v_initial_guess)[0]

v_values = np.linspace(0.1, 50, 500)

lhs_values = lhs(v_values)
rhs_values = rhs(v_values)

plt.figure(figsize=(10, 6))
plt.plot(v_values, lhs_values, label='$2mgv\\sin(\\alpha)$', color='blue')
plt.plot(v_values, rhs_values, label='$CdA\\rho(v^2 + v_{wind}^2 + 2vv_{wind}\\cos(\\phi))^{3/2}$', color='red')
plt.xlabel('Velocity (v)')
plt.xlim(0, 20)
plt.title('LHS and RHS vs. Velocity')
plt.legend()
plt.grid(True)
plt.ylim(0, 2000)
plt.show()

print(f"The intersection occurs at v = {v_intersection:.2f} m/s")

P_grav_check = m *g * sin_alpha * v_intersection
P_air_check = 0.5 * CdA * rho * (v_intersection ** 2 + v_w ** 2 + 2 * v_intersection * v_w * cos_phi) ** (3)

print(P_grav_check)
print(P_air_check)


Подробнее здесь: https://stackoverflow.com/questions/786 ... -in-python