I'm fairly new to programming and am trying to write a program using the inbuilt function 'integrate.solve_bvp' to determinine the trajectory of a projectile subject to boundary conditions.
I'm not by any means a programmer, so my knowledge and understanding is extremely limited. Please explain like I'm 5.
I need to be able to determine the launch and final velocity of a projectile, given its launch angle, and the time taken for it to return to the ground while considering drag.
I've written some code, but it doesn't work and I don't know why. I've tried reading the documentation, but it all goes over my head.
I started by considering the 1-D case (launch angle of 90 degrees), which seemed failly simple. I used SUVAT to gain a guess of the launch velocity (ie, the launch velocity if ignoring drag) and wrote the following code:
import numpy as np
from scipy import integrate
import matplotlib.pyplot as plt
drag_coef = 0.47 # average drag coefficient of a golf ball
air_density = 1.293
area = 1.45*10**-3 # cross sectionl area of a golf ball
mass = 45.9*10**-3 # mass of a golf ball
g = -9.80665
def function(time,height):
drag_factor = drag_coef * air_density * area / (2*mass)
return height[1], (g-drag_factor*height[1]*np.abs(height[1]))
def boundary_conditions(height_0,height_end):
return height_0[0], height_end[0]
time_scale = 10 # time at which you want projectile to hit the ground
velocity_0_guess = 49 # initial velocity guess (t=0)
time = np.linspace(0, time_scale, time_scale*1000+1)
height_0 = np.zeros(len(time))
velocity_0 = velocity_0_guess * np.ones(len(time))
height = np.array((height_0,velocity_0))
res = integrate.solve_bvp(function, boundary_conditions, time, height, max_nodes = time_scale*1000+1)
print(res.y[1] [0]) # calculated V_0
print(res.y[1] [time_scale*1000]) # Calculated V_end
print(res)
plt.plot(time, res.y[0], label="S_z") # caculated S_z
plt.xlabel("time [s]")
plt.ylabel("displacement [m]")
plt.show()
plt.plot(time, res.y[1], label="V_z") # calcuted V_z
plt.xlabel("time [s]")
plt.ylabel("velocity [m/s]")
plt.show()
However, even for this seemingly simple case, when I "print(res)"; I get a success result of 'False', along with the following statement:
message: 'The maximum number of mesh nodes is exceeded.'
And I don't know why as I think I've defined the number of nodes to equal the number of points in time that are being considered.
Clearly however this isn't the case as when I halve my 'time_scale' and 'velocity_0_guess' to 5 and 24.5 respecitvely, I get a successful result, even though both of these should equally valid:
message: 'The algorithm converged to the desired accuracy.'
I've tried to google the issue, but I haven't found anything that's been able to help me. I've looked through Reddit, and StackOverflow with no success. And I've even tried using ChatGPT to help fix my code, but that too was to no avail. So my step is posting a question.
I don't know if this is relevant, but I've been writing this program via the website: repl.it
The solver works iteratively. It alternates between computing an approximate solution of the DE on a fixed time grid and refining the time grid using an error estimate along the approximate solution.
The general idea is to start with an initial guess that roughly gives the shape of the desired solution. This usually should have no more points than necessary to define the shape. The solver then fills in the grid.
The resulting solution is found either in the function table res.x, res.y or with the "dense output" interpolation res.sol.
For your constant guess it is completely sufficient to have a minimal grid of two points
time = [0,time_scale]
height = [[0.0]*2, [velocity_0_guess]*2]
This finishes without complaint, giving res.y[1,[0,-1]] = [ 95.93681148, -30.32436139] and the number of grid points as len(res.x) = 27. Visibly, the time grid is no longer that of time, so you need to use res.x in the plots.
You can get a denser grid and more accurate solution by setting the error tolerance tol lower than its default value 1e-3
res = integrate.solve_bvp(function, boundary_conditions, time, height, tol=1e-6);
giving len(res.x) = 186 and res.y[1,[0,-1]] = [ 95.93702666, -30.32440457]