### 2.5 HW 5

2.5.1 Problem 1
2.5.2 Problem 2
2.5.3 Problem 3

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

#### 2.5.1 Problem 1 A polynomial $$f(x)$$ of order $$p$$ is integrated exactly with the Gaussian quadrature method using $$\frac{p+1}{2}$$ number of Gaussian points.

Hence $$f_1(x)=x^5$$ needs $$\frac{5+1}{2}=3$$ Gaussian points and $$f_2(x)=x^7$$ needs $$\frac{7+1}{2}=4$$ Gaussian points for exact result.

The integral $$\int _a^b f(x) \, dx$$ is ﬁrst converted to be in the domain $$\{-1,+1\}$$ as follows $\int _a^b f(x) \, dx = \frac{b-a}{2} \int _{-1}^{+1} f\left ( \frac{b-a}{2}t+\frac{b+a}{2} \right )\, dt$ For a polynomial $$f(x)$$ of order $$p=3$$, two Gaussian points are needed to evaluate the above integral exactly. Therefore the above integral simpliﬁes to $\int _a^b f(x) \, dx = A \left ( w_1 f(A t_1+B) + w_2 f(A t_2+B) \right )$ Where \begin{align*} A & =\frac{b-a}{2}\\ B & =\frac{b+a}{2} \end{align*}

And $$w_i$$ is the weight at location $$t_i$$. The weights and location of the weights are obtained from tables. For higher order polynomials, more points and weights are needed.

In general, using $$N$$ points the integral is \begin{align*} \int _a^b f(x) \, dx &= \frac{b-a}{2} \sum _{i=1}^N w_i f\left ( \frac{b-a}{2} t_i + \frac{a+b}{2} \right ) \\ &= A \sum _{i=1}^N w_i f(A t_i + B) \tag{1} \end{align*}

The program nma_EMA_471_HW5_problem_1.m integrates the above two polynomials $$f_1(x),f_2(x)$$ using Gaussian quadrature method (1) using $$N=3$$ and $$N=4$$ points respectively and compares the result of each to the analytical solution. The following table shows the result.

 function analytical result $$N$$ Gaussian quadrature result $$f_1(x)=x^4$$ $$\int _0^4 x^4\,dx = \frac{2048}{3} = 682.666666666667$$ 3 $$682.666666666667$$ $$f_1(x)=x^7$$ $$\int _0^4 x^7\,dx = 8192$$ 4 $$8192$$

The result is exact. Note: The above Matlab program used the exact weights and points for Gaussian quadrature as given in https://en.wikipedia.org/wiki/Gaussian_quadrature

#### 2.5.2 Problem 2

2.5.2.1 part a
2.5.2.2 part b  ##### 2.5.2.1 part a

The program nma_EMA_471_HW5_problem_2_part_a.m implements the ﬁrst part of this problem.

The following table shows the result of the computation. It shows the result of the integral using Gaussian quadrature for diﬀerent number of points with the relative error against Matlab’s Quad (integral) command.

 Number of points relative error (percentage) value of integral 2 $$79.845442$$ $$0.023377470178383$$ 3 $$23.892753$$ $$0.143704430262801$$ 4 $$3.060319$$ $$0.112441294428254$$ 5 $$4.692010$$ $$0.121433298541329$$ 6 $$0.019015$$ $$0.116013045399658$$ 7 $$0.011820$$ $$0.116004700249995$$ 8 $$0.001535$$ $$0.115989208171974$$

The following is a plot of the above data  ##### 2.5.2.2 part b

The program nma_EMA_471_HW5_problem_2_part_b.m implements the second part of this problem. By breaking the domain into 2 parts, the following table shows the result of the computation. It shows the result of the integral using Gaussian quadrature for $$5$$ and $$6$$ points with the relative error against Matlab’s Quad (integral) command. The integration was done on each subdomain and the results added.

 Number of points relative error (percentage) value of integral using Gaussian quadrature 5 $$1.231392$$ $$0.114562685084637$$ 6 $$0.000303$$ $$0.115990636860983$$

The above shows clearly that by breaking the domain into two smaller parts, and adding each result, the ﬁnal result of Gaussian quadrature improved compared to part(a) where one large domain was used. This makes sense. Because we have eﬀectively used more sampling points in part(b) compared to part(a) when looking at the whole domain.

This shows that, to obtain more accuracy using Gaussian quadrature, and still use the same number of points $$N$$, then we can break the domain into smaller regions, and use $$N$$ on each region, and add the result obtained from each region.

To see the diﬀerence between part(a) and (b) more clearly, the following table shows the result for 5 and 6 points side by side from part(a) and part(b). The table below shows the relative error is much smaller for part(b).

 Number of points relative error part(b) relative error part(a) 5 $$1.231392$$ $$4.692010$$ 6 $$0.000303$$ $$0.019015$$  #### 2.5.3 Problem 3

2.5.3.1 shape functions
2.5.3.2 $$x(s,t,r)$$ terms
2.5.3.3 $$y(r,s,t)$$ terms
2.5.3.4 $$z(r,s,t)$$ terms
2.5.3.5 results   ##### 2.5.3.1 shape functions

The following are the shape functions \begin{align*} f_I&=\frac{1}{8} (1-r) (1-s) (1-t) (-r-s-t-2)\\ f_J&=\frac{1}{8} (1-r) (s+1) (1-t) (-r+s-t-2)\\ f_K&=\frac{1}{8} (1-r) (s+1) (t+1) (-r+s+t-2)\\ f_L&=\frac{1}{8} (1-r) (1-s) (t+1) (-r-s+t-2)\\ f_M&=\frac{1}{8} (r+1) (1-s) (1-t) (r-s-t-2)\\ f_N&=\frac{1}{8} (r+1) (s+1) (1-t) (r+s-t-2)\\ f_O&=\frac{1}{8} (r+1) (s+1) (t+1) (r+s+t-2)\\ f_P&=\frac{1}{8} (r+1) (1-s) (t+1) (r-s+t-2)\\ f_Q&=\frac{1}{4} (1-r) \left (1-s^2\right ) (1-t)\\ f_R&=\frac{1}{4} (1-r) (s+1) \left (1-t^2\right )\\ f_S&=\frac{1}{4} (1-r) \left (1-s^2\right ) (t+1)\\ f_T&=\frac{1}{4} (1-r) (1-s) \left (1-t^2\right )\\ f_U&=\frac{1}{4} (r+1) \left (1-s^2\right ) (1-t)\\ f_V&=\frac{1}{4} (r+1) (s+1) \left (1-t^2\right )\\ f_W&=\frac{1}{4} (r+1) \left (1-s^2\right ) (t+1)\\ f_X&=\frac{1}{4} (r+1) (1-s) \left (1-t^2\right )\\ f_Y&=\frac{1}{4} \left (1-r^2\right ) (1-s) (1-t)\\ f_Z&=\frac{1}{4} \left (1-r^2\right ) (s+1) (1-t)\\ f_A&=\frac{1}{4} \left (1-r^2\right ) (s+1) (t+1)\\ f_B&=\frac{1}{4} \left (1-r^2\right ) (1-s) (t+1)\\ \end{align*}

To obtain the Jacobian, we need to obtain the determinant of $\begin{pmatrix} \frac{\partial{x}}{\partial{r}}&\frac{\partial{y}}{\partial{r}}&\frac{\partial{z}}{\partial{r}}\\ \frac{\partial{x}}{\partial{s}}&\frac{\partial{y}}{\partial{s}}&\frac{\partial{z}}{\partial{s}}\\ \frac{\partial{x}}{\partial{t}}&\frac{\partial{y}}{\partial{t}}&\frac{\partial{z}}{\partial{t}}\\ \end{pmatrix}$ Where \begin{align*} x(r,s,t) =& \sum _{i=I}^{I=B} x_i f_i(r,s,t)\\ y(r,s,t) =& \sum _{i=I}^{I=B} y_i f_i(r,s,t)\\ z(r,s,t) =& \sum _{i=I}^{I=B} z_i f_i(r,s,t) \end{align*}

Once we determine $$x(r,s,t),y(r,s,t),z(r,s,t)$$ from the above, then we can determine the Jacobian determinant at each Gaussian point.

##### 2.5.3.2 $$x(s,t,r)$$ terms

From the above sum, expanding $$x(r,s,t)$$ gives \begin{align*} x_ =&x_I\,f_I+x_J\,f_J+x_K\,f_K+x_L\,f_L+x_M\,f_M+x_N\,f_N+x_O\,f_O+x_P\,f_P+x_Q\,f_Q+x_R\,f_R+x_S\,f_S+x_T\,f_T+x_U\,f_U+x_V\,f_V+x_W\,f_W+x_X\,f_X+x_Y\,f_Y+x_Z\,f_Z+x_A\,f_A+x_B\,f_B \end{align*}

Substituting the values of $$f_i$$ into the above results in \begin{align*} x(r,s,t)=&x_I\,\frac{1}{8} (1-r) (1-s) (1-t) (-r-s-t-2)+\\ &x_J\,\frac{1}{8} (1-r) (s+1) (1-t) (-r+s-t-2)+\\ &x_K\,\frac{1}{8} (1-r) (s+1) (t+1) (-r+s+t-2)+\\ &x_L\,\frac{1}{8} (1-r) (1-s) (t+1) (-r-s+t-2)+\\ &x_M\,\frac{1}{8} (r+1) (1-s) (1-t) (r-s-t-2)+\\ &x_N\,\frac{1}{8} (r+1) (s+1) (1-t) (r+s-t-2)+\\ &x_O\,\frac{1}{8} (r+1) (s+1) (t+1) (r+s+t-2)+\\ &x_P\,\frac{1}{8} (r+1) (1-s) (t+1) (r-s+t-2)+\\ &x_Q\,\frac{1}{4} (1-r) \left (1-s^2\right ) (1-t)+\\ &x_R\,\frac{1}{4} (1-r) (s+1) \left (1-t^2\right )+\\ &x_S\,\frac{1}{4} (1-r) \left (1-s^2\right ) (t+1)+\\ &x_T\,\frac{1}{4} (1-r) (1-s) \left (1-t^2\right )+\\ &x_U\,\frac{1}{4} (r+1) \left (1-s^2\right ) (1-t)+\\ &x_V\,\frac{1}{4} (r+1) (s+1) \left (1-t^2\right )+\\ &x_W\,\frac{1}{4} (r+1) \left (1-s^2\right ) (t+1)+\\ &x_X\,\frac{1}{4} (r+1) (1-s) \left (1-t^2\right )+\\ &x_Y\,\frac{1}{4} \left (1-r^2\right ) (1-s) (1-t)+\\ &x_Z\,\frac{1}{4} \left (1-r^2\right ) (s+1) (1-t)+\\ &x_A\,\frac{1}{4} \left (1-r^2\right ) (s+1) (t+1)+\\ &x_B\,\frac{1}{4} \left (1-r^2\right ) (1-s) (t+1) \end{align*}

Taking partial derivative of the above w.r.t. $$r,s,t$$ in turn gives the following \begin{align*} \frac{\partial{x}}{\partial{r}} = &x_I\left (\frac{1}{8} (s-1) (t-1) (2 r+s+t+1)\right )+\\ &x_J\left (\frac{1}{8} (s+1) (t-1) (-2 r+s-t-1)\right )+\\ &x_K\left (-\frac{1}{8} (s+1) (t+1) (-2 r+s+t-1)\right )+\\ &x_L\left (-\frac{1}{8} (s-1) (t+1) (2 r+s-t+1)\right )+\\ &x_M\left (-\frac{1}{8} (s-1) (t-1) (-2 r+s+t+1)\right )+\\ &x_N\left (-\frac{1}{8} (s+1) (t-1) (2 r+s-t-1)\right )+\\ &x_O\left (\frac{1}{8} (s+1) (t+1) (2 r+s+t-1)\right )+\\ &x_P\left (\frac{1}{8} (s-1) (t+1) (-2 r+s-t+1)\right )+\\ &x_Q\left (-\frac{1}{4} \left (s^2-1\right ) (t-1)\right )+\\ &x_R\left (\frac{1}{4} (s+1) \left (t^2-1\right )\right )+\\ &x_S\left (\frac{1}{4} \left (s^2-1\right ) (t+1)\right )+\\ &x_T\left (-\frac{1}{4} (s-1) \left (t^2-1\right )\right )+\\ &x_U\left (\frac{1}{4} \left (s^2-1\right ) (t-1)\right )+\\ &x_V\left (-\frac{1}{4} (s+1) \left (t^2-1\right )\right )+\\ &x_W\left (-\frac{1}{4} \left (s^2-1\right ) (t+1)\right )+\\ &x_X\left (\frac{1}{4} (s-1) \left (t^2-1\right )\right )+\\ &x_Y\left (-\frac{1}{2} r (s-1) (t-1)\right )+\\ &x_Z\left (\frac{1}{2} r (s+1) (t-1)\right )+\\ &x_A\left (-\frac{1}{2} r (s+1) (t+1)\right )+\\ &x_B\left (\frac{1}{2} r (s-1) (t+1)\right ) \end{align*}

In the Matlab implementation,the terms $$r,s,t$$ in the above expression are the Gaussian integration points along the three directions. Similarly, we now ﬁnd $$\frac{\partial{x}}{\partial{s}}$$ as above. This results in \begin{align*} \frac{\partial{x}}{\partial{s}} = &x_I\left (\frac{1}{8} (r-1) (t-1) (r+2 s+t+1)\right )+\\ &x_J\left (-\frac{1}{8} (r-1) (t-1) (r-2 s+t+1)\right )+\\ &x_K\left (\frac{1}{8} (r-1) (t+1) (r-2 s-t+1)\right )+\\ &x_L\left (-\frac{1}{8} (r-1) (t+1) (r+2 s-t+1)\right )+\\ &x_M\left (\frac{1}{8} (r+1) (t-1) (r-2 s-t-1)\right )+\\ &x_N\left (-\frac{1}{8} (r+1) (t-1) (r+2 s-t-1)\right )+\\ &x_O\left (\frac{1}{8} (r+1) (t+1) (r+2 s+t-1)\right )+\\ &x_P\left (-\frac{1}{8} (r+1) (t+1) (r-2 s+t-1)\right )+\\ &x_Q\left (-\frac{1}{2} (r-1) s (t-1)\right )+\\ &x_R\left (\frac{1}{4} (r-1) \left (t^2-1\right )\right )+\\ &x_S\left (\frac{1}{2} (r-1) s (t+1)\right )+\\ &x_T\left (-\frac{1}{4} (r-1) \left (t^2-1\right )\right )+\\ &x_U\left (\frac{1}{2} (r+1) s (t-1)\right )+\\ &x_V\left (-\frac{1}{4} (r+1) \left (t^2-1\right )\right )+\\ &x_W\left (-\frac{1}{2} (r+1) s (t+1)\right )+\\ &x_X\left (\frac{1}{4} (r+1) \left (t^2-1\right )\right )+\\ &x_Y\left (-\frac{1}{4} \left (r^2-1\right ) (t-1)\right )+\\ &x_Z\left (\frac{1}{4} \left (r^2-1\right ) (t-1)\right )+\\ &x_A\left (-\frac{1}{4} \left (r^2-1\right ) (t+1)\right )+\\ &x_B\left (\frac{1}{4} \left (r^2-1\right ) (t+1)\right ) \end{align*}

Similarly, we now ﬁnd $$\frac{\partial{x}}{\partial{t}}$$ as above. This results in \begin{align*} \frac{\partial{x}}{\partial{t}} = &x_I\left (\frac{1}{8} (r-1) (s-1) (r+s+2 t+1)\right )+\\ &x_J\left (-\frac{1}{8} (r-1) (s+1) (r-s+2 t+1)\right )+\\ &x_K\left (\frac{1}{8} (r-1) (s+1) (r-s-2 t+1)\right )+\\ &x_L\left (-\frac{1}{8} (r-1) (s-1) (r+s-2 t+1)\right )+\\ &x_M\left (\frac{1}{8} (r+1) (s-1) (r-s-2 t-1)\right )+\\ &x_N\left (-\frac{1}{8} (r+1) (s+1) (r+s-2 t-1)\right )+\\ &x_O\left (\frac{1}{8} (r+1) (s+1) (r+s+2 t-1)\right )+\\ &x_P\left (-\frac{1}{8} (r+1) (s-1) (r-s+2 t-1)\right )+\\ &x_Q\left (-\frac{1}{4} (r-1) \left (s^2-1\right )\right )+\\ &x_R\left (\frac{1}{2} (r-1) (s+1) t\right )+\\ &x_S\left (\frac{1}{4} (r-1) \left (s^2-1\right )\right )+\\ &x_T\left (-\frac{1}{2} (r-1) (s-1) t\right )+\\ &x_U\left (\frac{1}{4} (r+1) \left (s^2-1\right )\right )+\\ &x_V\left (-\frac{1}{2} (r+1) (s+1) t\right )+\\ &x_W\left (-\frac{1}{4} (r+1) \left (s^2-1\right )\right )+\\ &x_X\left (\frac{1}{2} (r+1) (s-1) t\right )+\\ &x_Y\left (-\frac{1}{4} \left (r^2-1\right ) (s-1)\right )+\\ &x_Z\left (\frac{1}{4} \left (r^2-1\right ) (s+1)\right )+\\ &x_A\left (-\frac{1}{4} \left (r^2-1\right ) (s+1)\right )+\\ &x_B\left (\frac{1}{4} \left (r^2-1\right ) (s-1)\right ) \end{align*}

##### 2.5.3.3 $$y(r,s,t)$$ terms

We now repeat all the above to ﬁnd $$\frac{\partial{y}}{\partial{r}},\frac{\partial{y}}{\partial{s}}$$ and $$\frac{\partial{y}}{\partial{t}}$$. We ﬁrst need to expand $$y(r,s,t) = \sum _{i=I}^{I=B} y_i f_i(r,s,t)$$, which gives \begin{align*} y_ =&y_I\,f_I+y_J\,f_J+y_K\,f_K+y_L\,f_L+y_M\,f_M+y_N\,f_N+y_O\,f_O+y_P\,f_P+y_Q\,f_Q+y_R\,f_R+y_S\,f_S+y_T\,f_T+y_U\,f_U+y_V\,f_V+y_W\,f_W+y_X\,f_X+y_Y\,f_Y+y_Z\,f_Z+y_A\,f_A+y_B\,f_B \end{align*}

Expanding the above gives \begin{align*} y(r,s,t)=&y_I\,\frac{1}{8} (1-r) (1-s) (1-t) (-r-s-t-2)+\\ &y_J\,\frac{1}{8} (1-r) (s+1) (1-t) (-r+s-t-2)+\\ &y_K\,\frac{1}{8} (1-r) (s+1) (t+1) (-r+s+t-2)+\\ &y_L\,\frac{1}{8} (1-r) (1-s) (t+1) (-r-s+t-2)+\\ &y_M\,\frac{1}{8} (r+1) (1-s) (1-t) (r-s-t-2)+\\ &y_N\,\frac{1}{8} (r+1) (s+1) (1-t) (r+s-t-2)+\\ &y_O\,\frac{1}{8} (r+1) (s+1) (t+1) (r+s+t-2)+\\ &y_P\,\frac{1}{8} (r+1) (1-s) (t+1) (r-s+t-2)+\\ &y_Q\,\frac{1}{4} (1-r) \left (1-s^2\right ) (1-t)+\\ &y_R\,\frac{1}{4} (1-r) (s+1) \left (1-t^2\right )+\\ &y_S\,\frac{1}{4} (1-r) \left (1-s^2\right ) (t+1)+\\ &y_T\,\frac{1}{4} (1-r) (1-s) \left (1-t^2\right )+\\ &y_U\,\frac{1}{4} (r+1) \left (1-s^2\right ) (1-t)+\\ &y_V\,\frac{1}{4} (r+1) (s+1) \left (1-t^2\right )+\\ &y_W\,\frac{1}{4} (r+1) \left (1-s^2\right ) (t+1)+\\ &y_X\,\frac{1}{4} (r+1) (1-s) \left (1-t^2\right )+\\ &y_Y\,\frac{1}{4} \left (1-r^2\right ) (1-s) (1-t)+\\ &y_Z\,\frac{1}{4} \left (1-r^2\right ) (s+1) (1-t)+\\ &y_A\,\frac{1}{4} \left (1-r^2\right ) (s+1) (t+1)+\\ &y_B\,\frac{1}{4} \left (1-r^2\right ) (1-s) (t+1) \end{align*}

Taking partial derivatives of the above w.r.t. $$r,s,t$$ in turn, we see that it gives similar results to earlier ones, but the only diﬀerence is in the multipliers now being the $$y_i$$ values of coordinates instead of the $$x_i$$ coordinates. This is reproduced again for completion \begin{align*} \frac{\partial{y}}{\partial{r}} = &y_I\left (\frac{1}{8} (s-1) (t-1) (2 r+s+t+1)\right )+\\ &y_J\left (\frac{1}{8} (s+1) (t-1) (-2 r+s-t-1)\right )+\\ &y_K\left (-\frac{1}{8} (s+1) (t+1) (-2 r+s+t-1)\right )+\\ &y_L\left (-\frac{1}{8} (s-1) (t+1) (2 r+s-t+1)\right )+\\ &y_M\left (-\frac{1}{8} (s-1) (t-1) (-2 r+s+t+1)\right )+\\ &y_N\left (-\frac{1}{8} (s+1) (t-1) (2 r+s-t-1)\right )+\\ &y_O\left (\frac{1}{8} (s+1) (t+1) (2 r+s+t-1)\right )+\\ &y_P\left (\frac{1}{8} (s-1) (t+1) (-2 r+s-t+1)\right )+\\ &y_Q\left (-\frac{1}{4} \left (s^2-1\right ) (t-1)\right )+\\ &y_R\left (\frac{1}{4} (s+1) \left (t^2-1\right )\right )+\\ &y_S\left (\frac{1}{4} \left (s^2-1\right ) (t+1)\right )+\\ &y_T\left (-\frac{1}{4} (s-1) \left (t^2-1\right )\right )+\\ &y_U\left (\frac{1}{4} \left (s^2-1\right ) (t-1)\right )+\\ &y_V\left (-\frac{1}{4} (s+1) \left (t^2-1\right )\right )+\\ &y_W\left (-\frac{1}{4} \left (s^2-1\right ) (t+1)\right )+\\ &y_X\left (\frac{1}{4} (s-1) \left (t^2-1\right )\right )+\\ &y_Y\left (-\frac{1}{2} r (s-1) (t-1)\right )+\\ &y_Z\left (\frac{1}{2} r (s+1) (t-1)\right )+\\ &y_A\left (-\frac{1}{2} r (s+1) (t+1)\right )+\\ &y_B\left (\frac{1}{2} r (s-1) (t+1)\right ) \end{align*}

Similarly, we now ﬁnd $$\frac{\partial{y}}{\partial{s}}$$ as above. This results in \begin{align*} \frac{\partial{y}}{\partial{s}} = &y_I\left (\frac{1}{8} (r-1) (t-1) (r+2 s+t+1)\right )+\\ &y_J\left (-\frac{1}{8} (r-1) (t-1) (r-2 s+t+1)\right )+\\ &y_K\left (\frac{1}{8} (r-1) (t+1) (r-2 s-t+1)\right )+\\ &y_L\left (-\frac{1}{8} (r-1) (t+1) (r+2 s-t+1)\right )+\\ &y_M\left (\frac{1}{8} (r+1) (t-1) (r-2 s-t-1)\right )+\\ &y_N\left (-\frac{1}{8} (r+1) (t-1) (r+2 s-t-1)\right )+\\ &y_O\left (\frac{1}{8} (r+1) (t+1) (r+2 s+t-1)\right )+\\ &y_P\left (-\frac{1}{8} (r+1) (t+1) (r-2 s+t-1)\right )+\\ &y_Q\left (-\frac{1}{2} (r-1) s (t-1)\right )+\\ &y_R\left (\frac{1}{4} (r-1) \left (t^2-1\right )\right )+\\ &y_S\left (\frac{1}{2} (r-1) s (t+1)\right )+\\ &y_T\left (-\frac{1}{4} (r-1) \left (t^2-1\right )\right )+\\ &y_U\left (\frac{1}{2} (r+1) s (t-1)\right )+\\ &y_V\left (-\frac{1}{4} (r+1) \left (t^2-1\right )\right )+\\ &y_W\left (-\frac{1}{2} (r+1) s (t+1)\right )+\\ &y_X\left (\frac{1}{4} (r+1) \left (t^2-1\right )\right )+\\ &y_Y\left (-\frac{1}{4} \left (r^2-1\right ) (t-1)\right )+\\ &y_Z\left (\frac{1}{4} \left (r^2-1\right ) (t-1)\right )+\\ &y_A\left (-\frac{1}{4} \left (r^2-1\right ) (t+1)\right )+\\ &y_B\left (\frac{1}{4} \left (r^2-1\right ) (t+1)\right ) \end{align*}

We now ﬁnd $$\frac{\partial{y}}{\partial{t}}$$ as above. This results in \begin{align*} \frac{\partial{y}}{\partial{t}} = &y_I\left (\frac{1}{8} (r-1) (s-1) (r+s+2 t+1)\right )+\\ &y_J\left (-\frac{1}{8} (r-1) (s+1) (r-s+2 t+1)\right )+\\ &y_K\left (\frac{1}{8} (r-1) (s+1) (r-s-2 t+1)\right )+\\ &y_L\left (-\frac{1}{8} (r-1) (s-1) (r+s-2 t+1)\right )+\\ &y_M\left (\frac{1}{8} (r+1) (s-1) (r-s-2 t-1)\right )+\\ &y_N\left (-\frac{1}{8} (r+1) (s+1) (r+s-2 t-1)\right )+\\ &y_O\left (\frac{1}{8} (r+1) (s+1) (r+s+2 t-1)\right )+\\ &y_P\left (-\frac{1}{8} (r+1) (s-1) (r-s+2 t-1)\right )+\\ &y_Q\left (-\frac{1}{4} (r-1) \left (s^2-1\right )\right )+\\ &y_R\left (\frac{1}{2} (r-1) (s+1) t\right )+\\ &y_S\left (\frac{1}{4} (r-1) \left (s^2-1\right )\right )+\\ &y_T\left (-\frac{1}{2} (r-1) (s-1) t\right )+\\ &y_U\left (\frac{1}{4} (r+1) \left (s^2-1\right )\right )+\\ &y_V\left (-\frac{1}{2} (r+1) (s+1) t\right )+\\ &y_W\left (-\frac{1}{4} (r+1) \left (s^2-1\right )\right )+\\ &y_X\left (\frac{1}{2} (r+1) (s-1) t\right )+\\ &y_Y\left (-\frac{1}{4} \left (r^2-1\right ) (s-1)\right )+\\ &y_Z\left (\frac{1}{4} \left (r^2-1\right ) (s+1)\right )+\\ &y_A\left (-\frac{1}{4} \left (r^2-1\right ) (s+1)\right )+\\ &y_B\left (\frac{1}{4} \left (r^2-1\right ) (s-1)\right ) \end{align*}

##### 2.5.3.4 $$z(r,s,t)$$ terms

We now repeat all the above to ﬁnd $$\frac{\partial{z}}{\partial{r}},\frac{\partial{z}}{\partial{s}}$$ and $$\frac{\partial{z}}{\partial{t}}$$. These produce similar results to the above, but will have $$z_i$$ as multipliers. We ﬁrst need to expand $$z(r,s,t) = \sum _{i=I}^{I=B} z_i f_i(r,s,t)$$, which gives \begin{align*} z_ =&z_I\,f_I+z_J\,f_J+z_K\,f_K+z_L\,f_L+z_M\,f_M+z_N\,f_N+z_O\,f_O+z_P\,f_P+z_Q\,f_Q+z_R\,f_R+z_S\,f_S+z_T\,f_T+z_U\,f_U+z_V\,f_V+z_W\,f_W+z_X\,f_X+z_Y\,f_Y+z_Z\,f_Z+z_A\,f_A+z_B\,f_B \end{align*}

Expanding the above gives \begin{align*} z(r,s,t)=&z_I\,\frac{1}{8} (1-r) (1-s) (1-t) (-r-s-t-2)+\\ &z_J\,\frac{1}{8} (1-r) (s+1) (1-t) (-r+s-t-2)+\\ &z_K\,\frac{1}{8} (1-r) (s+1) (t+1) (-r+s+t-2)+\\ &z_L\,\frac{1}{8} (1-r) (1-s) (t+1) (-r-s+t-2)+\\ &z_M\,\frac{1}{8} (r+1) (1-s) (1-t) (r-s-t-2)+\\ &z_N\,\frac{1}{8} (r+1) (s+1) (1-t) (r+s-t-2)+\\ &z_O\,\frac{1}{8} (r+1) (s+1) (t+1) (r+s+t-2)+\\ &z_P\,\frac{1}{8} (r+1) (1-s) (t+1) (r-s+t-2)+\\ &z_Q\,\frac{1}{4} (1-r) \left (1-s^2\right ) (1-t)+\\ &z_R\,\frac{1}{4} (1-r) (s+1) \left (1-t^2\right )+\\ &z_S\,\frac{1}{4} (1-r) \left (1-s^2\right ) (t+1)+\\ &z_T\,\frac{1}{4} (1-r) (1-s) \left (1-t^2\right )+\\ &z_U\,\frac{1}{4} (r+1) \left (1-s^2\right ) (1-t)+\\ &z_V\,\frac{1}{4} (r+1) (s+1) \left (1-t^2\right )+\\ &z_W\,\frac{1}{4} (r+1) \left (1-s^2\right ) (t+1)+\\ &z_X\,\frac{1}{4} (r+1) (1-s) \left (1-t^2\right )+\\ &z_Y\,\frac{1}{4} \left (1-r^2\right ) (1-s) (1-t)+\\ &z_Z\,\frac{1}{4} \left (1-r^2\right ) (s+1) (1-t)+\\ &z_A\,\frac{1}{4} \left (1-r^2\right ) (s+1) (t+1)+\\ &z_B\,\frac{1}{4} \left (1-r^2\right ) (1-s) (t+1) \end{align*}

Taking partial derivative of the above w.r.t. $$r,s,t$$ in turns gives the following \begin{align*} \frac{\partial{z}}{\partial{r}} = &z_I\left (\frac{1}{8} (s-1) (t-1) (2 r+s+t+1)\right )+\\ &z_J\left (\frac{1}{8} (s+1) (t-1) (-2 r+s-t-1)\right )+\\ &z_K\left (-\frac{1}{8} (s+1) (t+1) (-2 r+s+t-1)\right )+\\ &z_L\left (-\frac{1}{8} (s-1) (t+1) (2 r+s-t+1)\right )+\\ &z_M\left (-\frac{1}{8} (s-1) (t-1) (-2 r+s+t+1)\right )+\\ &z_N\left (-\frac{1}{8} (s+1) (t-1) (2 r+s-t-1)\right )+\\ &z_O\left (\frac{1}{8} (s+1) (t+1) (2 r+s+t-1)\right )+\\ &z_P\left (\frac{1}{8} (s-1) (t+1) (-2 r+s-t+1)\right )+\\ &z_Q\left (-\frac{1}{4} \left (s^2-1\right ) (t-1)\right )+\\ &z_R\left (\frac{1}{4} (s+1) \left (t^2-1\right )\right )+\\ &z_S\left (\frac{1}{4} \left (s^2-1\right ) (t+1)\right )+\\ &z_T\left (-\frac{1}{4} (s-1) \left (t^2-1\right )\right )+\\ &z_U\left (\frac{1}{4} \left (s^2-1\right ) (t-1)\right )+\\ &z_V\left (-\frac{1}{4} (s+1) \left (t^2-1\right )\right )+\\ &z_W\left (-\frac{1}{4} \left (s^2-1\right ) (t+1)\right )+\\ &z_X\left (\frac{1}{4} (s-1) \left (t^2-1\right )\right )+\\ &z_Y\left (-\frac{1}{2} r (s-1) (t-1)\right )+\\ &z_Z\left (\frac{1}{2} r (s+1) (t-1)\right )+\\ &z_A\left (-\frac{1}{2} r (s+1) (t+1)\right )+\\ &z_B\left (\frac{1}{2} r (s-1) (t+1)\right ) \end{align*}

Similarly, $$\frac{\partial{z}}{\partial{s}}$$ results in \begin{align*} \frac{\partial{z}}{\partial{s}} = &z_I\left (\frac{1}{8} (r-1) (t-1) (r+2 s+t+1)\right )+\\ &z_J\left (-\frac{1}{8} (r-1) (t-1) (r-2 s+t+1)\right )+\\ &z_K\left (\frac{1}{8} (r-1) (t+1) (r-2 s-t+1)\right )+\\ &z_L\left (-\frac{1}{8} (r-1) (t+1) (r+2 s-t+1)\right )+\\ &z_M\left (\frac{1}{8} (r+1) (t-1) (r-2 s-t-1)\right )+\\ &z_N\left (-\frac{1}{8} (r+1) (t-1) (r+2 s-t-1)\right )+\\ &z_O\left (\frac{1}{8} (r+1) (t+1) (r+2 s+t-1)\right )+\\ &z_P\left (-\frac{1}{8} (r+1) (t+1) (r-2 s+t-1)\right )+\\ &z_Q\left (-\frac{1}{2} (r-1) s (t-1)\right )+\\ &z_R\left (\frac{1}{4} (r-1) \left (t^2-1\right )\right )+\\ &z_S\left (\frac{1}{2} (r-1) s (t+1)\right )+\\ &z_T\left (-\frac{1}{4} (r-1) \left (t^2-1\right )\right )+\\ &z_U\left (\frac{1}{2} (r+1) s (t-1)\right )+\\ &z_V\left (-\frac{1}{4} (r+1) \left (t^2-1\right )\right )+\\ &z_W\left (-\frac{1}{2} (r+1) s (t+1)\right )+\\ &z_X\left (\frac{1}{4} (r+1) \left (t^2-1\right )\right )+\\ &z_Y\left (-\frac{1}{4} \left (r^2-1\right ) (t-1)\right )+\\ &z_Z\left (\frac{1}{4} \left (r^2-1\right ) (t-1)\right )+\\ &z_A\left (-\frac{1}{4} \left (r^2-1\right ) (t+1)\right )+\\ &z_B\left (\frac{1}{4} \left (r^2-1\right ) (t+1)\right ) \end{align*}

And $$\frac{\partial{z}}{\partial{t}}$$ gives \begin{align*} \frac{\partial{z}}{\partial{t}} = &z_I\left (\frac{1}{8} (r-1) (s-1) (r+s+2 t+1)\right )+\\ &z_J\left (-\frac{1}{8} (r-1) (s+1) (r-s+2 t+1)\right )+\\ &z_K\left (\frac{1}{8} (r-1) (s+1) (r-s-2 t+1)\right )+\\ &z_L\left (-\frac{1}{8} (r-1) (s-1) (r+s-2 t+1)\right )+\\ &z_M\left (\frac{1}{8} (r+1) (s-1) (r-s-2 t-1)\right )+\\ &z_N\left (-\frac{1}{8} (r+1) (s+1) (r+s-2 t-1)\right )+\\ &z_O\left (\frac{1}{8} (r+1) (s+1) (r+s+2 t-1)\right )+\\ &z_P\left (-\frac{1}{8} (r+1) (s-1) (r-s+2 t-1)\right )+\\ &z_Q\left (-\frac{1}{4} (r-1) \left (s^2-1\right )\right )+\\ &z_R\left (\frac{1}{2} (r-1) (s+1) t\right )+\\ &z_S\left (\frac{1}{4} (r-1) \left (s^2-1\right )\right )+\\ &z_T\left (-\frac{1}{2} (r-1) (s-1) t\right )+\\ &z_U\left (\frac{1}{4} (r+1) \left (s^2-1\right )\right )+\\ &z_V\left (-\frac{1}{2} (r+1) (s+1) t\right )+\\ &z_W\left (-\frac{1}{4} (r+1) \left (s^2-1\right )\right )+\\ &z_X\left (\frac{1}{2} (r+1) (s-1) t\right )+\\ &z_Y\left (-\frac{1}{4} \left (r^2-1\right ) (s-1)\right )+\\ &z_Z\left (\frac{1}{4} \left (r^2-1\right ) (s+1)\right )+\\ &z_A\left (-\frac{1}{4} \left (r^2-1\right ) (s+1)\right )+\\ &z_B\left (\frac{1}{4} \left (r^2-1\right ) (s-1)\right ) \end{align*}

Finally now we can determine the Jacobian and its determinant using the above expressions. This is done in the Matlab code provided. The following Jacobian Matrix is evaluated at each Gaussian integration point then its determinant is found using det() command. \begin{align*} \begin{pmatrix} \frac{\partial{x}}{\partial{r}}&\frac{\partial{y}}{\partial{r}}&\frac{\partial{z}}{\partial{r}}\\ \frac{\partial{x}}{\partial{s}}&\frac{\partial{y}}{\partial{s}}&\frac{\partial{z}}{\partial{s}}\\ \frac{\partial{x}}{\partial{t}}&\frac{\partial{y}}{\partial{t}}&\frac{\partial{z}}{\partial{t}}\\ \end{pmatrix} \end{align*}

##### 2.5.3.5 results

The ﬁrst step was to obtain estimate of the volume in order to verify that the volume calculation was valid and that the Jacobian was correct.

An independent small piece of code was written to plot the 3D shape and obtain its volume using a build-in function in the computer algebra program Mathematica. This is a plot of the 3D shape generated and below it is the code used to generate the plot, with the volume found shown in the title. We now know the volume should be $$0.0042514 \text{cm}^3$$ from the above independent veriﬁcation. The Matlab code was now implemented, and the volume was veriﬁed to be the same. Also, a separate test was run to verify that $$\|J\|=1$$ for a test 3D volume which was aligned along the same orientation as the natural coordinates as shown below. The code used to plot the above is

This test also passed in Matlab and gave a volume of $$8 \text{cm}^3$$ as expected. Here is the small code segment in Matlab which veriﬁes the above.

The output of the above is given below, with the rest of the program output.

The program nma_EMA_471_HW5_problem_3.m implements the main solution to this problem and included in the zip ﬁle. It runs both the Jacobian veriﬁcation and the load calculations after that.

The main loop of the Matlab function iterates over three indices $$i,j,k$$ from 1 to 3 each. In the inner most loop, it ﬁnds the determinant of the Jacobian, the mass density, and evaluates the shape function at the Gaussian points, then sums the result. At the end it prints the work-equivalent conversion for each of the twenty nodes. The following shows the main core of the program

The ﬁnal result is in the following table

 Corner work-equivalent load ($$z$$ direction) $$\frac{\text{gram-cm}}{\text{sec}^2}$$ Newton units I $$-0.1097835$$ $$-0.000001934$$ J $$-0.1336398$$ $$-0.000001990$$ K $$-0.1159590$$ $$-0.000002$$ L $$-0.1143340$$ $$-0.000001928$$ M $$-0.1373666$$ $$-0.000001927$$ N $$-0.1611689$$ $$-0.000001984$$ O $$-0.1418730$$ $$-0.00000198$$ P $$-0.1160910$$ $$-0.000001924$$ Q $$0.1701952$$ $$0.000002614$$ R $$0.1812529$$ $$0.000002739$$ S $$0.0808510$$ $$0.000002592$$ T $$0.0741811$$ $$0.000002513$$ U $$0.2242564$$ $$0.000002589$$ V $$0.2374809$$ $$0.000002728$$ W $$0.1905987$$ $$0.000002589$$ X $$0.1802555$$ $$0.000002476$$ Y $$0.1723161$$ $$0.000002482$$ Z $$0.2334560$$ $$0.000002734$$ A $$0.1929484$$ $$0.000002710$$ B $$0.0986489$$ $$0.000002487$$

We also see that the load on the corners is negative while on the middle nodes it is positive. This agrees with what one would expect as per class notes on the 8-node element. Only diﬀerence is that this is a 3D element.

The following is the console output from running the above program. It is implemented using Matlab 2016a. It starts with the Jacobian veriﬁcation then it will run the main task next only if the veriﬁcation passes.

>>nma_EMA_471_HW5_problem_3()

starting verification of Jacobian....
|J|= 1.000 at Gaussian point [r=-0.775,s=-0.775,t=-0.775]
|J|= 1.000 at Gaussian point [r=0.000,s=-0.775,t=-0.775]
|J|= 1.000 at Gaussian point [r=0.775,s=-0.775,t=-0.775]
|J|= 1.000 at Gaussian point [r=-0.775,s=-0.775,t=0.000]
|J|= 1.000 at Gaussian point [r=0.000,s=-0.775,t=0.000]
|J|= 1.000 at Gaussian point [r=0.775,s=-0.775,t=0.000]
|J|= 1.000 at Gaussian point [r=-0.775,s=-0.775,t=0.775]
|J|= 1.000 at Gaussian point [r=0.000,s=-0.775,t=0.775]
|J|= 1.000 at Gaussian point [r=0.775,s=-0.775,t=0.775]
|J|= 1.000 at Gaussian point [r=-0.775,s=0.000,t=-0.775]
|J|= 1.000 at Gaussian point [r=0.000,s=0.000,t=-0.775]
|J|= 1.000 at Gaussian point [r=0.775,s=0.000,t=-0.775]
|J|= 1.000 at Gaussian point [r=-0.775,s=0.000,t=0.000]
|J|= 1.000 at Gaussian point [r=0.000,s=0.000,t=0.000]
|J|= 1.000 at Gaussian point [r=0.775,s=0.000,t=0.000]
|J|= 1.000 at Gaussian point [r=-0.775,s=0.000,t=0.775]
|J|= 1.000 at Gaussian point [r=0.000,s=0.000,t=0.775]
|J|= 1.000 at Gaussian point [r=0.775,s=0.000,t=0.775]
|J|= 1.000 at Gaussian point [r=-0.775,s=0.775,t=-0.775]
|J|= 1.000 at Gaussian point [r=0.000,s=0.775,t=-0.775]
|J|= 1.000 at Gaussian point [r=0.775,s=0.775,t=-0.775]
|J|= 1.000 at Gaussian point [r=-0.775,s=0.775,t=0.000]
|J|= 1.000 at Gaussian point [r=0.000,s=0.775,t=0.000]
|J|= 1.000 at Gaussian point [r=0.775,s=0.775,t=0.000]
|J|= 1.000 at Gaussian point [r=-0.775,s=0.775,t=0.775]
|J|= 1.000 at Gaussian point [r=0.000,s=0.775,t=0.775]
|J|= 1.000 at Gaussian point [r=0.775,s=0.775,t=0.775]
volume for test is [8.000] (is it 8?)
!! passed Jacobian test. Will run main program now

volume is 0.001426 cm^3

load at corner I = -0.1934002 [gram-cm/sec^2] = -0.000001934 N
load at corner J = -0.1990333 [gram-cm/sec^2] = -0.000001990 N
load at corner K = -0.2000363 [gram-cm/sec^2] = -0.000002000 N
load at corner L = -0.1928213 [gram-cm/sec^2] = -0.000001928 N
load at corner M = -0.1927113 [gram-cm/sec^2] = -0.000001927 N
load at corner N = -0.1984265 [gram-cm/sec^2] = -0.000001984 N
load at corner O = -0.1980087 [gram-cm/sec^2] = -0.000001980 N
load at corner P = -0.1923718 [gram-cm/sec^2] = -0.000001924 N
load at corner Q = 0.2614284 [gram-cm/sec^2] = 0.000002614 N
load at corner R = 0.2739156 [gram-cm/sec^2] = 0.000002739 N
load at corner S = 0.2592383 [gram-cm/sec^2] = 0.000002592 N
load at corner T = 0.2513418 [gram-cm/sec^2] = 0.000002513 N
load at corner U = 0.2589067 [gram-cm/sec^2] = 0.000002589 N
load at corner V = 0.2727979 [gram-cm/sec^2] = 0.000002728 N
load at corner W = 0.2588535 [gram-cm/sec^2] = 0.000002589 N
load at corner X = 0.2475648 [gram-cm/sec^2] = 0.000002476 N
load at corner Y = 0.2481743 [gram-cm/sec^2] = 0.000002482 N
load at corner Z = 0.2733760 [gram-cm/sec^2] = 0.000002734 N
load at corner A = 0.2710102 [gram-cm/sec^2] = 0.000002710 N
load at corner B = 0.2487254 [gram-cm/sec^2] = 0.000002487 N
>>