How to solve complex equations numerically in Python?

Question

Originally, I had the following equation:

2mgv×sin(\alpha) = CdA×\rho(v^2 + v_{wind}^2 + 2vv_{wind}cos(\phi))^(3/2)

which I could express as the following non-linear equation:

(K × v)^(2/3) = v^2 + v_{wind} + 2vv_{wind}cos(\phi)

To solve this, I need to use a numerical approach.

I tried writing a code for this in Python using fsolve from scipy.optimize. However, the results I got are not too promising. What else should I try? Should I use a different approach/package or just my code needs to be improved? I also experienced that the result is highly dependent on v_initial_guess.

Please note, that I would consider myself as a beginner in programming.

I also tried writing a code to solve the equation for v using the Newton-Raphson method but I wasnt too successful.

Here is my code:

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

EDIT: This is what my code looks like now:

    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)

Answer 1

The solution, using the secant method is the following:

import numpy as np
import matplotlib.pyplot as plt


m = 60  
g = 9.81  
alpha = np.radians(2)  # incline 
CdA = 0.321  
rho = 1.225
v_wind = 4
phi = np.radians(60)  # wind angle with the direction of motion
mu = 0.005

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


lhs = m * g * sin_alpha     # gravitational force

def rhs(v):
    return 0.5 * CdA * rho * (v + v_wind *cos_phi)**2    # drag force

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

def secant_method(func, v0, v1, tolerance=1e-6, max_iterations=1000):
    for _ in range(max_iterations):
        f_v0 = func(v0)
        f_v1 = func(v1)
        
        v_new = v1 - f_v1 * (v1 - v0) / (f_v1 - f_v0)
        if abs(v_new - v1) < tolerance:
            return v_new
        v0, v1 = v1, v_new
    raise ValueError("Secant method did not converge")

v_initial_guess_1 = 15
v_initial_guess_2 = 25
v_intersection = secant_method(difference, v_initial_guess_1, v_initial_guess_2)

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


P_grav_check = m * g * sin_alpha * v_intersection
P_air_check = 0.5 * CdA * rho * (v_intersection + v_wind * cos_phi)**2 * v_intersection
P_friction = mu * m * g * np.cos(alpha) * v_intersection

alpha_degree = np.degrees(alpha)
alpha_percentage = np.tan(alpha) *100

print(f"\nopposing component of v_wind: ",v_wind *cos_phi, "m/s")
print(f"The intersection occurs at v = {v_intersection:.4f} m/s")
print(f" at a {alpha_degree} degrees incline which equals to {alpha_percentage:.2f} %")
print(f"\nP_grav: {P_grav_check:.2f} W")
print(f"P_air_check: {P_air_check:.2f} W")
print(f"P_friction:  {P_friction:.2f} W")
print(f"P = {P_grav_check + P_air_check + P_friction}")


plt.figure(figsize=(10, 6))
plt.plot(v_values, [lhs]*len(v_values), label='$mg\\sin(\\alpha)$', color='blue')
plt.plot(v_values, rhs_values, label='$0.5CdA\\rho ({v + v_{wind}\\cos(\\phi)})^2$', color='red')
plt.xlabel('Velocity (v)')
plt.xlim(0, 20)
plt.title('Gravitational vs. Drag Force')
plt.legend()
plt.grid(True)
plt.ylim(0, 200)
plt.show()