### 2.2 HW 2

2.2.1 Problem 1
2.2.2 Problem 2
2.2.3 Problem 3

1.
This HW in one PDF (letter size)
2.
This HW in one PDF (legal size)

#### 2.2.1 Problem 1

2.2.1.1 Output
2.2.1.2 Source code

problem description

solution

The ﬁrst 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 diﬃculty with this problem was to guess the correct initial solution to use. If the wrong guess was used, then Matlab gives an error

Unable to solve the collocation equations -- a singular Jacobian encountered.

##### 2.2.1.1 Output

Four plots were generated, with diﬀerent 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

#### 2.2.2 Problem 2

2.2.2.1 Results
2.2.2.2 Source code

problem description

solution

The ﬁrst 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 ﬁnd the constants from initial and boundary conditions. At $$y\left ( 0\right ) =0$$ hence$$0=A+B+C \tag{1}$$ From $$y\left ( 1\right ) =1\,$$ we obtain$$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}$$ 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$$0=B+C+D\frac{\sqrt{3}}{2}+\frac{1}{3} \tag{3}$$ And from $$y^{\prime }\left ( 1\right ) =0$$ we ﬁnd$$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}$$ 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 coeﬃcients 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 ﬁnd a solution using ﬁnite diﬀerences as well as per problem statement. The ﬁrst 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 ﬁnite diﬀerence 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 ﬁlled using an iteration process as shown below.

The diﬀerential equation $y^{\prime \prime \prime \prime }\left ( x\right ) -y^{\prime }\left ( x\right ) =e^{t}$ In ﬁnite diﬀerences 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 diﬀerence with $$O\left ( h^{2}\right )$$ local truncation error for the approximation of $$y^{\prime }\left ( x\right )$$ and $$5$$ points centered diﬀerence 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 ﬁnd $$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 ﬁll 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

$\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 & 1 & \left ( -4+\frac{1}{2}h^{3}\right ) & 6 & \left ( -4+\frac{1}{2}h^{3}\right ) & 1 & 0\\ \cdots & \cdots & \ddots & \cdots & \ddots & \cdots & \cdots & \cdots \\ 0 & \cdots & \left ( -4+\frac{1}{2}h^{3}\right ) & 6 & \left ( -4+\frac{1}{2}h^{3}\right ) & 1 & 0 & \cdots \\ & & & 1 & \left ( -4+\frac{1}{2}h^{3}\right ) & 6 & \left ( -4-\frac{1}{2}h^{3}\right ) & 1\\ 0 & 0 & \cdots & 0 & 1 & \left ( -4+\frac{1}{2}h^{3}\right ) & 6 & \left ( -4-\frac{1}{2}h^{3}\right ) \\ 0 & 0 & 0 & \cdots & 0 & 1 & \left ( -4+\frac{1}{2}h^{3}\right ) & 7 \end{bmatrix}\begin{bmatrix} y_{1}\\ y_{2}\\ y_{3}\\ y_{4}\\ \vdots \\ \vdots \\ \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\left ( 0\right ) \\ h^{4}e^{3h}\\ h^{4}e^{4h}\\ \vdots \\ \vdots \\ \vdots \\ h^{4}e^{x_{N-2}}\\ h^{4}e^{x_{N-1}}-y\left ( 1\right ) \\ h^{4}e^{x_{N}}-2hy^{\prime }\left ( 1\right ) +\left ( \frac{1}{2}h^{3}+4\right ) y\left ( 1\right ) \end{bmatrix}$

The system is now solved for $$x$$, which is the $$y_{i}$$ solution and plotted. The Matlab code is given below.

##### 2.2.2.1 Results

The program nma_EMA_471_HW2_prob_2.m generates 4 result for diﬀerent grid sizes. It uses $$8,15,30,100$$ grid points each time, to compare the result. bvp4c and the ﬁnite diﬀerence 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 ﬁnite diﬀerence was also plotted against the bvp4c solution. The diﬀerences 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 diﬀerence, the absolute diﬀerence between bvp4c and ﬁnite diﬀerence solution was plotted, by plotting $$\left | y_{\text{FDM}} - y_{\text{bvp4c}} \right |$$ at each grid point

As more grid points added, the diﬀerence between bvp4c and the analytical solution became smaller. The same for ﬁnite diﬀerence solution. At 100 grid points, the largest absolute diﬀerence between bvp4c and FDM was $$0.00275$$, near the middle of the range. This is compared to the diﬀerence being $$0.008$$ when using 8 grid points.

This plot shows the diﬀerence 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 diﬀerence 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.

8 grid points plots

15 grid points plots

30 grid points plots

100 grid points plots

#### 2.2.3 Problem 3

2.2.3.1 Finite diﬀerence method
2.2.3.2 Results
2.2.3.3 source code

problem description

solution

The ﬁrst 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)$$r\frac{d^{2}\sigma }{dr^{2}}+3\frac{d\sigma }{dr}=-\left ( 3+\upsilon \right ) \rho \omega ^{2}r \tag{1}$$ 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 ﬁrst 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 ﬁrst,\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 ﬁnd 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$$\sigma _{r}\left ( r\right ) =\frac{-C_{1}}{2}\frac{1}{r^{2}}+\frac{r^{2}}{8}k+C_{2} \tag{3}$$ Where $$C_{1},C_{2}$$ are constants of integration, which we will now ﬁnd 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$$\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}$$ Now we will ﬁnd 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$$\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}$$ 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 ﬁnite diﬀerence method. We need ﬁrst 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\,$$.

##### 2.2.3.1 Finite diﬀerence method

$\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 diﬀerence 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 diﬀerence, 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 ﬁnite diﬀerence method. Using ﬁnite diﬀerence, 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$$$$\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}$$ 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 ﬁnite diﬀerence solution). This was implemented in Matlab.

##### 2.2.3.2 Results

The program was run for $$10,15,20,25,30,100$$ grid points. Each time, 4 plots were generated:

1.
Plot comparing the analytical and FDM result for $$\sigma _{r}$$ Title contains location and value of maximum $$\sigma _{r}$$ reported by FDM based method.
2.
Plot comparing the analytical and bvp4c result for $$\sigma _{r}$$ Title contains location and value of maximum $$\sigma _{r}$$ reported by bvp4c based method
3.
Plot comparing bvp4c and FDM result for $$\sigma _{r}$$
4.
Plot comparing analytical result and FDM for $$\sigma _{\theta }$$ Title contains location and value of maximum $$\sigma _{\theta }$$ reported by FDM based method.

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.

10 grid points

15 grid points

20 grid points

25 grid points

30 grid points

100 grid points