3a
import numpy as np
def trajectory(horizontal_speed, initial_height, xmax, nmax):
    # Define physical parameters
    g = 980  # cm/s^2
    kappa = 0.5  # 1/s
    a = 1  # cm
    b = 0.25  # cm

    # Define step size
    dx = xmax / (nmax - 1)

    # Initialise arrays
    x = np.linspace(0, xmax, nmax)
    y = np.zeros(nmax)
    y[0] = initial_height
    y[1] = y[0] - (g * (1 - np.exp(-(y[0] - a) / b)) * (dx ** 2)) / (2 * (horizontal_speed ** 2))

    # Iterate to compute the rest of the y values
    for n in range(1, nmax - 1):
        y[n + 1] = (-2 * dx ** 2 * g * (1 - np.exp(-(y[n] - a) / b)) + dx * kappa * horizontal_speed * y[
            n - 1] + 4 * horizontal_speed ** 2 * y[n] - 2 * horizontal_speed ** 2 * y[n - 1]) / (
                               2 * horizontal_speed ** 2 + kappa * dx * horizontal_speed)
    return x, y

3b
def residual(horizontal_speed, initial_height, x_teacup, y_teacup, nmax):
    # Compute the trajectory of the ball using the given parameters
    x, y = trajectory(horizontal_speed, initial_height, x_teacup, nmax)

    # Compute the residual as the difference between the final y-value and the target y-value
    residual = y[-1] - y_teacup

    return residual


3c
import math
from scipy.optimize import bisect

def find_horizontal_speed(initial_height, x_teacup, y_teacup, nmax, bounce=False):
    # Define physical parameters
    g = 980  # cm/s^2

    # Compute the initial guess for u0 based on the parabolic trajectory
    h = initial_height - y_teacup
    u0_guess = x_teacup * math.sqrt(g / (2 * h))

    # Define the bracketing interval for the root-finding algorithm
    lower_bound = 0.8 * u0_guess
    upper_bound = 1.2 * u0_guess

    # Define the residual function F(u0) for a given horizontal speed u0
    def F(u0):
        return residual(u0, initial_height, x_teacup, y_teacup, nmax)

    # Find the root of F(u0) using the bisection method with a tolerance of 0.1cm
    u0_root = bisect(F, lower_bound, upper_bound, xtol=0.1)

    return u0_root

Modify your function find_horizontal_speed() from Task 3c, so that it computes   with exactly 1 bounce of the marble before it lands into the cup when bounce is set to True.
Use a tolerance of no more than 0.1cm.
Strategy
In the bounce case, the interaction with the floor produces unpredictable results. However, we know that friction always reduces the motion, and the interaction with the floor produces most of the damping (also true with real objects).
However, we can pretend there is exactly zero damping, which gives us an estimate that will undershoot the mark when we include friction.
In this case, we estimate that the marble falls to the floor with a parabolic trajectory. It then bounces back and returns to its original height with a mirror-image trajectory. After that, the marble continues to fall toward the target but with a reduced distance.
Putting everything together, produces a lower bound for the speed estimate. The zero-bounce estimate from Task 3c now gives us an upper bound, and we can bracket our solution with:
((xteacup*sqrt(g/(2h)))/(1+2*sqrt(y0/h))) < u0 < xteacup*sqrt(g/(2h)) where h=y0-yteacup
Starting from this interval [umin,umax], it's possible to use a modified version of the bisection method,which also counts the total number of bounces and peaks for the trajectory with the midpoint umidpoint=(umin+umax)/2.
You can then choose the next half-interval based on this information.The constraints are as follows:
 - The marble must bounce exactly once before getting into the cup.
 - The marble must fall downwards into the cup.

下面是测试代码
initial_height = 150
x_teacup = 100
y_teacup = 75
nmax = 2 ** 12
assert round(find_horizontal_speed(initial_height, x_teacup, y_teacup, nmax, bounce=True), 1) == 77.4
print('Test passed.')