I'm a professor of Modeling & Simulation for engineering undergraduates and one of the examples I used to teach while talking about boundary value problems is the following:
A ball is thrown from a height of 1.5 metres. The initial speed of the ball is unknown, but it reaches its maximum height at 0.5 seconds. How fast should we throw the ball?
We could easily use the algebraic equation y(t) = y0 + v0t - gt²/2 to solve this problem, but that is supposed to be an introductory problem that will lead us to more interesting boundary-value differential equations further on.
I've written a version of the following program a couple years ago and it used to run fine while using fsolve to find the unknown initial condition (velocity of the ball at t=0).
import numpy as np
from scipy.integrate import odeint
from scipy.optimize import fsolve, bisect
def func(Y,t):
# Parameters:
g = 9.8 # m/s^2
# ODE
y, dy = Y # position, velocity
d2y = -g # acceleration = - gravity
return dy, d2y # velocity, acceleration
def shooting(dy0,t): # dy0 is a guess for the initial condition that
Y0 = [1.5,dy0] # ... also satisfies the boundary condition
Y = odeint(func, Y0, t)
#y = Y[:,0]
dy = Y[:,1]
zero = dy[-1] # the boundary condition is: maximum height at the final time
return zero
tf = 0.5 # s
t = np.linspace(0,tf)
dy0 = fsolve(shooting, 0, args=(t,))
#dy0 = bisect( shooting, 0, 10, args = (t,) )
print(f"Initial velocity of the ball: v0 = {dy0:.2f} m/s")
print(f"Final result of function shooting: zero = {shooting(dy0, t)}")
However, I've run it recently and it started to result in a ValueError.
ValueError: setting an array element with a sequence. The requested array has an inhomogeneous shape after 1 dimensions. The detected shape was (2,) + inhomogeneous part.
If I use the bisection method to find the unknown, the problem works just fine, though. Can anyone explain me what causes the error?
