I am self-studying numerical methods and this is a problem from the book I am using:
Let $x$ be the distance along the axis of the beam with $x= 0$ at the built-in end. The variable $y$ measures the deflection of the beam. It can be shown that the elastic deflections satisfy the differential equation \begin{align*} y′′=\frac{F(L−x)[1 + (y′)^{2}]^{3/2}}{EI} \end{align*}
where $E$ is Young’s modulus for the material and $I$ is the moment of inertia of a cross section of the beam. Since both the deflection and its slope vanish at $x= 0$, we have the initial conditions $y(0) = 0, y′(0) = 0$. Use a numerical method of your choice to approximate $y(100)$ using $h= 2$ for the case of a beam of high-strength aluminum alloy with rectangular cross section. The beam is 100 in.long and $I= 0.128 in^4$. The value of Young’s modulus is $E= 1.0×10^7$ psi. Solve the IVP for $F= 60,80,100,150$,and $200$ lbs to obtain some feeling about the deflection of the beam.
I wrote the given IVP as a system: I let $x=y'$, so that $x'=y''=\frac{F(L−x)[1 + x^{2}]^{3/2}}{EI}$, yielding the system (after plugging in the given values for $E, I, L$) \begin{align*} &y'=x, \qquad y(0)=0\\ &x'=\frac{F(100−x)[1 + x^{2}]^{3/2}}{10^{7}\cdot 0.128}, \qquad x(0)=0. \end{align*} I wrote a Runge-Kutta (order 4) code to solve the IVP for the given values of $F$. The output seems reasonable; since as the $x$ increases (i.e., the force is applied further along on the beam), the deflection increases, which is intuitive. Similarly, the deflection increases with force. I wanted to just double-check that my system is correct, etc. (and that my code seems fine; I'm new to it). The first few iterations of my MATLAB code with F=60:
x y
0.00 0.00000
2.00 0.03882
4.00 0.15439
6.00 0.34566
8.00 0.61183
10.00 0.95240
12.00 1.36709
14.00 1.85587
16.00 2.41895
18.00 3.05674
20.00 3.76988
22.00 4.55920
24.00 5.42580
26.00 6.37095
28.00 7.39620
30.00 8.50331
If someone could double-check that this seems correct, I would appreciate it.