Python - Solving Bernoulli's Beam Equation With Scipy

The process of answering the question has already started in the question on the link bellow, but that topic was specifically about integrating a function, which was answered. So I

Solution 1:

You are integrating a differential equation, your approach of computing in a loop the definite integrals is, let's say, sub-optimal.

The standard approach in Scipy is the use of scipy.integrate.solve_ivp, that uses a suitable integration method (by default, Runge-Kutta 45) to provide the solution in terms of a special object.

As usual in the field of numerical integration of ordinary differential equations, the method is limited to a system of 1st order differential equation, but your 2nd degree equation can be transformed to a system of 1st degree equations introducing an helper function

    Y" = M ⇒ {y' = h, h' = M} 

While this sounds complicated, its implementation is quite simple

In [51]: #########################################################################
    ...: # L, EJ = 1.0
    ...: #########################################################################
    ...: # exact solution
    ...: from numpy import linspace
    ...: x = linspace(0, 1, 51)
    ...: y1, y0 = (x**2-1)/2, (x**3-3*x+2)/6
    ...: #########################################################################
    ...: # numerical solution
    ...: from scipy.integrate import solve_ivp
    ...: s = solve_ivp(
    ...:     lambda x, Y: [Y[1], x], # Y[0]=y0, Y[1]=y1, dy0dx=y1, dy1dx=x
    ...:     [1.0, 0.0], # interval of integration (NB: reversed, because...)
    ...:     [0.0, 0.0], # initial conditions (at the 1st point of integ interval)
    ...:     t_eval=np.linspace(1, 0, 101) # where we want the solution to be known
    ...:     )
    ...: #########################################################################
    ...: # plotting
    ...: from matplotlib.pyplot import grid, legend, plot
    ...: plot(x, y0, label="Exact y")
    ...: plot(x, y1, label="Exact y'")
    ...: plot(s.t[::2], s.y[0][::2], label="Numeric y",  linestyle='', marker='.')
    ...: plot(s.t[::2], s.y[1][::2], label="Numeric y'", linestyle='', marker='.')
    ...: legend() ; grid() ;

In [52]: 

---

---

The OP reported an issue understanding solve_ivp(lambda x, Y: [Y[1], x], ....

We have a system of 1st order ODEs in normal form

y₁' = f₁(x, y₁(x), …, yₙ(x))
…   = …
yₙ' = f₁(x, y₁(x), …, yₙ(x))

that can be written, using capital letters to signify vector quantities

Y' = F(x, Y(x))

to solve the system of differential equations solve_ipv needs exactly this F(x, Y) function.

In place of the lambda expression one could write a function definition like the following, that is possibly more self explanatory

defF(x, Y):
    # unpack the vector of function values
    y = Y[0]
    h = Y[1]
    # compute the derivatives
    dy_over_dx = h
    dh_over_dx = x
    # return the vector of derivatives valuesreturn [dy_over_dx, dh_over_dx]

s = solve_ivp(F, …)

that in the answer was succinctly (too much succinctly?) was expressed as lambda x,Y:[Y[1],x] …