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 first 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 simplifies 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
The program nma_EMA_471_HW5_problem_2_part_a.m implements the first part of this problem.
The following table shows the result of the computation. It shows the result of the integral using Gaussian quadrature for different 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
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 final result of Gaussian quadrature improved compared to part(a) where one large domain was used. This makes sense. Because we have effectively 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 difference 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\) |
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.
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 find \(\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 find \(\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*}
We now repeat all the above to find \(\frac{\partial{y}}{\partial{r}},\frac{\partial{y}}{\partial{s}}\) and \(\frac{\partial{y}}{\partial{t}}\). We first 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 difference 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 find \(\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 find \(\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*}
We now repeat all the above to find \(\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 first 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*}
The first 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 verification. The Matlab code was now implemented, and the volume was verified 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 verifies 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 file. It runs both the Jacobian verification 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 finds 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 final 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 difference 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 verification then it will run the main task next only if the verification passes.