problem description
solution
The first step is to convert the system to state space.\[ y^{\prime \prime \prime }+10y^{\prime \prime }-5y^{3}+ty=t^{2}\] Let \(x_{1}=y,x_{2}=y^{\prime },x_{3}=y^{\prime \prime }\), therefore\begin{align*} \dot{x}_{1} & =x_{2}\\ \dot{x}_{2} & =x_{3}\\ \dot{x}_{3} & =-10y^{\prime \prime }+5y^{3}-ty+t^{2}\\ & =-10x_{3}+5x_{1}^{3}-tx_{1}+t^{2} \end{align*}
The above \(\dot{x}\) vector is what returned back in the RHS call used by bvp4c. The following shows the solution obtained and the code used. One difficulty with this problem was to guess the correct initial solution to use. If the wrong guess was used, then Matlab gives an error
Four plots were generated, with different number of grid points, using \(N=20,40,60,80\) grid point, to see how the solution improves with more grid points added. At \(N=80\), the solution was smooth. Here are the results
problem description
solution
The first step is to convert the system to state space. (I used \(t\) below as the independent variable, instead of \(x\) as given in the problem statement, to reduce confusion with the \(x_{i}\) used for state space setup).\[ y^{\left ( 4\right ) }-y^{\prime }=e^{t}\] Let \(x_{1}=y,x_{2}=y^{\prime },x_{3}=y^{\prime \prime },x_{4}=y^{\prime \prime \prime }\), therefore\begin{align*} \dot{x}_{1} & =x_{2}\\ \dot{x}_{2} & =x_{3}\\ \dot{x}_{3} & =x_{4}\\ \dot{x}_{4} & =y^{\prime }+e^{t}\\ & =x_{1}+e^{t} \end{align*}
The above \(\dot{x}\) vector is what returned back in the RHS call used by bvp4c. Now we will solve it analytically in order to compare the solutions. The homogenous ODE is \(y^{\prime \prime \prime \prime }-y^{\prime }=0\), hence the characteristic equation is \(\lambda ^{4}-\lambda =0\), which has solutions\begin{align*} \lambda & =0\\ \lambda & =1^{\frac{1}{3}}\\ & =\left \{ 1,-\frac{1}{2}-\frac{\sqrt{3}}{2}i,-\frac{1}{2}+\frac{\sqrt{3}}{2}i\right \} \end{align*}
Therefore the homogenous solution is\begin{align*} y_{h} & =Ae^{\lambda _{1}}+Be^{\lambda _{2}}+Ce^{\lambda _{3}}+De^{\lambda _{4}}\\ & =A+Be^{t}+Ce^{\left ( -\frac{1}{2}-\frac{\sqrt{3}}{2}i\right ) t}+De^{\left ( -\frac{1}{2}+\frac{\sqrt{3}}{2}i\right ) t} \end{align*}
Since the homogenous solution contains \(e^{t}\) solution, and the forcing function is also \(e^{t}\), we can’t guess \(e^{t}\) as particular solution. We try \(y_{p}=cte^{t}\). Hence\begin{align*} y_{p}^{\prime } & =\left ( cte^{t}+ce^{t}\right ) \\ y_{p}^{\prime \prime } & =\left ( cte^{t}+ce^{t}\right ) +ce^{t}\\ y_{p}^{\prime \prime \prime } & =\left ( \left ( cte^{t}+ce^{t}\right ) +ce^{t}\right ) +ce^{t}\\ y_{p}^{\prime \prime \prime \prime } & =\left ( \left ( \left ( cte^{t}+ce^{t}\right ) +ce^{t}\right ) +ce^{t}\right ) +ce^{t} \end{align*}
Substituting this in the original ODE we obtain\begin{align*} \left ( \left ( \left ( cte^{t}+ce^{t}\right ) +ce^{t}\right ) +ce^{t}\right ) +ce^{t}-\left ( cte^{t}+ce^{t}\right ) & =e^{t}\\ 3ce^{t} & =e^{t} \end{align*}
Hence \(3c=1\) or \(c=\frac{1}{3}\), therefore \(y_{p}=\frac{1}{3}te^{t}\) and the full solution is\begin{align} y & =y_{h}+y_{p}\nonumber \\ & =A+Be^{t}+Ce^{\left ( -\frac{1}{2}-\frac{3}{2}i\right ) t}+De^{\left ( -\frac{1}{2}+\frac{3}{2}i\right ) t}+\frac{1}{3}te^{t}\nonumber \\ & =A+Be^{t}+Ce^{^{\frac{-1}{2}t}}e^{\frac{-3}{2}it}+De^{^{\frac{-1}{2}t}}e^{\frac{3}{2}it}+\frac{1}{3}te^{t}\nonumber \\ & =A+Be^{t}+e^{^{\frac{-1}{2}t}}\left ( C\cos \left ( \frac{\sqrt{3}}{2}t\right ) +D\sin \left ( \frac{\sqrt{3}}{2}t\right ) \right ) +\frac{1}{3}te^{t} \tag{A} \end{align}
Now we find the constants from initial and boundary conditions. At \(y\left ( 0\right ) =0\) hence\begin{equation} 0=A+B+C \tag{1} \end{equation} From \(y\left ( 1\right ) =1\,\) we obtain\begin{equation} 1=A+Be+e^{\frac{-1}{2}}\left ( C\cos \left ( \frac{\sqrt{3}}{2}\right ) +D\sin \left ( \frac{\sqrt{3}}{2}\right ) \right ) +\frac{1}{3}e \tag{2} \end{equation} To apply the other conditions, we need \(y^{\prime }\left ( t\right ) \). Taking derivative of (A) gives\[ y^{\prime }=Be^{t}+e^{^{\frac{-1}{2}t}}\left ( C\cos \left ( \frac{\sqrt{3}}{2}t\right ) +D\sin \left ( \frac{\sqrt{3}}{2}t\right ) \right ) +e^{^{\frac{-1}{2}t}}\left ( -C\frac{\sqrt{3}}{2}\sin \left ( \frac{\sqrt{3}}{2}t\right ) +D\frac{\sqrt{3}}{2}\cos \left ( \frac{\sqrt{3}}{2}t\right ) \right ) +\frac{1}{3}e^{t}+\frac{1}{3}te^{t}\] Using \(y^{\prime }\left ( 0\right ) =1\) gives\begin{equation} 0=B+C+D\frac{\sqrt{3}}{2}+\frac{1}{3} \tag{3} \end{equation} And from \(y^{\prime }\left ( 1\right ) =0\) we find\begin{equation} 0=Be+e^{^{\frac{-1}{2}}}\left ( C\cos \left ( \frac{\sqrt{3}}{2}\right ) +D\sin \left ( \frac{\sqrt{3}}{2}\right ) \right ) +e^{^{\frac{-1}{2}}}\left ( -C\frac{\sqrt{3}}{2}\sin \left ( \frac{\sqrt{3}}{2}\right ) +D\frac{\sqrt{3}}{2}\cos \left ( \frac{\sqrt{3}}{2}\right ) \right ) +\frac{2}{3}e \tag{4} \end{equation} Now we solve Eq (1,2,3,4) for \(A,B,C,D\). With the help of the computer, the above were now solved and the coefficients substituted back into (A) to give the analytical solution below\[ y\left ( t\right ) =\frac{1}{3}te^{t}-3.956e^{t}-16.1397e^{\frac{-1}{2}t}\cos \left ( \frac{\sqrt{3}}{2}t\right ) -6.2898e^{\frac{-1}{2}t}\sin \left ( \frac{\sqrt{3}}{2}t\right ) +20.096 \] Now that we have the analytical solution, we now need to find a solution using finite differences as well as per problem statement. The first step is to set up the grid. The following diagram shows the grid used
Total of \(N\) grid points is used. Since the solution is known at \(i=0\) and \(i=N-1\), the solution at the remaining only \(N-2\,\ \)points needs to be determined using finite difference scheme. We will now derive the FD equations for grid points \(i=1,2,3\). From grid point \(i=3\) to \(i=N\), the same pattern repeats, and the matrix \(A_{N\times N}\) will be filled using an iteration process as shown below.
The differential equation \[ y^{\prime \prime \prime \prime }\left ( x\right ) -y^{\prime }\left ( x\right ) =e^{t}\] In finite differences form is\[ \frac{y_{i-2}-4y_{i-1}+6y_{i}-4y_{i+1}+y_{i+2}}{h^{4}}-\frac{y_{i+1}-y_{i-1}}{2h}=e^{x_{i}}\] Where we used centered difference with \(O\left ( h^{2}\right ) \) local truncation error for the approximation of \(y^{\prime }\left ( x\right ) \) and \(5\) points centered difference for the approximation of \(y^{\prime \prime \prime \prime }\left ( x\right ) \)
At \(i=1\)\begin{align*} \frac{y_{-1}-4y_{0}+6y_{1}-4y_{2}+y_{3}}{h^{4}}-\frac{y_{2}-y_{0}}{2h} & =e^{x_{1}}\\ y_{-1}-4y_{0}+6y_{1}-4y_{2}+y_{3}-\frac{1}{2}h^{3}\left ( y_{2}-y_{0}\right ) & =h^{4}e^{x_{1}} \end{align*}
To find \(y_{-1}\), since \(y^{\prime }\left ( 0\right ) \) is known, using \(y^{\prime }\left ( 0\right ) =y_{0}^{\prime }=\frac{y_{1}-y_{-1}}{2h}\) given \(y_{-1}=y_{1}-2hy_{0}^{\prime }\) and the above becomes\begin{align*} \left ( y_{1}-2hy_{0}^{\prime }\right ) -4y_{0}+6y_{1}-4y_{2}+y_{3}-\frac{1}{2}h^{3}\left ( y_{2}-y_{0}\right ) & =h^{4}e^{x_{1}}\\ y_{0}\left ( -4+\frac{1}{2}h^{3}\right ) +7y_{1}+y_{2}\left ( -4-\frac{1}{2}h^{3}\right ) +y_{3} & =h^{4}e^{x_{1}}+2hy_{0}^{\prime }\\ 7y_{1}+y_{2}\left ( -4-\frac{1}{2}h^{3}\right ) +y_{3} & =h^{4}e^{x_{1}}+2hy_{0}^{\prime }+y_{0}\left ( 4-\frac{1}{2}h^{3}\right ) \end{align*}
At \(i=2\)\begin{align*} \frac{y_{0}-4y_{1}+6y_{2}-4y_{3}+y_{4}}{h^{4}}-\frac{y_{3}-y_{1}}{2h} & =e^{x_{i}}\\ y_{0}-4y_{1}+6y_{2}-4y_{3}+y_{4}-\frac{1}{2}h^{3}\left ( y_{3}-y_{1}\right ) & =h^{4}e^{x_{2}}\\ y_{1}\left ( -4+\frac{1}{2}h^{3}\right ) +6y_{2}+y_{3}\left ( -4-\frac{1}{2}h^{3}\right ) +y_{4} & =h^{4}e^{x_{2}}-y_{0} \end{align*}
At \(i=3\)\begin{align*} \frac{y_{1}-4y_{2}+6y_{3}-4y_{4}+y_{5}}{h^{4}}-\frac{y_{4}-y_{2}}{2h} & =e^{x_{3}}\\ y_{1}-4y_{2}+6y_{3}-4y_{4}+y_{5}-\frac{1}{2}h^{3}\left ( y_{4}-y_{2}\right ) & =h^{4}e^{x_{3}}\\ y_{1}+y_{2}\left ( -4+\frac{1}{2}h^{3}\right ) +6y_{3}+y_{4}\left ( -4-\frac{1}{2}h^{3}\right ) +y_{5} & =h^{4}e^{x_{3}} \end{align*}
The rest will now be repeated with \(i\) being increased by one for each new row, and shifted to the right by one for each row. For example for \(i=4\)\begin{align*} \frac{y_{2}-4y_{3}+6y_{4}-4y_{5}+y_{6}}{h^{4}}-\frac{y_{5}-y_{3}}{2h} & =e^{x_{4}}\\ y_{2}-4y_{3}+6y_{4}-4y_{5}+y_{6}-\frac{1}{2}h^{3}\left ( y_{5}-y_{3}\right ) & =h^{4}e^{x_{4}}\\ y_{2}+y_{3}\left ( -4+\frac{1}{2}h^{3}\right ) +6y_{4}+y_{5}\left ( -4+\frac{1}{2}h^{3}\right ) +y_{6} & =h^{4}e^{x_{4}} \end{align*}
Therefore, the \(Ax=b\) system to solve is
\[ \makebox [\textwidth ]{\fbox{ \scriptsize $ \begin{bmatrix} 7 & \left ( -4-\frac{1}{2}h^{3}\right ) & 1 & 0 & \cdots & 0 & 0 & 0\\ \left ( -4+\frac{1}{2}h^{3}\right ) & 6 & \left ( -4-\frac{1}{2}h^{3}\right ) & 1 & 0 & \cdots & 0 & 0\\ 1 & \left ( -4+\frac{1}{2}h^{3}\right ) & 6 & \left ( -4-\frac{1}{2}h^{3}\right ) & 1 & 0 & \cdots & 0\\ 0 & 1 & \left ( -4+\frac{1}{2}h^{3}\right ) & 6 & \left ( -4+\frac{1}{2}h^{3}\right ) & 1 & 0 & \cdots \\ 0 & 0 & \left ( -4+\frac{1}{2}h^{3}\right ) & 6 & \left ( -4+\frac{1}{2}h^{3}\right ) & 1 & 0 & \cdots \\ & & & & & & & \\ & & & & & & & \\ & & & & & & & \end{bmatrix}\begin{bmatrix} y_{1}\\ y_{2}\\ y_{3}\\ y_{4}\\ \vdots \\ y_{N-2}\\ y_{N-1}\\ y_{N}\end{bmatrix} =\begin{bmatrix} h^{4}e^{h}-2hy_{0}^{\prime }+y_{0}\left ( 4-\frac{1}{2}h^{3}\right ) \\ h^{4}e^{2h}-y_{0}\\ h^{4}e^{3h}\\ h^{4}e^{4h}\\ \vdots \\ \\ \\ \end{bmatrix} $}} \]
Now to fill in the last 3 rows. At \(i=N\)\[ \frac{y_{N-2}-4y_{N-1}+6y_{N}-4y_{N+1}+y_{N+2}}{h^{4}}-\frac{y_{N+1}-y_{N-1}}{2h}=e^{x_{i}}\] But we do not know \(y_{N+2}\) but since we know \(y^{\prime }\left ( 1\right ) =0\), then using \[ y^{\prime }\left ( 1\right ) =y_{N+1}^{\prime }=\frac{y_{N+2}-y_{N}}{2h}\] Hence \(y_{N+2}=2hy^{\prime }\left ( 1\right ) +y_{N}\) and the above becomes (by also replacing \(y_{N+1}=y\left ( 1\right ) \), which is known).\begin{align*} \frac{y_{N-2}-4y_{N-1}+6y_{N}-4y_{N+1}+2hy^{\prime }\left ( 1\right ) +y_{N}}{h^{4}}-\frac{y_{N+1}-y_{N-1}}{2h} & =e^{x_{N}}\\ y_{N-2}-4y_{N-1}+6y_{N}-4y\left ( 1\right ) +2hy^{\prime }\left ( 1\right ) +y_{N}-\frac{1}{2}h^{3}\left ( y\left ( 1\right ) -y_{N-1}\right ) & =h^{4}e^{x_{N}}\\ y_{N-2}+y_{N-1}\left ( -4+\frac{1}{2}h^{3}\right ) +7y_{N} & =h^{4}e^{x_{N}}-2hy^{\prime }\left ( 1\right ) +\left ( \frac{1}{2}h^{3}+4\right ) y\left ( 1\right ) \end{align*}
At \(i=N-1\)\begin{align*} \frac{y_{N-3}-4y_{N-2}+6y_{N-1}-4y_{N}+y_{N+1}}{h^{4}}-\frac{y_{N}-y_{N-2}}{2h} & =e^{x_{N-1}}\\ y_{N-3}-4y_{N-2}+6y_{N-1}-4y_{N}+y\left ( 1\right ) -\frac{1}{2}h^{3}\left ( y_{N}-y_{N-2}\right ) & =h^{4}e^{x_{N-1}}\\ y_{N-3}+y_{N-2}\left ( -4+\frac{1}{2}h^{3}\right ) +6y_{N-1}+y_{N}\left ( -4-\frac{1}{2}h^{3}\right ) & =h^{4}e^{x_{N-1}}-y\left ( 1\right ) \end{align*}
At \(i=N-2\)\begin{align*} \frac{y_{N-4}-4y_{N-3}+6y_{N-2}-4y_{N-1}+y_{N}}{h^{4}}-\frac{y_{N-1}-y_{N-3}}{2h} & =e^{x_{N-2}}\\ y_{N-4}-4y_{N-3}+6y_{N-2}-4y_{N-1}+y_{N}-\frac{1}{2}h^{3}\left ( y_{N-1}-y_{N-3}\right ) & =h^{4}e^{x_{N-2}}\\ y_{N-4}+y_{N-3}\left ( -4+\frac{1}{2}h^{3}\right ) +6y_{N-2}+y_{N-1}\left ( -4-\frac{1}{2}h^{3}\right ) +y_{N} & =h^{4}e^{x_{N-2}} \end{align*}
The \(Ax=b\) system becomes
The system is now solved for \(x\), which is the \(y_{i}\) solution and plotted. The Matlab code is given below.
The program nma_EMA_471_HW2_prob_2.m generates 4 result for different grid sizes. It uses \(8,15,30,100\) grid points each time, to compare the result. bvp4c and the finite difference method, both used the same grid size each time. Each time, the result is compare with the analytical solution (which used a much smaller grid than both).
At small number of grid points (large \(h\)), bvp4c seems to be more accurate at the boundary. Here is side by side showing the result for grid size of 8 points over the whole range.
The finite difference was also plotted against the bvp4c solution. The differences between them show up near where the solution changes most rapidly, around \(x=0.2\) and near the right boundary also. Here is the plot when using 8 grid points
To see more clearly the difference, the absolute difference between bvp4c and finite difference solution was plotted, by plotting \(\left | y_{\text{FDM}} - y_{\text{bvp4c}} \right | \) at each grid point
As more grid points added, the difference between bvp4c and the analytical solution became smaller. The same for finite difference solution. At 100 grid points, the largest absolute difference between bvp4c and FDM was \(0.00275\), near the middle of the range. This is compared to the difference being \(0.008\) when using 8 grid points.
This plot shows the difference at 100 grid points
A total of 16 plots generated (4 for each test case) as described above. Below all 16 plots are given, with description title of each plot. Then the source code used to generate these plots is listed.
In conclusion: As more grid points added, both bvp4c and FDM approached the analytical solution. There remains difference between bvp4c and FDM in the middle range. Both converged to the same result at the boundaries. This is expected, since the solution there is given from the problem boundary conditions.
problem description
solution
The first step is to determine the analytical solution, in order to compare with the numerical solutions. The ODE to solve is (below, \(\sigma \) is used instead of \(\sigma _{r}\) to simplify the notation)\begin{equation} r\frac{d^{2}\sigma }{dr^{2}}+3\frac{d\sigma }{dr}=-\left ( 3+\upsilon \right ) \rho \omega ^{2}r \tag{1} \end{equation} Since \(\upsilon ,\rho ,\omega \) are all constants, the above can be written as\[ r\frac{d^{2}\sigma }{dr^{2}}+3\frac{d\sigma }{dr}=kr \] Where \(k=-\left ( 3+\upsilon \right ) \rho \omega ^{2}\). To solve the above, we introduce \(f=\frac{d\sigma }{dr}\) and it becomes a first order ODE\begin{align} r\frac{df}{dr}+3f & =kr\nonumber \\ \frac{df}{dr}+\frac{3}{r}f & =k \tag{2} \end{align}
This is separable now, and can be solved for \(f\). Looking at the homogenous ODE first,\begin{align*} \frac{df}{dr}+\frac{3}{r}f & =0\\ \frac{df}{f} & =-\frac{3}{r}dr \end{align*}
Integrating both sides\begin{align*} \ln f & =-3\ln r+c_{1}\\ f & =e^{-3\ln r+c_{1}}\\ & =c_{2}e^{-3\ln r} \end{align*}
Hence the homogenous solution is\[ f_{h}=\frac{c_{2}}{r^{3}}\] Where \(c_{2}\) is constant of integration. To find the particular solution \(f_{p}\), and since the RHS of (2) is constant \(k\), then we guess \(f_{p}=crk\), where \(c\) is constant. Hence (2) becomes\begin{align*} ck+\frac{3}{r}rck & =k\\ 4c & =1\\ c & =\frac{1}{4} \end{align*}
Hence \[ f_{p}=\frac{r}{4}k \] And the complete solution is\begin{align*} f & =f_{h}+f_{p}\\ & =\frac{c_{2}}{r^{3}}+\frac{r}{4}k \end{align*}
Now that we found \(f\left ( r\right ) \), and since \(f=\frac{d\sigma }{dr}\) then we have\[ \frac{d\sigma }{dr}=\frac{c_{2}}{r^{3}}+\frac{r}{4}k \] Which is separable. Hence\begin{align*} d\sigma & =\left ( \frac{c_{2}}{r^{3}}+\frac{r}{4}k\right ) dr\\ \sigma & ={\displaystyle \int } \left ( \frac{c_{2}}{r^{3}}+\frac{r}{4}k\right ) dr\\ & =\frac{-c_{2}}{2}\frac{1}{r^{2}}+\frac{r^{2}}{8}k+c_{3} \end{align*}
Hence the analytical solution is\begin{equation} \sigma _{r}\left ( r\right ) =\frac{-C_{1}}{2}\frac{1}{r^{2}}+\frac{r^{2}}{8}k+C_{2} \tag{3} \end{equation} Where \(C_{1},C_{2}\) are constants of integration, which we will now find from boundary conditions. At \(r=r_{i},\sigma _{r}=0\) and at \(r=r_{o},\sigma _{r}=0\), hence we obtain two equations\begin{align*} 0 & =\frac{-C_{1}}{2}\frac{1}{r_{i}^{2}}+\frac{r_{i}^{2}}{8}k+C_{2}\\ 0 & =\frac{-C_{1}}{2}\frac{1}{r_{o}^{2}}+\frac{r_{o}^{2}}{8}k+C_{2} \end{align*}
Solving these gives\begin{align*} C_{1} & =\frac{-k}{4}r_{i}^{2}r_{o}^{2}\\ C_{2} & =\frac{-k}{8}\left ( r_{i}^{2}+r_{o}^{2}\right ) \end{align*}
Therefore (3) becomes\[ \sigma _{r}\left ( r\right ) =\frac{k}{8}\left ( r_{i}^{2}r_{o}^{2}\frac{1}{r^{2}}+r^{2}-\left ( r_{i}^{2}+r_{o}^{2}\right ) \right ) \] But \(k=-\left ( 3+\upsilon \right ) \rho \omega ^{2}\), hence the above becomes\begin{equation} \fbox{$\sigma _r\left ( r\right ) =\frac{\left ( 3+\upsilon \right ) \rho \omega ^2}{8}\left ( \left ( r_i^2+r_o^2\right ) -r_i^2r_o^2\frac{1}{r^2}-r^2\right ) $} \tag{4} \end{equation} Now we will find the analytical solution for the hoop stress\[ \sigma _{\theta }=r\frac{d\sigma }{dr}+\sigma _{r}+\rho \omega ^{2}r^{2}\] From (4), we see that\[ \frac{d\sigma }{dr}=\frac{\left ( 3+\upsilon \right ) \rho \omega ^{2}}{8}\left ( \frac{2}{r^{3}}r_{i}^{2}r_{o}^{2}-2r\right ) \] Hence\begin{equation} \fbox{$\sigma _\theta =\frac{\left ( 3+\upsilon \right ) \rho \omega ^2}{8}\left ( \frac{2}{r^2}r_i^2r_o^2-2r^2\right ) +\sigma _r+\rho \omega ^2r^2$} \tag{5} \end{equation} Now that we have found the analytical solutions, we can implement the numerical solution and compare. We start with bvp4c to verify the analytical solution with, then implement the finite difference method. We need first to convert the ODE to state space.\[ \frac{d^{2}\sigma }{dr^{2}}+\frac{3}{r}\frac{d\sigma }{dr}=-\left ( 3+\upsilon \right ) \rho \omega ^{2}\] Let \(x_{1}=\sigma _{r}\), and let \(x_{2}=\frac{d\sigma }{dr}\), hence\begin{align*} \dot{x}_{1} & =x_{2}\\ \dot{x}_{2} & =-\left ( 3+\upsilon \right ) \rho \omega ^{2}-\frac{3}{r}x_{2} \end{align*}
The boundary conditions are (using bvp4c notations) \(y_{a}\left ( 1\right ) =0,y_{b}\left ( 1\right ) =0\,\).
\[ \frac{d^{2}\sigma }{dr^{2}}+3\frac{1}{r}\frac{d\sigma }{dr}=-\left ( 3+\upsilon \right ) \rho \omega ^{2}\] The following grid is used
The grid is only over \(r=0.01\cdots 0.1\) meters since this is the distance between the internal radius and the outer radius. For \(\frac{d^{2}\sigma }{dr^{2}}\) we use 3 points centered difference approximation, for \(i=1\cdots N\) which has \(O\left ( h^{2}\right ) \) local truncation error\[ \frac{d^{2}\sigma }{dr^{2}}=\frac{\sigma _{i+1}-2\sigma _{i}+\sigma _{i-1}}{h^{2}}\] Where \(\sigma _{0},\sigma _{N+1}\) are given as the boundary conditions (both are zero). For \(\frac{d\sigma }{dr}\) we use two points centered difference, which also has \(O\left ( h^{2}\right ) \) local truncation error\[ \frac{d\sigma }{dr}=\frac{\sigma _{i+1}-\sigma _{i-1}}{2h}\] Therefore, the ODE becomes \begin{align*} r_{i}\frac{\sigma _{i+1}-2\sigma _{i}+\sigma _{i-1}}{h^{2}}+3\frac{\sigma _{i+1}-\sigma _{i-1}}{2h} & =-\left ( 3+\upsilon \right ) \rho \omega ^{2}r_{i}\\ \sigma _{i+1}-2\sigma _{i}+\sigma _{i-1}+\frac{3}{2r_{i}}h\left ( \sigma _{i+1}-\sigma _{i-1}\right ) & =-h^{2}\left ( 3+\upsilon \right ) \rho \omega ^{2}\\ \sigma _{i-1}\left ( 1-\frac{3}{2r_{i}}h\right ) -2\sigma _{i}+\sigma _{i+1}\left ( 1+\frac{3}{2r_{i}}h\right ) & =-h^{2}\left ( 3+\upsilon \right ) \rho \omega ^{2} \end{align*}
Where \(r_{i}=r_{i}+ih=0.01+ih\), since are starting from \(0.01\), the above becomes\[ \sigma _{i-1}\left ( 1-\frac{3}{2\left ( r_{i}+ih\right ) }h\right ) -2\sigma _{i}+\sigma _{i+1}\left ( 1+\frac{3}{2\left ( r_{i}+ih\right ) }h\right ) =-h^{2}\left ( 3+\upsilon \right ) \rho \omega ^{2}\] For \(i=1\), \[ \sigma _{0}\left ( 1-\frac{3}{2\left ( r_{i}+h\right ) }h\right ) -2\sigma _{1}+\sigma _{2}\left ( 1+\frac{3}{2\left ( r_{i}+h\right ) }h\right ) =-h^{2}\left ( 3+\upsilon \right ) \rho \omega ^{2}\] But \(\sigma _{0}\) is the boundary conditions, which we move to the RHS, hence\[ -2\sigma _{1}+\sigma _{2}\left ( 1+\frac{3}{2\left ( r_{i}+h\right ) }h\right ) =-h^{2}\left ( 3+\upsilon \right ) \rho \omega ^{2}-\sigma _{0}\left ( 1-\frac{3}{2\left ( r_{i}+h\right ) }h\right ) \] For \(i=2\)\[ \sigma _{1}\left ( 1-\frac{3}{2\left ( r_{i}+2h\right ) }h\right ) -2\sigma _{2}+\sigma _{3}\left ( 1+\frac{3}{2\left ( r_{i}+2h\right ) }h\right ) =-h^{2}\left ( 3+\upsilon \right ) \rho \omega ^{2}\] For \(i=3\)\[ \sigma _{2}\left ( 1-\frac{3}{2\left ( r_{i}+3h\right ) }h\right ) -2\sigma _{3}+\sigma _{4}\left ( 1+\frac{3}{2\left ( r_{i}+3h\right ) }h\right ) =-h^{2}\left ( 3+\upsilon \right ) \rho \omega ^{2}\] The same pattern repeats. For \(i=N\)\[ \sigma _{N-1}\left ( 1-\frac{3}{2\left ( r_{i}+Nh\right ) }h\right ) -2\sigma _{N}+\sigma _{N+1}\left ( 1+\frac{3}{2\left ( r_{i}+Nh\right ) }h\right ) =-h^{2}\left ( 3+\upsilon \right ) \rho \omega ^{2}\] But \(\sigma _{N+1}\) is given (the right side boundary conditions, hence\[ \sigma _{N-1}\left ( 1-\frac{3}{2\left ( r_{i}+Nh\right ) }h\right ) -2\sigma _{N}=-h^{2}\left ( 3+\upsilon \right ) \rho \omega ^{2}-\sigma _{N+1}\left ( 1+\frac{3}{2\left ( r_{i}+Nh\right ) }h\right ) \] And for \(i=N-1\)\[ \sigma _{N-2}\left ( 1-\frac{3}{2\left ( r_{i}+\left ( N-1\right ) h\right ) }h\right ) -2\sigma _{N-1}+\sigma _{N}\left ( 1+\frac{3}{2\left ( r_{i}+\left ( N-1\right ) h\right ) }h\right ) =-h^{2}\left ( 3+\upsilon \right ) \rho \omega ^{2}\] Therefore the \(Ax=b\) system is below. The \(A\) matrix is\[\begin{bmatrix} -2 & 1+\frac{3h}{2\left ( r_{i}+h\right ) } & 0 & \cdots & \cdots & \cdots & \cdots & \cdots \\ 1-\frac{3h}{2\left ( r_{i}+2h\right ) } & -2 & 1+\frac{3h}{2\left ( r_{i}+2h\right ) } & 0 & \cdots & \cdots & \cdots & \cdots \\ 0 & 1-\frac{3h}{2\left ( r_{i}+3h\right ) } & -2 & 1+\frac{3h}{2\left ( r_{i}+3h\right ) } & 0 & \cdots & \cdots & \cdots \\ \vdots & 0 & 0 & 0 & 0 & \ddots & \cdots & \cdots \\ \vdots & \vdots & \vdots & \vdots & \vdots & \cdots & \cdots & \cdots \\ \vdots & \vdots & \vdots & \vdots & \vdots & \cdots & \cdots & \cdots \\ \vdots & \vdots & \vdots & \vdots & 0 & 1-\frac{3h}{2\left ( r_{i}+\left ( N-1\right ) h\right ) } & -2 & 1+\frac{3h}{2\left ( r_{i}+\left ( N-1\right ) h\right ) }\\ 0 & 0 & 0 & 0 & 0 & 0 & 1-\frac{3h}{2\left ( r_{i}+Nh\right ) } & -2 \end{bmatrix} \] And the \(x\) vector is \(\begin{bmatrix} \sigma _{1} & \sigma _{1} & \cdots & \sigma _{N-2} & \sigma _{N-1} & \sigma _{N}\end{bmatrix} ^{T}\) and the \(b\) vector is\[\begin{bmatrix} -h^{2}\left ( 3+\upsilon \right ) \rho \omega ^{2}-\sigma _{0}\left ( 1-\frac{3h}{2\left ( r_{i}+h\right ) }\right ) \\ -h^{2}\left ( 3+\upsilon \right ) \rho \omega ^{2}\\ -h^{2}\left ( 3+\upsilon \right ) \rho \omega ^{2}\\ \vdots \\ \vdots \\ -h^{2}\left ( 3+\upsilon \right ) \rho \omega ^{2}\\ -h^{2}\left ( 3+\upsilon \right ) \rho \omega ^{2}\\ -h^{2}\left ( 3+\upsilon \right ) \rho \omega ^{2}-\sigma _{N+1}\left ( 1+\frac{3h}{2\left ( r_{i}+Nh\right ) }\right ) \end{bmatrix} \]
Finding numerical solution for hoop stress To compare the analytical solution for \(\sigma _{\theta }\) with the numerical one, we need to evaluate \(\sigma _{\theta }=r\frac{d\sigma }{dr}+\sigma _{r}+\rho \omega ^{2}r^{2}\) numerically from the numerical solution we found about using finite difference method. Using finite difference, this expression becomes\[ \sigma _{\theta _{i}}=\left ( r_{int}+ih\right ) \frac{\sigma _{i+1}-\sigma _{i-1}}{2h}+\sigma _{i}+\rho \omega ^{2}\left ( r_{int}+ih\right ) ^{2}\] Where \(r_{int}=0.01\) meter. Since we do not know \(\sigma _{-1}\) and can not obtain derivative at the edge to approximate using phantom grid point, we will start from \(i=1\cdots N\), so that \(\sigma _{i-1}=\sigma _{0}\) which is known.
Hence for \(i=1\) we have\[ \sigma _{\theta _{1}}=\left ( r_{int}+h\right ) \frac{\sigma _{2}-\sigma _{0}}{2h}+\sigma _{1}+\rho \omega ^{2}\left ( r_{int}+h\right ) ^{2}\] And for \(i=N\)\begin{equation} \sigma _{\theta _{N}}=\left ( r_{int}+Nh\right ) \frac{\sigma _{N+1}-\sigma _{N-1}}{2h}+\sigma _{N}+\rho \omega ^{2}\left ( r_{int}+Nh\right ) ^{2} \tag{6} \end{equation} In the above, \(\sigma _{N+1}\) and \(\sigma _{0}\) are the boundary conditions we are given. All other \(\sigma _{i}\) values are obtained from the numerical solution we did in the last section (the finite difference solution). This was implemented in Matlab.
The program was run for \(10,15,20,25,30,100\) grid points. Each time, 4 plots were generated:
The result shows that bvp4c was more accurate than FDM for small number of grid points (large \(h\)) since the bvp4c curve was closer to the analytical solution than FDM curve. To get more accurate numerical hoop stress, the number of grid points needed to be over \(50\). At \(N=100\) grid points, the result of \(\sigma _{\theta }\) matched the analytical solution extremely well as can be seen from the plots below. The following table gives the maximum radial stress found and the location as reported by bvp4c and FDM based methods. Units for stress is \(N/m^{2}\) and units for \(r\) is meters. This shows the maximum radial stress occurs at left edge \(r=0.01\) meter as expected. This is where the inner hole edge starts. The maximum radial stress is located at \(r=0.032\) meter. About \(2\) cm after the inner hole edge.
| \(N\) (grid points) | Max \(\sigma _{r}\) FDM based, and \(r\) location | Max \(\sigma _{r}\) bvp4c based, and location |
| \(10\) | \(\sigma _{r}=\)\(1,783,600\) at \(r=0.03\) | \(\sigma _{r}=\)\(1,733,752\) at \(r=0.03\) |
| \(15\) | \(\sigma _{r}=\)\(1,731,188\) at \(r=0.29\) | \(\sigma _{r}=\)\(1,752,483\) at \(r=0.029\) |
| \(20\) | \(\sigma _{r}=\)\(1,732,767\) at \(r=0.034\) | \(\sigma _{r}=\)\(1,741,523\) at \(r=0.034\) |
| \(25\) | \(\sigma _{r}=\)\(1,735,518\) at \(r=0.033\) | \(\sigma _{r}=\)\(1,741,592\) at \(r=0.033\) |
| \(30\) | \(\sigma _{r}=\)\(1,735,984\) at \(r=0.032\) | \(\sigma _{r}=\)\(1,740,593\) at \(r=0.032\) |
| \(100\) | \(\sigma _{r}=\)\(1,736,401\) at \(r=0.032\) | \(\sigma _{r}=\)\(1,736,759\) at \(r=0.032\) |
The following table gives the maximum hoop stress \(\sigma _{\theta }\) found and the location as reported FDM based method. Units for stress is \(N/m^{2}\) and units for \(r\) is meters. Each plot below also display this information in the title.
| \(N\) (grid points) | Max \(\sigma _{_{\theta }}\) value and location |
| \(10\) | \(\sigma _{_{\theta }}=\)\(3,626,000\) at \(r=0.01\) |
| \(15\) | \(\sigma _{_{\theta }}=\)\(3,639,120\) at \(r=0.01\) |
| \(20\) | \(\sigma _{_{\theta }}=\)\(3,665,659\) at \(r=0.01\) |
| \(25\) | \(\sigma _{_{\theta }}=\)\(3,697,797\) at \(r=0.01\) |
| \(30\) | \(\sigma _{_{\theta }}=\)\(3,730,811\) at \(r=0.01\) |
| \(100\) | \(\sigma _{_{\theta }}=\)\(4,004,798\) at \(r=0.01\) |
There are total of 24 plots below. Four plots for each \(N\). This helps show that more grid points resulted in more accurate numerical solution.
