[SOLVED] Solve second order PDEs with `scipy.integrate`, with both initial and final position known

Solve second order PDEs with `scipy.integrate`, with both initial and final position known

I'm confused reading scipy.integrate.solve_ivp documentation.

I'm interested in a ballistic problem with drag and Magnus effect, but I'm focussing first on the simpler problem, considering only gravitational force. The corresponding PDE is

Transforming to a first order PDE, we can write

I have the initial and final 3D positions and time of the ball $\big(x(t),y(t),z(t)\big)$ for $t\in[T_0,T_1]$ , but I don't understand how to provide this information to the solver. It expects y0 which, in my notations, is $X(T_0)$ , but I don't know the velocity at $T_0$ . (note: I know I can infer it from the second degree solution, but I don't want to since the solution will get very much more complex once I integrate the other forces).

How should the problem be transformed to add the other initial condition on the position $X(T_1)$ ?

note: I also looked at the documentation of solve_bvp, but from my understanding, it doesn't fit the problem I try to solve…

Solution

The function scipy.integrate.solve_bvp is indeed appropriate, with p0 and p1: (x(t), y(t), z(t)) at t=T0, T1. For posterity, here are the two functions required by solve_bvp:

def bc(X0, X1, args=None):
    x0, y0, z0, vx0, vy0, vz0 = X0
    x1, y1, z1, vx1, vy1, vz1 = X1

    return np.array([
        x0 - p0.x,
        y0 - p0.y,
        z0 - p0.z,
        x1 - p1.x,
        y1 - p1.y,
        z1 - p1.z,
    ])

def pde(t, X, args=None):
    x, y, z, vx, vy, vz = X
    dXdt = np.vstack([
        vx,
        vy,
        vz,
        np.ones_like(t)*0,
        np.ones_like(t)*0,
        np.ones_like(t)*g,
    ])
    return dXdt