First we will consider equation
(4.4.1). Suppose that the beam is a W12x22 structural steel I-beam. Then
\(L = 120\) in.,
\(E = 29 \times 10^{6} \textrm{lb.}/\textrm{in.}^{2}\) and
\(I = 121 \textrm{in.}^{4}\text{.}\) Suppose that the beam is carrying a uniform load of
\(100,000\) lb. so that
\(w = 120,000/120 = 10,000\) and a tension of
\(T = 10,000\) lb.. We calculate from
(4.4.3) \(\alpha = 2.850 \times 10^{-6}\) and
\(\beta = 1.425 \times 10^{-6}\text{.}\) Thus we have the following BVP:
\begin{equation}
y'' = 2.850 \times 10^{-6}y + 1.425 \times 10^{-6}x(120-x), \qquad y(0) = y(120) = 0\text{.}\tag{4.4.5}
\end{equation}
First subdivide the interval \([0,120]\) into four equal subintervals. The nodes of this subdivion are \(x_{0} = 0\text{,}\) \(x_{1} = 30\text{,}\) \(x_{2} = 60\text{,}\) …, \(x_{4} = 120\text{.}\) We will then let \(y_{0}, y_{1}, \ldots, y_{4}\) denote the deflections at the nodes. From the boundary conditions we have immediately:
\begin{equation*}
y_{0} = y_{4} = 0\text{.}
\end{equation*}
To determine the deflections at the interior points we will rely on the differential equation. Recall the central difference formula
\begin{equation*}
y''(x_{i}) \approx \frac{y_{i+1}- 2 y_{i} + y_{i-1}}{h^{2}}\text{.}
\end{equation*}
In this case we have
\(h = (b-a)/n =(120-0)/4= 30\text{.}\) Replacing all the variables in the equation
(4.4.2) by their discrete versions we get
\begin{equation*}
y_{i+1}- 2 y_{i} + y_{i-1}= h^{2} \alpha y_{i} + h^{2} \beta x_{i} (L-x_{i})\text{.}
\end{equation*}
Substituting in for \(\alpha\text{,}\) \(\beta\) and \(h\) we obtain:
\begin{equation*}
\begin{split}y_{i+1}- 2 y_{i} + y_{i-1}&= 900 \times 2.850 \times 10^{-6}y_{i} + 900 \times 1.425 \times 10^{-6}x_{i}(120 - x_{i}) \\&= 2.565 \times 10^{-3}y_{i} + 1.2825 \times 10^{-3}x_{i}(120 - x_{i})\end{split}\text{.}
\end{equation*}
This equation makes sense for \(i = 1,2,3\text{.}\) At \(x_{1} = 30\text{,}\) the equation becomes:
\begin{equation*}
\begin{split}y_{2} - 2y_{1} + y_{0}&= 2.565 \times 10^{-3}y_{1}+ 1.2825 \times 10^{-3}\times 30(90)\\ \Leftrightarrow y_{2} - 2.002565 y_{1}&= 3.46275\end{split}\text{.}
\end{equation*}
Note that this equation is linear in the unknowns \(y_{1}\) and \(y_{2}\text{.}\) At \(x_{2} = 60\) we have:
\begin{equation*}
\begin{split}y_{3} - 2 y_{2} + y_{1}&= .002565 y_{2} + 1.2825 \times 10^{-3}\times 60^{2}\\ \Leftrightarrow y_{3} - 2.002565 y_{2} + y_{1}&= 4.617\end{split}\text{.}
\end{equation*}
At \(x_{3} = 90\) we have (since \(y_{4}=0\))
\begin{equation*}
\begin{split} - 2.002565 y_{3} + y_{2}&= 3.46275\end{split}\text{.}
\end{equation*}
Thus \((y_{1},y_{2},y_{3})\) is the solution of the linear system:
\begin{equation*}
\left( \begin{array}{ccc|c}-2.002565 & 1 & 0 & 3.46275 \\ 1 & -2.002565 & 1 & 4.617 \\ 0 & 1 & -2.002565 & 3.46275\end{array} \right)\text{.}
\end{equation*}
We can easily find the solution of this system in MATLAB:
>> A = [ -2.002565 1 0 ; 1 -2.002565 1 ; 0 1 -2.002565]
>> b = [ 3.46275 4.617 3.46275 ]'
>> y = A
To graph the solution, we need define the
\(x\) values and add on the values at the endpoints:
>> x = 0:30:120
>> y = [0 ; y ; 0]
>> plot(x,y,'d')
Adding a spline will result in an excellent graph.
The exact solution of this BVP is given in
(4.4.4). That equation, with the parameter values for the W12x22 I-beam as in the example, is in the program
myexactbeam.m (
Program A.2.2). We can plot the true solution on the same graph:
>> hold on
>> myexactbeam
Thus our numerical solution is extremely good considering how few subintervals we used and how very large the deflection is.
