text:
# Problem 5.1: Deriving the Weak Form From the Strong Form (Data 1b)

## Given: Strong Form, Boundary Conditions, and Basis Functions

The Strong Form from lecture 9-1 and Fish and Belytschko p. 73:

 $\frac{\partial }{{\partial x}}\left[ {J G \frac{{\partial \phi}}{{\partial x}}} \right] + m = 0$ (1.1)

 ${\text{JG}} = \left( {{\text{2}} + {\text{3x}}} \right)$ (1.2)

 ${\text{m}}\left( {\text{x}} \right) = 5x$ (1.3)

 ${\Gamma _g} = \left\{ 0 \right\} \to u(x = 0) = 4$ (1.4)

 ${\Gamma _g} = \left\{ 1 \right\} \to \frac{{du}}{{dx}}\left( {x = 1} \right) = - 6$ (1.5)

## Find: The Weak Form in Order to Solve for $u^h(x)$ with Following Basis Functions

Solve for $\displaystyle u(x)$ and $\displaystyle u^h(x)$ using the following basis functions until the error $\displaystyle e(n) =|u^h(0.5)-u(0.5)| \leq 10^{-6}$.

1. Polynomial Basis Function $\displaystyle \mathcal{F}_P :=\{ x^j,j=0,1,2,...\}$
2. Exponential Basis Function $\displaystyle \mathcal{F}_e :=\{ e^{jx},j=0,1,2,...\}$
3. Fourier Basis Function $\qquad \displaystyle \mathcal{F}_F :=\{cos(jx),j=0,1,2,...,sin(kx),k=1,2,...\}$

On the domain defined by,

 $\Omega = \left] {\alpha ,\beta } \right[ = \left] { 0,1} \right[$

## Solution: Using Boundary Conditions and Constraint Breaking Solutions

### 1. Approximating Solution Using Polynomial Basis Function

Here is the matlab code used to obtain the approximating solutions: Matlab Code

 $\frac{\partial}{\partial x} \left[ {\left( {{\text{2}} + {\text{3x}}} \right) \frac{{\partial u}}{{\partial x}}} \right] + 5x = 0$ (1.21)

Multiplying Equ 2.20 by w(x) and integrating within the domain:

 $\displaystyle \int_0^{1} {w(x) \cdot \left[\frac{\partial}{\partial x}\left [ (2+3x)\frac{\partial u}{\partial x} \right ]+5x\right] } dx = 0$ (1.22)

Next we will integrate by parts to obtain,

 $\displaystyle \left. {w(x)(2+3x)\frac{\partial u}{\partial x}} \right|_0^{1} - \int_0^{1} {\frac{\partial w}{\partial x}(2+3x)\frac{\partial u}{\partial x}} dx + \int_0^{1} 5wx dx= 0$ (1.23)

Since $\displaystyle w(x=0) = 0$ we can simplify (2.22) in the following way to obtain the weak form,

 $-30w(1) - \int_0^1 {\left[ {\frac{{dw}}{{dx}}\left( {2 + 3x} \right)\frac{{du}}{{dx}}} \right]} dx + \int_0^1 {5xw(x)} dx = 0$ (1.24)

Now we can use the polynomial basis functions,

 $w(x) = {b_i}{(x - 1)^i} \to \frac{{dw}}{{dx}} = i{b_i}{(x - 1)^{i - 1}}$ (1.25)
 $\displaystyle u(x) = {b_j}{(x - 1)^j} \to \frac{{du}}{{dx}}j{b_j}{(x - 1)^{j - 1}}$ (1.26)
 $\displaystyle such \quad that \quad w(0) = 0$ (1.27)

Which allows us to write the weak form (2.25) in terms of our matrices $\displaystyle {\color{red}\bold{K}_{FF}}$, $\displaystyle {\color{red}\bold{d}_{F}}$ and $\displaystyle {\color{red}\bold{F}_{F}}$.

 $\underbrace{\left[ \int_0^1 {{\frac{{\partial b_i(x)}}{{\partial x}} (2 + 3x)\frac{{\partial b_j(x)}}{{\partial x}}} } dx \right]}_{{\color{red}\bold{K}_{FF}}} \underbrace{d_i}_{{\color{red}\bold{d}_{F}}} = \underbrace{\left[ \int_0^1 {5xb_i(x)} dx - 30b_i(x)|_{x=1} \right]}_{{\color{red}\bold{F}_{F}}}$ (1.28)

Next we wish to obtain the solution to $\displaystyle \bold{d_F}$ by solving the following system of equations,

 $\displaystyle \left[ \bold{K_{FF}} \right]_{{\color{red}n \times n}} \left[ \bold{d_{F}} \right]_{{\color{red}n \times 1}}= \left[ \bold{F_{F}} \right]_{{\color{red}n \times 1}}$ (1.29)

Solving for the approximating solution, the d matrix can transposed and multiplied by the basis function matrix:

 $U_p^h = {d_0}{b_0} + {d_1}{b_1} + {d_2}{b_2}$ (1.30)

Using the provided matlab code to obtain the d coefficients and applying to the basis function:

#### n=8

We found that for $\displaystyle n=8$ the error was of the order $\displaystyle 10^{-6}$. Using the Matlab code and rounding off to the nearest hundredths place gives the matrices $\displaystyle \left[ \bold{K_{FF}} \right]_{{\color{red}n \times n}}$ and $\displaystyle \left[ \bold{F_{F}} \right]_{{\color{red}n \times 1}}$

 $\displaystyle \left[ \bold{K_{FF}} \right]_{{\color{red}n \times n}}=\left( \begin{matrix} 7/2 & 4 & 17/4 & 22/5 & 9/2 & 32/7 & 37/8 & 14/3\\ 4 & 17/3 & 33/5 & 36/5 & 160/21 & 111/14 & 49/6 & 376/45\\ 17/4 & 33/5 & 81/10 & 64/7 & 555/56 & 21/2 & 329/30 & 624/55\\ 22/5 & 36/5 & 64/7 & 74/7 & 35/3 & 188/15 & 728/55 & 152/11\\ 9/2 & 160/21 & 555/56 & 35/3 & 235/18 & 156/11 & 665/44 & 620/39\\ 32/7 & 111/14 & 21/2 & 188/15 & 156/11 & 171/11 & 217/13 & 1608/91\\ 37/8 & 49/6 & 329/30 & 728/55 & 665/44 & 217/13 & 469/26 & 96/5\\ 14/3 & 376/45 & 624/55 & 152/11 & 620/39 & 1608/91 & 96/5 & 308/15\\ \end{matrix} \right)$

 $\displaystyle \left[ \bold{F_{F}} \right]_{{\color{red}n \times 1}}={{\left( \begin{matrix} 41/3 \\ 53/4 \\ 13 \\ 77/6 \\ 89/7 \\ 101/8 \\ 113/9 \\ 25/2\\ \\ \end{matrix} \right)}^{T}}$

Using our Matlab code we obtain the following solution for $\displaystyle \left[ \bold{d_{F}} \right]_{{\color{red}n \times 1}}$. We have chosen to report these values in the long format to demonstrate that rounding error has not occurred.

 $\displaystyle \left[ \bold{d_{F}} \right]_{{\color{red}n \times 1}}= {{\left( \begin{matrix} & 4.0000000 \\ & 7.2497811 \\ & -5.4294546 \\ & 4.9208409 \\ & -5.0302374 \\ & 4.5191594 \\ & -3.0192656 \\ & 1.2522298 \\ & -0.23479307 \\ \end{matrix} \right)}}$

 $U_p^h(x) = 4 - 7.2497811(x - 1) + -5.4294546{(x - 1)^2} +4.9208409{(x - 1)^3} + -5.0302374{(x - 1)^4} + 4.5191594{(x - 1)^5} + -3.0192656{(x - 1)^6} 1.2522298 (x - 1)^7-0.23479307(x - 1)^8$

(1.31)

### Matlab Code

#### Plot: Comparing the exact solution $u(x)$ to $u_h(x)$ for $n=2,4,6,8$

The plots for solutions and error at x=.5 are as follows:

#### Plot: Showing the error $e(n) = |u(x)-u_h(x)|$

For the polynomial approximate solution, with the increasing number of nodes, the error decreases. For this particular instance, the error XXXXXXX

### 2. Approximating Soultion Using Fourier Basis Function

Here is the matlab code used to obtain the approximating solutions: Matlab Code

Using the same process as the polynomial solution and using the general equation from equation 2.21, the fourier solution is obtained the same way. W(x) is now equation 5.35 from HW4

Using the provided matlab code to obtain the d coefficients and applying to the basis function:

#### n=8

Using the Matlab code and rounding off to the nearest hundredths place gives the matrices $\displaystyle \left[ \bold{K_{FF}} \right]_{{\color{red}n \times n}}$ and $\displaystyle \left[ \bold{F_{F}} \right]_{{\color{red}n \times 1}}$

 $\displaystyle \left[ \bold{K_{FF}} \right]_{{\color{red}n \times n}}=\left( \begin{matrix} 1.1444 & -1.3612 & 3.2564 & -0.5408 & 3.7008 & 2.3208 & 1.2335 & 4.7249 \\ -1.3612 & 2.3556 & -4.1866 & 2.4001 & -5.8462 & -0.0674 & -4.7299 & -2.9460 \\ 3.2564 & -4.1866 & 9.5121 & -2.3503 & 11.6194 & 5.4649 & 5.7892 & 12.8497 \\ -0.5408 & 2.4001 & -2.3503 & 4.4879 & -5.4725 & 5.3501 & -8.9290 & 3.8999 \\ 3.7008 & -5.8462 & 11.6194 & -5.4725 & 16.8127 & 2.2054 & 14.2443 & 12.2128 \\ 2.3208 & -0.0674 & 5.4649 & 5.3501 & 2.2054 & 14.6873 & -9.6620 & 19.6948 \\ 1.2335 & -4.7299 & 5.7892 & -8.9290 & 14.2443 & -9.6620 & 23.4828 & -3.0985 \\ 4.7249 & -2.9460 & 12.8497 & 3.8999 & 12.2128 & 19.6948 & -3.0985 & 32.5172 \\ \end{matrix} \right)$

 $\displaystyle \left[ \bold{F_{F}} \right]_{{\color{red}n \times 1}}={{\left( \begin{matrix} 4 & -6.1075 & 11.6035 & -18.9907 & 13.0886 & -27.2503 & 3.4218 & -23.8065 & -8.5011 \\ \end{matrix} \right)}^{T}}$

Using our Matlab code we obtain the following solution for $\displaystyle \left[ \bold{d_{F}} \right]_{{\color{red}n \times 1}}$. We have chosen to report these values in the long format to demonstrate that rounding error has not occurred.

 $\displaystyle \left[ \bold{d_{F}} \right]_{{\color{red}n \times 1}}= {{\left( \begin{matrix} 4.0000 & 47.4091 & 48.5755 & -3.4993 & -30.6212 & -4.1185 & 7.2053 & 0.9062 & -0.4248 \\ \end{matrix} \right)^T}}$

### Matlab Code

#### Plot: Comparing the exact solution $u(x)$ to $u_h(x)$ for $n=2,4,6,8$

The plots for solutions and error at x=.5 are as follows:

### 3. Approximating Soultion Using Exponential Basis Function

Here is the matlab code used to obtain the approximating solutions: Matlab Code

Using the same process as the polynomial solution and using the general equation from equation 2.21, the exponential solution is obtained the same way. W(x) is now equation 5.41 from HW4

Using the provided matlab code to obtain the d coefficients and applying to the basis function:

#### n=8

Using the Matlab code and rounding off to the nearest hundredths place gives the matrices $\displaystyle \left[ \bold{K_{FF}} \right]_{{\color{red}n \times n}}$ and $\displaystyle \left[ \bold{F_{F}} \right]_{{\color{red}n \times 1}}$

 $\displaystyle \left[ \bold{K_{FF}} \right]_{{\color{red}n \times n}}= 10^8\left( \begin{matrix} 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0001 & 0.0003 \\ 0.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0001 & 0.0002 & 0.0006 & 0.0017\\ 0.0000 & 0.0000 & 0.0000 & 0.0001 & 0.0003 & 0.0008 & 0.0022 & 0.0062\\ 0.0000 & 0.0000 & 0.0001 & 0.0003 & 0.0008 & 0.0025 & 0.0072 & 0.0206\\ 0.0000 & 0.0001 & 0.0003 & 0.0008 & 0.0026 & 0.0077 & 0.0225 & 0.0649\\ 0.0000 & 0.0002 & 0.0008 & 0.0025 & 0.0077 & 0.0232 & 0.0682 & 0.1973\\ 0.0001 & 0.0006 & 0.0022 & 0.0072 & 0.0225 & 0.0682 & 0.2014 & 0.5858\\ 0.0003 & 0.0017 & 0.0062 & 0.0206 & 0.0649 & 0.1973 & 0.5858 & 1.7106\\ \end{matrix} \right)$

 $\displaystyle \left[ \bold{F_{F}} \right]_{{\color{red}n \times 1}}={{10^4\left( \begin{matrix} 4 & 0.0023 & 0.0085 & 0.0249 & 0.0692 & 0.1885 & 0.5107 & 1.3817 & 3.7387 \\ \end{matrix} \right)}^{T}}$

Using our Matlab code we obtain the following solution for $\displaystyle \left[ \bold{d_{F}} \right]_{{\color{red}n \times 1}}$. We have chosen to report these values in the long format to demonstrate that rounding error has not occurred.

 $\displaystyle \left[ \bold{d_{F}} \right]_{{\color{red}n \times 1}}= {{\left( \begin{matrix} & 4.000000000000 \\ & -416.206015822685 \\ & 1802.371470713180 \\ & -4682.097146605770 \\ & 7738.463302468000 \\ & -8239.102282958180 \\ & 5487.180599937650 \\ & -2083.316654302480 \\ & 344.581494129590 \\ \end{matrix} \right)}}$

### Matlab Code

#### Plots

The plots for solutions and error at x=.5 are as follows:

# Problem 5.2: Deriving the Weak Form From the Strong Form (Data 1)

## Given: Strong Form, Boundary Conditions, and Basis Functions

The Strong Form from lecture 9-1 and Fish and Belytschko p. 73:

 $\frac{\partial }{{\partial x}}\left[ {J G \frac{{\partial \phi}}{{\partial x}}} \right] + m = 0$ (2.1)

 ${\text{JG}} = \left( {{\text{2}} + {\text{3x}}} \right)$ (2.2)

 ${\text{m}}\left( {\text{x}} \right) = 5x$ (2.3)

 ${\Gamma _g} = \left\{ 1 \right\} \to u(x = 1) = 4$ (2.4)

 ${\Gamma _g} = \left\{ 0 \right\} \to \frac{{du}}{{dx}}\left( {x = 0} \right) = - 6$ (2.5)

## Find: The Weak Form in Order to Solve for $u^h(x)$ with Following Basis Functions

Solve for $\displaystyle u(x)$ and $\displaystyle u^h(x)$ using the following basis functions until the error $\displaystyle e(n) =|u^h(0.5)-u(0.5)| \leq 10^{-6}$.

1. Polynomial Basis Function $\displaystyle \mathcal{F}_P :=\{ x^j,j=0,1,2,...\}$
2. Exponential Basis Function $\displaystyle \mathcal{F}_e :=\{ e^{jx},j=0,1,2,...\}$
3. Fourier Basis Function $\qquad \displaystyle \mathcal{F}_F :=\{cos(jx),j=0,1,2,...,sin(kx),k=1,2,...\}$

On the domain defined by,

 $\Omega = \left] {\alpha ,\beta } \right[ = \left] { 0,1} \right[$

## Solution: Using Boundary Conditions and Constraint Breaking Solutions

### 1. Approximating Solution Using Polynomial Basis Function

Here is the matlab code used to obtain the approximating solutions: Matlab Code

 $\frac{\partial}{\partial x} \left[ {\left( {{\text{2}} + {\text{3x}}} \right) \frac{{\partial u}}{{\partial x}}} \right] + 5x = 0$ (2.21)

Multiplying Equ 2.20 by w(x) and integrating within the domain:

 $\displaystyle \int_0^{1} {w(x) \cdot \left[\frac{\partial}{\partial x}\left [ (2+3x)\frac{\partial u}{\partial x} \right ]+5x\right] } dx = 0$ (2.22)

Next we will integrate by parts to obtain,

 $\displaystyle \left. {w(x)(2+3x)\frac{\partial u}{\partial x}} \right|_0^{1} - \int_0^{1} {\frac{\partial w}{\partial x}(2+3x)\frac{\partial u}{\partial x}} dx + \int_0^{1} 5wx dx= 0$ (2.23)

Since $\displaystyle w(x=1) = 0$ we can simplify (2.22) in the following way to obtain the weak form,

 $12w(0) - \int_0^1 {\left[ {\frac{{dw}}{{dx}}\left( {2 + 3x} \right)\frac{{du}}{{dx}}} \right]} dx + \int_0^1 {5xw(x)} dx = 0$ (2.24)

Now we can use the polynomial basis functions,

 $w(x) = {b_i}{(x - 1)^i} \to \frac{{dw}}{{dx}} = i{b_i}{(x - 1)^{i - 1}}$ (2.25)
 $\displaystyle u(x) = {b_j}{(x - 1)^j} \to \frac{{du}}{{dx}}j{b_j}{(x - 1)^{j - 1}}$ (2.26)
 $\displaystyle such \quad that \quad w(1) = 0$ (2.27)

Which allows us to write the weak form (2.25) in terms of our matrices $\displaystyle {\color{red}\bold{K}_{FF}}$, $\displaystyle {\color{red}\bold{d}_{F}}$ and $\displaystyle {\color{red}\bold{F}_{F}}$.

 $\underbrace{\left[ \int_0^1 {{\frac{{\partial b_i(x)}}{{\partial x}} (2 + 3x)\frac{{\partial b_j(x)}}{{\partial x}}} } dx \right]}_{{\color{red}\bold{K}_{FF}}} \underbrace{d_i}_{{\color{red}\bold{d}_{F}}} = \underbrace{\left[ \int_0^1 {5xb_i(x)} dx + (12+18x)b_i(x)|_{x=0} \right]}_{{\color{red}\bold{F}_{F}}}$ (2.28)

Next we wish to obtain the solution to $\displaystyle \bold{d_F}$ by solving the following system of equations,

 $\displaystyle \left[ \bold{K_{FF}} \right]_{{\color{red}n \times n}} \left[ \bold{d_{F}} \right]_{{\color{red}n \times 1}}= \left[ \bold{F_{F}} \right]_{{\color{red}n \times 1}}$ (2.29)

Solving for the approximating solution, the d matrix can transposed and multiplied by the basis function matrix:

 $U_p^h = {d_0}{b_0} + {d_1}{b_1} + {d_2}{b_2}$ (2.30)

Using the provided matlab code to obtain the d coefficients and applying to the basis function:

#### n=2

 $U_p^h(x) = 4 - 2.558(x - 1) + 1.293{(x - 1)^2}$ (2.31)

#### n=4

 $U_p^h(x) = 4 - 2.876(x - 1) + .602{(x - 1)^2} + .310{(x - 1)^3} + .698{(x - 1)^4}$ (2.32)

#### n=6

 $U_p^h(x) = 4 - 2.899(x - 1) + .400{(x - 1)^2} - .121{(x - 1)^3} + .672{(x - 1)^4} + .604{(x - 1)^5} + .378{(x - 1)^6}$ (2.33)

#### n=8

We found that for $\displaystyle n=8$ the error was of the order $\displaystyle 10^{-6}$. Using the Matlab code and rounding off to the nearest hundredths place gives the matrices $\displaystyle \left[ \bold{K_{FF}} \right]_{{\color{red}n \times n}}$ and $\displaystyle \left[ \bold{F_{F}} \right]_{{\color{red}n \times 1}}$

 $\displaystyle \left[ \bold{K_{FF}} \right]_{{\color{red}n \times n}}=\left( \begin{matrix} 3.50 & -3.00 & 2.75 & -2.60 & 2.50 & -2.49 & 2.38 & -2.33 \\ -3.00 & 3.67 & -3.90 & 4.00 & -4.05 & 4.07 & -4.08 & 4.09 \\ 2.75 & -3.90 & 4.50 & -4.86 & 5.09 & -5.25 & 5.37 & -5.45 \\ -2.60 & 4.00 & -4.86 & 5.43 & -5.83 & 6.13 & -6.36 & 6.54 \\ 2.50 & -4.05 & 5.09 & -5.83 & 6.39 & -6.82 & 7.16 & -7.44 \\ -2.43 & 4.07 & -5.25 & 6.13 & -6.82 & 7.36 & -7.81 & 8.18 \\ 2.38 & -4.08 & 5.37 & -6.36 & 7.16 & -7.81 & 8.35 & -8.80 \\ -2.33 & 4.09 & -5.45 & 6.54 & -7.44 & 8.18 & -8.80 & 9.33 \\ \end{matrix} \right)$

 $\displaystyle \left[ \bold{F_{F}} \right]_{{\color{red}n \times 1}}={{\left( \begin{matrix} -12.83 & 12.42 & -12.15 & 12.17 & -12.12 & 12.09 & -12.07 & 12.06 \\ \end{matrix} \right)}^{T}}$

Using our Matlab code we obtain the following solution for $\displaystyle \left[ \bold{d_{F}} \right]_{{\color{red}n \times 1}}$. We have chosen to report these values in the long format to demonstrate that rounding error has not occurred.

 $\displaystyle \left[ \bold{d_{F}} \right]_{{\color{red}n \times 1}}= {{\left( \begin{matrix} & 4.000000000000 \\ & -2.899914248271 \\ & 0.373061337131 \\ & -0.279875467206 \\ & 0.324030513018 \\ & 0.438725044951 \\ & 0.810687654673 \\ & 0.613125135156 \\ & 0.229921935886 \\ \end{matrix} \right)}}$

 $U_p^h(x) = 4 - 2.899(x - 1) + .400{(x - 1)^2} - .121{(x - 1)^3} + .672{(x - 1)^4} + .604{(x - 1)^5} + .378{(x - 1)^6}$

(2.34)

#### Plot: Comparing the exact solution $u(x)$ to $u_h(x)$ for $n=2,4,6,8$

The plots for solutions and error at x=.5 are as follows:

#### Plot: Showing the error $e(n) = |u(x)-u_h(x)|$

For the polynomial approximate solution, with the increasing number of nodes, the error decreases. For this particular instance, the error XXXXXXX

### 2. Approximating Soultion Using Fourier Basis Function

Here is the matlab code used to obtain the approximating solutions: Matlab Code

Using the same process as the polynomial solution and using the general equation from equation 2.21, the fourier solution is obtained the same way. W(x) is now equation 5.35 from HW4

Using the provided matlab code to obtain the d coefficients and applying to the basis function:

#### n=2

 $U_f^h(x) = 4 - 4.125(\cos (x - 1) - 1) - 2.297(\sin (x - 1))$ (2.34)

#### n=4

 $U_f^h(x) = 4 - 11.309(\cos (x - 1) - 1) - 2.066(\sin (x - 1))- 2.422(\cos (2x - 2) - 1) - 0.394(\sin (2x - 2))$ (2.35)

etc...

#### n=8

Using the Matlab code and rounding off to the nearest hundredths place gives the matrices $\displaystyle \left[ \bold{K_{FF}} \right]_{{\color{red}n \times n}}$ and $\displaystyle \left[ \bold{F_{F}} \right]_{{\color{red}n \times 1}}$

#### Plot: Comparing the exact solution $u(x)$ to $u_h(x)$ for $n=2,4,6,8$

The plots for solutions and error at x=.5 are as follows:

### 3. Approximating Soultion Using Exponential Basis Function

Here is the matlab code used to obtain the approximating solutions: Matlab Code

Using the same process as the polynomial solution and using the general equation from equation 2.21, the exponential solution is obtained the same way. W(x) is now equation 5.41 from HW4

Using the provided matlab code to obtain the d coefficients and applying to the basis function:

#### n=2

 $U_e^h(x) = 4 - 14.678({\operatorname{e} ^{(x - 1)}} - 1) + 6.434({e^{(2x - 2)}} - 1)$ (2.36)

#### n=4

 $U_e^h(x) = 4 - 58.869({\operatorname{e} ^{(x - 1)}} - 1) + 93.9({e^{(2x - 2)}} - 1) - 72.659({e^{(3x - 3)}} - 1) + 21.596({e^{(4x - 4)}} - 1)$ (2.37)

#### n=8

Using the Matlab code and rounding off to the nearest hundredths place gives the matrices $\displaystyle \left[ \bold{K_{FF}} \right]_{{\color{red}n \times n}}$ and $\displaystyle \left[ \bold{F_{F}} \right]_{{\color{red}n \times 1}}$

 $\displaystyle \left[ \bold{K_{FF}} \right]_{{\color{red}n \times n}}=\left( \begin{matrix} 1.72 & 2.63 & 3.17 & 3.51 & 3.75 & 3.92 & 4.05 & 4.15 \\ 2.63 & 4.23 & 5.27 & 6.00 & 6.53 & 6.94 & 7.26 & 7.52 \\ 3.17 & 5.27 & 6.74 & 7.83 & 8.67 & 9.33 & 9.87 & 10.31 \\ 3.51 & 6.00 & 7.83 & 9.25 & 10.37 & 11.28 & 12.03 & 12.67 \\ 3.75 & 6.53 & 8.67 & 10.37 & 11.75 & 12.89 & 13.85 & 14.67 \\ 3.92 & 6.94 & 9.33 & 11.28 & 12.89 & 14.25 & 15.41 & 16.41 \\ 4.05 & 7.26 & 9.87 & 12.03 & 13.85 & 15.41 & 16.75 & 17.92 \\ 4.15 & 7.52 & 10.31 & 12.67 & 14.67 & 16.41 & 17.92 & 19.25 \\ \end{matrix} \right)$

 $\displaystyle \left[ \bold{F_{F}} \right]_{{\color{red}n \times 1}}={{\left( \begin{matrix} -8.25 & -11.46 & -12.76 & -13.34 & -13.62 & -13.78 & -13.88 & -13.95 \\ \end{matrix} \right)}^{T}}$

Using our Matlab code we obtain the following solution for $\displaystyle \left[ \bold{d_{F}} \right]_{{\color{red}n \times 1}}$. We have chosen to report these values in the long format to demonstrate that rounding error has not occurred.

 $\displaystyle \left[ \bold{d_{F}} \right]_{{\color{red}n \times 1}}= {{\left( \begin{matrix} & 4.000000000000 \\ & -416.206015822685 \\ & 1802.371470713180 \\ & -4682.097146605770 \\ & 7738.463302468000 \\ & -8239.102282958180 \\ & 5487.180599937650 \\ & -2083.316654302480 \\ & 344.581494129590 \\ \end{matrix} \right)}}$

The plots for solutions and error at x=.5 are as follows:

# Problem 5.3: Torsion of a Bar

Modified version of HW 4.4 defined on page 2 of Lecture 21 from Fish & Belytschko Page 73 Problem 3.7

## Given

 $\frac{\partial }{{\partial x}}\left[ {J G \frac{{\partial \phi}}{{\partial x}}} \right] + m = 0$ (5.3.1)

 ${\text{JG}} = \left( {{\text{2}} + {\text{3x}}} \right)$ (5.3.2)

 ${\text{m}}\left( {\text{x}} \right) = 5x$ (5.3.3)
 ${\Gamma _g} = \left\{ 0 \right\} \to u(x = 0) = 4$ (5.3.4)

 ${\Gamma _g} = \left\{ 1 \right\} \to \frac{{du}}{{dx}}\left( {x = 1} \right) = - 6$ (5.3.5)

Where the Weak Form for a bar in torsion can be defined as

 $\int\limits_{0}^{1} \frac{\partial w}{\partial x} a_2(x) \frac{\partial \phi}{{\partial x}}dx = w(1) h + \int\limits_{0}^{1} w(x)f(x)dx ~~~~~~~\forall w ~ with ~ w(0)=0 ~ and ~ u(x=0) = 4$ (5.3.2)

Make use of the Lagrange Interpolation Basis Function (LIBF) with uniform discretization (equidistant nodes) using $\displaystyle m = 4,6,8$ where the LIBF can defined as follows

 $L_{i,m}(x) = \prod_{\begin{smallmatrix}{k=1}\\ k\neq i\end{smallmatrix}}^{m} \frac{x-x_k}{x_i-x_k} = \frac{(x-x_1)}{(x_i-x_1)} \cdots \frac{(x-x_{k-1})}{(x_i-x_{k-1})} \frac{(x-x_{k})}{(x_i-x_{k})}$. (5.3.3)

## Find

a) Explain how Langrange Interpolation Basis Functions (LIBF) are used as Constraint Breaking Solutions
b) Plot all LIBF used
c) Use matlab quad function or WolframAlpha.com to integrate
d) Plot approximate vs exact solution and convergence error versus equidistant nodes (m)

## Solution

### (a)

Langrange Interpolation Basis Functions satisfy the condition of the constraint breaking solution by assuming a value of zero at every node except for the node where the LIBF is defined there it assumes a value of 1. It acts like the Kronecker Delta

 $L_{i,m}(x) = \begin{cases} 1, & \text{if } x = x_i \\ 0, & \text{if } x = x_k \cdots x_m \end{cases}$ (5.3.4)

### (b)

The following is the LIBF for the $\displaystyle m = 4$

 Figure 1: LIBF m = 4([1])

The following is the LIBF for the $\displaystyle m = 6$

 Figure 1: LIBF m = 6([2])

The following is the LIBF for the $\displaystyle m = 8$

 Figure 1: LIBF m = 8([3])

### (c)

#### Plot: Comparing the exact solution $u(x)$ to $u_h(x)$ for $n=2,4,6,8$

The plots for solutions and error at x=.5 are as follows:

# Problem 5.4: Torsion of a Bar

Modified version of HW 4.4 defined on page 2 of Lecture 21 from Fish & Belytschko Page 73 Problem 3.7

## Given

 $\frac{\partial }{{\partial x}}\left[ {J G \frac{{\partial \phi}}{{\partial x}}} \right] + m = 0$ (5.3.1)

 ${\text{JG}} = \left( {{\text{2}} + {\text{3x}}} \right)$ (5.3.2)

 ${\text{m}}\left( {\text{x}} \right) = 5x$ (5.3.3)
 ${\Gamma _g} = \left\{ 0 \right\} \to u(x = 0) = 4$ (5.3.4)

 ${\Gamma _g} = \left\{ 1 \right\} \to \frac{{du}}{{dx}}\left( {x = 1} \right) = - 6$ (5.3.5)

Where the Weak Form for a bar in torsion can be defined as

 $\int\limits_{0}^{1} \frac{\partial w}{\partial x} a_2(x) \frac{\partial \phi}{{\partial x}}dx = w(1) h + \int\limits_{0}^{1} w(x)f(x)dx ~~~~~~~\forall w ~ with ~ w(0)=0 ~ and ~ u(x=0) = 4$ (5.3.6)

Make use of the Lagrange Interpolation Basis Function (LIBF) with uniform discretization (equidistant nodes) using $\displaystyle m = 4,6,8$ where the LIBF can defined as follows

 $L_{i,m}(x) = \prod_{\begin{smallmatrix}{k=1}\\ k\neq i\end{smallmatrix}}^{m} \frac{x-x_k}{x_i-x_k} = \frac{(x-x_1)}{(x_i-x_1)} \cdots \frac{(x-x_{k-1})}{(x_i-x_{k-1})} \frac{(x-x_{k})}{(x_i-x_{k})}$. (5.3.6)

## Find

a) Explain how Langrange Interpolation Basis Functions (LIBF) are used as Constraint Breaking Solutions
b) Plot all LIBF used
c) Use matlab quad function or WolframAlpha.com to integrate
d) Plot approximate vs exact solution and convergence error versus equidistant nodes (m)

## Solution

### (a)

Langrange Interpolation Basis Functions satisfy the condition of the constraint breaking solution by assuming a value of zero at every node except for the node where the LIBF is defined there it assumes a value of 1. It acts like the Kronecker Delta

 $L_{i,m}(x) = \begin{cases} 1, & \text{if } x = x_i \\ 0, & \text{if } x = x_k \cdots x_m \end{cases}$ (5.4.4)

### (b)

The following is the LIBF for the $\displaystyle m = 4$

 Figure 1: LIBF m = 4([4])

The following is the LIBF for the $\displaystyle m = 6$

 Figure 1: LIBF m = 6([5])

The following is the LIBF for the $\displaystyle m = 8$

 Figure 1: LIBF m = 8([6])

### (c)

Our solution $u \left ( x \right )$ and our weight functions can be approximated using the following definitions,

 \displaystyle \begin{align} & u^{h}(x)=\sum_{j} d_{j} b_{j} \\ &\\ & w^{h}(x)=\sum_{i} c_{i} b_{i} \\ \end{align} $(5.3.3)$

To approximate our solution u, we divide our x domain into m nodal subdivisions. Additionally, we will use m LIBF to approximate our solution. These LIBF will have the following form,

 \displaystyle \begin{align} & L_{i,m}(x)=\prod_{k=1}^{m}\frac{(x-x_{k})}{(x_{i}-x_{k})} \\ \end{align} $(5.4.5)$

Therefore, for m=4, we take the first 4 Lagrange Interpolation Basis Functions, or more specifically,

 \displaystyle \begin{align} & L_{1,4}=\frac{(x-x_{2})(x-x_{3})(x-x_{4})}{(x_{1}-x_{2})(x_{1}-x_{3})(x_{1}-x_{4})} \\ \end{align} $(5.4.6)$
 \displaystyle \begin{align} & L_{2,4}=\frac{(x-x_{1})(x-x_{3})(x-x_{4})}{(x_{2}-x_{1})(x_{2}-x_{3})(x_{2}-x_{4})} \\ \end{align} $(5.4.7)$
 \displaystyle \begin{align} & L_{3,4}=\frac{(x-x_{1})(x-x_{2})(x-x_{4})}{(x_{3}-x_{1})(x_{3}-x_{2})(x_{3}-x_{4})} \\ \end{align} $(5.4.8)$
 \displaystyle \begin{align} & L_{4,4}=\frac{(x-x_{1})(x-x_{2})(x-x_{3})}{(x_{4}-x_{1})(x_{4}-x_{2})(x_{4}-x_{3})}\\ \end{align} $(5.4.9)$

Here we want to develop equations for $\tilde{K_{KK}}$ and $\tilde{F_K}$ and solve for $\tilde{d_{K}}$, similar to the previous problems.

### (d)

#### Plot: Comparing the exact solution $u(x)$ to $u_h(x)$ for $n=2,4,6,8$

The plots for solutions and error at x=.5 are as follows:

# Problem 5.5 : Calculix-Using CCX module

## Disk Problem

After building disc in cgx module then we might want to know that properties of nodes,lines, surfaces, and bodies so now lets look at how we can get more information what we have done for our model.

First we have to open our program. To do this run the following ;

 $\displaystyle [All Programs]\Rightarrow [bConverged]\Rightarrow [Calculix]\Rightarrow [Calculix Command]$

Open calculix and go to directory where the basic examples stays in your harddrive.

 $\displaystyle ...\Rightarrow [Program Files]\Rightarrow [bConverged] \Rightarrow [Calculix] \Rightarrow [Cgx] \Rightarrow [Examples] \Rightarrow [Basic]$

Then type to open disc example:

 $\displaystyle [cgx]\Rightarrow [disc]$

Now we want to learn about point properties. To see point properties in our graph secren type

 $\displaystyle \begin{matrix} plot & pa & all \end{matrix}$

We need to mesh disc to see the mesh properties to mesh the disc please type following steps

 $\displaystyle \begin{matrix} elty & all & qu8 \end{matrix}$
 $\displaystyle \begin{matrix} mesh & all \end{matrix}$

What we have done so far is that we opened disc and meshed it with quadratic elements. So we are ready to see element properties on our graph screen. To see properties of elements type,

 $\displaystyle \begin{matrix} plot & ma & all \end{matrix}$

Also we can get more specific properties form the special editor of calculix called 'SciTE.' We can go to examples directory and just basically right click on disc example and open with by editing. SciTE automatically will be opened.

From below we can easly see the properties of all components of disc. Especially point properties are boxed.

## Generating Different Meshes

We will do the same procedures as we have done before. Our model is prepared. So we will mesh our model and take properties of different kinds of meshes for 2-D. Our meshes for 2-D are:

 $\displaystyle \begin{matrix} tr3 & 3 node & 2D element \\ gu4 & 4 node & 2D element\\ qu8 & 8 node & 2D element \end{matrix}$ (N)
 $\displaystyle \begin{matrix} Mesh & with & tr3 \end{matrix}$
 $\displaystyle \begin{matrix} Mesh & with & qu4 \end{matrix}$
 $\displaystyle \begin{matrix} Mesh & with & qu8 \end{matrix}$

## CCX module

In this module we will analyze for unknowns that we are interested in. This unknowns can be displacements or temperatures. Also we can find fluxes between these nodes.

For the first step go to directory where beam example stays in your harddrive.

 $\displaystyle ....[Program Files]\Rightarrow [bConverged]\Rightarrow [Calculix]\Rightarrow [ccx]\Rightarrow [test]$

Find beampl.inp and right click on it. Then edit it by selecting open with. SciTe automatically will be opened. Every command inside is ready so we have assumed the we have completed by the mesh and ready for solving.So click tools and click 'solve'.

You have to see that $\displaystyle Job Finished$ Then click tools and click 'post prosess'.

Results

 $\displaystyle Result 1$
 $\displaystyle Result2$

# Problem 5.6: The Quadratic Lagrangian Element Basis Function

Lecture 30-6

## Given:

Given the Linear Lagrangian Element Basis Function,

from Lecture 30-6

## Find:

The Quadratic Lagrangian Element Basis Function

## Solution:

The above is the depiction of the basis function satisfying the CBS for element node i-1

The above is the depiction of the basis function satisfying the CBS for element node i

The above is the depiction of the basis function satisfying the CBS for element node i+1

Next we overlay the three functions to show the shape function for the quadratic Lagrangian

# Problem 5.7 Weak form solution using LLEIF

## Given

Solve the following discrete weak form with Linear Lagrange Element Interpolation Basis.

 $\displaystyle \sum {{c_i}\left[ {\sum {{{\tilde K}_{ij}}{d_j} - {{\tilde F}_i}} } \right]} = 0$ (3.7.1)

Where

 $\displaystyle {{\tilde K}_{ij}} = \int\limits_0^1 {{{b'}_i}(2 + 3x){{b'}_j}dx}$ (3.7.2)
 $\displaystyle {F_i} = {b_i}\left( 1 \right)\left( {12} \right) + \int_0^1 {{b_i}\left( x \right)5x} dx$ (3.7.3)

Such that

 $\displaystyle {u^h}\left( 0 \right) = \sum {{d_j}{b_j}\left( 0 \right) = 4}$ (3.7.4)
 $\displaystyle {w^h}\left( 0 \right) = \sum {{c_i}{b_i}\left( 0 \right) = 0}$ (3.7.5)

## Find

### d. Plot ${u^h}$ vs ${u}$

Plot ${u^h}\left( {0.5} \right) - u\left( {0.5} \right)$ vs n (total no. of dof's)

## Solution

### a. Using LLEBF in constraint breaking solution(CBS)

We have following constraint to be satisfied to use weak form

 $\displaystyle {w^h}\left( 0 \right) = \sum {{c_i}{b_i}\left( 0 \right) = 0}$ (3.7.6)

At node 1 (i.e. at x=0) all the basis function are zero except 1.

 $\displaystyle {w^h} = {c_1}{b_1} + {c_2}{b_2} + {c_3}{b_3} + {c_4}{b_4} + \cdots$ (3.7.7)

 $\displaystyle {w^h} = {c_1}(1) + {c_2}(0) + {c_3}(0) + {c_4}(0) + \cdots$ (3.7.8)
 $\displaystyle {w^h} = {c_1}$ (3.7.9)

For

 $\displaystyle {w^h}\left( 0 \right) = 0$ (3.7.10)
 $\displaystyle {c_1} = 0$ (3.7.11)

To satisfy homogeneous boundary condition while using Linear Lagrange Element Interpolation Functions, we just have to select our $w^h$ such that ${c_1} = 0$.

### c Formulation of K and F on element level

We can divide the total domain into elemental domain.

 $\displaystyle \Omega = \bigcup\limits_{e = 1}^{nel} {{w^e}}$ (3.7.12)

Values of K and F will be summation of K and F over the elements

 $\displaystyle {{\tilde K}_{ij}} = \sum\limits_{e = 1}^{nel} {K_{ij}^e}$ (3.7.13)

 $\displaystyle K_{ij}^e = \int\limits_{{w_e}} {{{b'}^e}_i(2 + 3x){{b'}^e}_jdx}$ (3.7.14)

 $\displaystyle {{\tilde F}_i} = \sum\limits_{e = 1}^{nel} {F_i^e}$ (3.7.15)
 $\displaystyle {F_i}^e = \int_{{w_e}} {b_i^e5xdx}$ (3.7.16)

There will be two basis function for each element. Generalised values of these basis function will be as follow

 $\displaystyle {b_e} = \frac{{x - {x_{e + 1}}}}{{{x_e} - {x_{e + 1}}}}$ (3.7.17)
 $\displaystyle {b_{e + 1}} = \frac{{x - {x_e}}}{{{x_{e + 1}} - {x_e}}}$ (3.7.18)

### d. Solution using Matlab

Plot of exact solution and approximate solution with 2 elements

Figure 1: For elements = 2 [8]

Plot of exact solution and approximate solution with 4 elements

Figure 1: For elements = 4 [9]

Plot of exact solution and approximate solution with 6 elements

Figure 1: For elements = 6 [10]

Plot of exact solution and approximate solution with 8 elements

Figure 1: For elements = 8 [11]

Plot of error vs dof's

 $\displaystyle err = u\left( {0.5} \right) - {u^h}\left( {0.5} \right)$ (3.7.19)
 $\displaystyle dof's = nel + 1$ (3.7.20)

where nel= number of elements.

Figure 1: Error vs dof's [12]

for 16 elements

dof's=16+1=17

Figure 1: Error vs dof's [13]

## comment

Solution is converging with increasing no. of elements. For n=16 error is less than 0.0001.

# Problem 5.8:Using LLEBF to solve G1DM1.0/D1

## Given

G1DM1.0/D1 : Model for elastic bar problem was given in meeting 9.

 $\displaystyle \frac{\mathrm{d} }{\mathrm{d} x}\left [ (2+3x)\frac{\mathrm{d} u}{\mathrm{d} x} \right ]+5x=0$ (8.1)
 $\displaystyle u(1)=4$ (8.2)
 $\displaystyle \frac{\mathrm{d} u(x=0)}{\mathrm{d} x}=-3$ (8.3)

Weak Form

 $\displaystyle \int_{\Omega }\frac{\mathrm{d} w}{\mathrm{d} x}AE\frac{\mathrm{d} u}{\mathrm{d} x}dx-\left ( wAt \right )_{\Gamma _{h}}-\int_{\Omega }wbdx=0 \Rightarrow \Omega =]0,1[$ (8.4)

We have to manipulate weak form in order to use for element matrices. As to element mesh we use weak form for each elements and then we can take summation for all elements to get global matrices.

 $\displaystyle \sum\limits_{e=1}^{nel}{{{w}^{eT}}\left\{ \int\limits_{{{x}_{1}}^{e}}^{{{x}_{2}}^{e}}{{{B}^{eT}}{{A}^{e}}{{E}^{e}}{{B}^{e}}dx{{d}^{e}}-\int\limits_{{{x}_{1}}^{e}}^{x_{2}^{e}}{{{N}^{eT}}bdx-{{({{N}^{eT}}{{A}^{e}}\overline{t})}_{x=0}}}} \right\}=0}$ (8.5)

Approximate functions gives as in F&B Book,

 $\displaystyle u^{h}=\sum_{e=1}^{nel}N^{e}d^{e}=\left ( \sum_{e=1}^{nel}N^{e}L^{e}\right )\mathbf{d}$ (8.6)
 $\displaystyle w^{h}=\sum_{e=1}^{nel}N^{e}w^{e}=\left ( \sum_{e=1}^{nel}N^{e}L^{e}\right )\mathbf{w}$ (8.7)

Where N is called as 'shape functions' and L is called 'gather matrices'.

 $\displaystyle N^{e}=\frac{1}{l^{e}}[\begin{matrix} x_{2}^{e}-x & x-x_{1}^{e} \end{matrix}]$ (8.8)
 $\displaystyle B^{e}=\frac{1}{l^{e}}[\begin{matrix} -1 & 1 \end{matrix}]$ (8.9)

## Solution

### Analytical Solution

 $\displaystyle \int \left ( (2+3x)\frac{\mathrm{d} u}{\mathrm{d} x} \right )dx=-\int (5x)dx+c_{1}$ (8.10)
 $\displaystyle 2\frac{\mathrm{d} u}{\mathrm{d} x}(x=0)=0+c_{1} \Rightarrow c_{1}=-6$ (8.11)
 $\displaystyle u=-\frac{5}{2}\int \frac{x^2}{2+3x}dx-\int \frac{6}{2+3x}dx+c_{2}$ (8.12)
 $\displaystyle -\frac{5}{2}\left [ \frac{1}{54}(9x^2-12x+8log(3x+2)-12) \right ]-\left [ \frac{1}{3}log(3x+2) \right ]+c_{2}$ (8.13)

Using other boundry condition ;

 $\displaystyle c_{2}=\frac{119}{36}+\frac{64}{27}log(5)$ (8.14)
 $\displaystyle u(x)=-\frac{5}{12}x^2+\frac{5}{9}x-\frac{64}{27}log(2+3x)+\frac{139}{36}+\frac{64}{27}log(5)$ (8.15)

### FEA solutions

When evaluating unknowns on distributed nodes our shape functions must meet some conditions.Our approximate functions must be complete in each element. As saying completeness we are saying that by using our nodal values we must be able to identify coefficients in our approximate solutions. After some manipulation we conclude with shape functions. For each nodes we have shape functions and for middle nodes shape functions come from adjacent elements.

For example two element mesh ;

 $\displaystyle N=N^{(1)}L^{(1)}+N^{(2)}L^{(2)}=[\begin{matrix} N_{1}^{(1)} & N_{2}^{(1)}+N_{1}^{(2)} & N_{2}^{(2)} \end{matrix}]$ (8.16)

For this example at x1 (first node) only N1 will be 1. Others will be zero. So we verify CBS.[1]

#### Using 4 elements

 $\displaystyle \Omega =]0,1[$

Lets find shape functions for 4 elements.

 \displaystyle\begin{align} &N^{(1)}=\frac{1}{0.25}[\begin{matrix} 0.25-x & x \end{matrix}]\\ &N^{(2)}=\frac{1}{0.25}[\begin{matrix} 0.5-x & x-0.25 \end{matrix}]\\ &N^{(3)}=\frac{1}{0.25}[\begin{matrix} 0.75-x & x-0.5 \end{matrix}]\\ &N^{(4)}=\frac{1}{0.25}[\begin{matrix} 1-x & x-0.75 \end{matrix}] \end{align} (8.17)

These shape functions are useful for us to make connection between nodes. After defining these tools now we can find stiffness matrices for all elements. Then after finding all of stiffness matrices for elements we can construct global stiffness matrix.

 $\displaystyle K^{e}=\int_{x_{1}^{e}}^{x_{2}^{e}}B^{eT}A^{e}E^{e}B^{e}dx$ (8.18)
 \displaystyle\begin{align} {{K}^{e}}&=\int\limits_{x_{1}^{e}}^{x_{2}^{e}}{\frac{1}{{{l}^{e}}}\left[ \begin{matrix} -1 \\ 1 \\ \end{matrix} \right]}(2+3x)\frac{1}{{{l}^{e}}}[\begin{matrix} -1 & 1 \\ \end{matrix}]dx \\ &=\frac{1}{{{({{l}^{e}})}^{2}}}\left[ \begin{matrix} 1 & -1 \\ -1 & 1 \\ \end{matrix} \right]\left[ 2(x_{2}^{e}-x_{1}^{e})+\frac{3}{2}(x_{2}^{e2}-x_{1}^{e2}) \right] \end{align} (8.19)

 $\displaystyle {{K}^{e}}=\frac{1}{{{l}^{e}}}\left[ \begin{matrix} 1 & -1 \\ -1 & 1 \\ \end{matrix} \right]\left[ 2+\frac{3}{2}(x_{2}^{e}+x_{1}^{e}) \right]$ (8.20)
 $\displaystyle K^{(1)}=\begin{bmatrix} 9.5 & -9.5\\ -9.5& 9.5 \end{bmatrix}$ (8.21)
 $\displaystyle K^{(2)}=\begin{bmatrix} 12.5 & -12.5\\ -12.5& 12.5 \end{bmatrix}$ (8.22)
 $\displaystyle K^{(3)}=\begin{bmatrix} 15.5 & -15.5\\ -15.5& 15.5 \end{bmatrix}$ (8.23)
 $\displaystyle K^{(4)}=\begin{bmatrix} 18.5 & -18.5\\ -18.5& 18.5 \end{bmatrix}$ (8.24)

We can construct global stiffness matrix now. This will be distribution of all four element stiffness matrices into appropriate places. For our convenience i will use direct method. But it is useful using below equation for computing purposes.

 $\displaystyle K=\sum_{e=1}^{nel}L^{eT}K^{e}L^{e}=L^{1T}K^{1}L^{1}+L^{2T}K^{2}L^{2}+L^{3T}K^{3}L^{3}+L^{4T}K^{4}L^{4}$ (8.25)
 $\displaystyle K=\left[ \begin{matrix} \begin{matrix} \begin{matrix} \begin{matrix} 9.5 \\ -9.5 \\ \end{matrix} \\ 0 \\ 0 \\ 0 \\ \end{matrix} & \begin{matrix} \begin{matrix} -9.5 \\ 22 \\ \end{matrix} \\ -12.5 \\ 0 \\ 0 \\ \end{matrix} \\ \end{matrix} & \begin{matrix} \begin{matrix} 0 \\ -12.5 \\ \end{matrix} \\ 28 \\ \begin{matrix} -15.5 \\ 0 \\ \end{matrix} \\ \end{matrix} & \begin{matrix} \begin{matrix} 0 \\ 0 \\ \end{matrix} \\ -15.5 \\ 34 \\ -18.5 \\ \end{matrix} & \begin{matrix} \begin{matrix} 0 \\ 0 \\ \end{matrix} \\ 0 \\ -18.5 \\ 18.5 \\ \end{matrix} \\ \end{matrix} \right]$ (8.26)

The next step we have to find force matrices. Which ha two components. Lets first find force matrix at natural boundary condition and second distributed force matrix on our domain.

 $\displaystyle f_{\Gamma _{h}}=-(N^{eT}A^{e}b)_{\Gamma _{h}^{e}}=-N^{eT}(x_{1})(2+3x)_{x=0}(-3)=6N^{eT}(x_{1})$ (8.27)
 $\displaystyle f_{{\Gamma }^{(1)}}=6\left[ \begin{matrix} N_{1}^{(1)}({{x}_{1}}) \\ N_{2}^{(1)}({{x}_{2}}) \\ \end{matrix} \right]=\left[ \begin{matrix} 6 \\ 0 \\ \end{matrix} \right]$ (8.28)
 \displaystyle\begin{align} &f_{{\Gamma }^{(2)}}=f_{{\Gamma }^{(3)}}=f_{{\Gamma }^{(4)}}=6\left[ \begin{matrix} N_{1}^{(e)}({{x}_{1}}) \\ N_{2}^{(e)}({{x}_{2}}) \\ \end{matrix} \right]=\left[ \begin{matrix} 0 \\ 0 \\ \end{matrix} \right] &for & e=2,3,4 \end{align} (8.29)
 \displaystyle\begin{align} f_{\Gamma _{h}} &=\sum_{e=1}^{nel}L^{eT}f_{\Gamma _{h}}^{e} \\ &=\left[ \begin{matrix} \begin{matrix} 6 \\ 0 \\ \end{matrix} \\ 0 \\ 0 \\ 0 \\ \end{matrix} \right] \end{align} (N)
 \displaystyle\begin{align} f_{\Omega }^{e} & =\int_{x_{1}^{e}}^{x_{2}^{e}}N^{eT}b(x)dx\\ &=\int\limits_{x_{1}^{e}}^{x_{2}^{e}}{\frac{1}{{{l}^{e}}}\left[ \begin{matrix} x_{2}^{e}-x \\ x-x_{1}^{e} \\ \end{matrix} \right]}5xdx \\ &=\frac{5}{{{l}^{e}}}\left[ \begin{matrix} \frac{x_{2}^{e}(x_{2}^{e2}-x_{1}^{e2})}{2}-\frac{(x_{2}^{e3}-x_{1}^{e3})}{3} \\ \frac{(x_{2}^{e3}-x_{1}^{e3})}{3}-\frac{x_{1}^{e}(x_{2}^{e2}-x_{1}^{e2})}{2} \\ \end{matrix} \right] \end{align} (8.30)
 $\displaystyle f_{\Omega }^{(1)}=\left[ \begin{matrix} 0.05208334 \\ 0.10416666 \\ \end{matrix} \right]$ (8.31)
 $\displaystyle f_{\Omega }^{(2)}=\left[ \begin{matrix} 0.2083334 \\ 0.2604166 \\ \end{matrix} \right]$ (8.32)
 $\displaystyle f_{\Omega }^{(3)}=\left[ \begin{matrix} 0.36458334 \\ 0.41666666 \\ \end{matrix} \right]$ (8.33)
 $\displaystyle f_{\Omega }^{(4)}=\left[ \begin{matrix} 0.52083334 \\ 0.57291666 \\ \end{matrix} \right]$ (8.34)
 \displaystyle\begin{align} f_{\Omega }&=\sum_{e=1}^{4}L^{eT}f_{\Omega }^{e} &\left[ \begin{matrix} \begin{matrix} 0.05208334 \\ 0.3125 \\ \end{matrix} \\ 0.625 \\ 0.9375 \\ 0.57291666 \\ \end{matrix} \right] \end{align} (8.35)
 $\displaystyle \left[ \begin{matrix} \begin{matrix} 9.5 \\ -9.5 \\ \end{matrix} & \begin{matrix} -9.5 \\ 22 \\ \end{matrix} & \begin{matrix} 0 \\ -12.5 \\ \end{matrix} & \begin{matrix} \begin{matrix} 0 & 0 \\ \end{matrix} \\ \begin{matrix} 0 & 0 \\ \end{matrix} \\ \end{matrix} \\ 0 & -12.5 & 28 & \begin{matrix} -15.5 & 0 \\ \end{matrix} \\ 0 & 0 & -15.5 & \begin{matrix} 34 & -18.5 \\ \end{matrix} \\ 0 & 0 & 0 & \begin{matrix} -18.5 & 18.5 \\ \end{matrix} \\ \end{matrix} \right]\left[ \begin{matrix} \begin{matrix} {{d}_{1}} \\ {{d}_{2}} \\ \end{matrix} \\ {{d}_{3}} \\ {{d}_{4}} \\ 4 \\ \end{matrix} \right]=\left[ \begin{matrix} \begin{matrix} 0.05208334 \\ 0.3125 \\ \end{matrix} \\ 0.625 \\ 0.9375 \\ 0.57291666 \\ \end{matrix} \right]+\left[ \begin{matrix} \begin{matrix} 6 \\ 0 \\ \end{matrix} \\ 0 \\ 0 \\ 0 \\ \end{matrix} \right]+\left[ \begin{matrix} \begin{matrix} 0 \\ 0 \\ \end{matrix} \\ 0 \\ 0 \\ {{r}_{1}} \\ \end{matrix} \right]$ (8.36)
 $\displaystyle \left[ \begin{matrix} {{d}_{1}} \\ {{d}_{2}} \\ {{d}_{3}} \\ {{d}_{4}} \\ \end{matrix} \right]=\left[ \begin{matrix} 2.0257 \\ 1.3886 \\ 0.8794 \\ 0.4285 \\ \end{matrix} \right]$ (8.37)

Comparing results with exact ones at the nodes;

 \displaystyle\begin{align} &\left[ \begin{matrix} 2.0257 \\ 1.3886 \\ 0.8794 \\ 0.4285 \\ \end{matrix} \right]+\left[ \begin{matrix} 4 \\ 4 \\ 4 \\ 4 \\ \end{matrix} \right]=\left[ \begin{matrix} 6.0257 \\ 5.3886 \\ 4.8794 \\ 4.4285 \\ \end{matrix} \right] &{{d}_{exact}}=\left[ \begin{matrix} 6.0331 \\ 5.3911 \\ 4.8802 \\ 4.4286 \\ \end{matrix} \right] \end{align} (8.38)
 $\displaystyle 4 elements$

#### Using 6 elements

Shape functions for 6 elements

 \displaystyle\begin{align} & N^{(1)}=\frac{1}{0.1666}[\begin{matrix} 0.1666-x & x \end{matrix}]\\ &N^{(2)}=\frac{1}{0.1666}[\begin{matrix} 0.3333-x & x-0.16666 \end{matrix}]\\ & N^{(3)}=\frac{1}{0.1666}[\begin{matrix} 0.5-x & x-0.3333 \end{matrix}]\\ &N^{(4)}=\frac{1}{0.1666}[\begin{matrix} 0.6666-x & x-0.5 \end{matrix}]\\ & N^{(5)}=\frac{1}{0.1666}[\begin{matrix} 0.8332-x & x-0.6666 \end{matrix}]\\ & N^{(6)}=\frac{1}{0.1666}[\begin{matrix} 1-x & x-0.8332 \end{matrix}]\\ \end{align} (8.39)
 $\displaystyle \left[ \begin{matrix} \begin{matrix} \begin{matrix} \begin{matrix} 13.5048 & -13.5048 & 0 & 0 \\ \end{matrix} & 0 & 0 & 0 \\ \end{matrix} \\ \begin{matrix} \begin{matrix} -13.5048 & 30.0105 & -16.5057 & 0 \\ \end{matrix} & 0 & 0 & 0 \\ \end{matrix} \\ \begin{matrix} \begin{matrix} 0 & -16.5057 & 36.0132 & -19.5075 \\ \end{matrix} & 0 & 0 & 0 \\ \end{matrix} \\ \begin{matrix} \begin{matrix} 0 & 0 & -19.5075 & 42.0157 \\ \end{matrix} & -22.5084 & 0 & 0 \\ \end{matrix} \\ \end{matrix} \\ \begin{matrix} \begin{matrix} 0 & 0 & 0 & -22.5084 \\ \end{matrix} & 48.0164 & -25.508 & 0 \\ \end{matrix} \\ \begin{matrix} \begin{matrix} 0 & 0 & 0 & 0 \\ \end{matrix} & -25.508 & 54.008 & -28.5 \\ \end{matrix} \\ \begin{matrix} \begin{matrix} 0 & 0 & 0 & 0 \\ \end{matrix} & 0 & -28.5 & 28.5 \\ \end{matrix} \\ \end{matrix} \right]\left[ \begin{matrix} \begin{matrix} {{d}_{1}} \\ {{d}_{2}} \\ {{d}_{3}} \\ {{d}_{4}} \\ \end{matrix} \\ {{d}_{5}} \\ {{d}_{6}} \\ 4 \\ \end{matrix} \right]=\left[ \begin{matrix} \begin{matrix} 6.0232199942 \\ 0.13890282 \\ 0.277972087 \\ 0.416641681 \\ \end{matrix} \\ 0.555277743 \\ 0.694970317 \\ 0.394287575 \\ \end{matrix} \right]+\left[ \begin{matrix} \begin{matrix} 0 \\ 0 \\ 0 \\ 0 \\ \end{matrix} \\ 0 \\ 0 \\ {{r}_{1}} \\ \end{matrix} \right]$ (8.40)
 $\displaystyle \left[ \begin{matrix} \begin{matrix} {{d}_{1}} \\ {{d}_{2}} \\ {{d}_{3}} \\ \end{matrix} \\ {{d}_{4}} \\ {{d}_{5}} \\ {{d}_{6}} \\ \end{matrix} \right]=\left[ \begin{matrix} \begin{matrix} 2.0292 \\ 1.5831 \\ 1.2098 \\ \end{matrix} \\ 0.8797 \\ 0.575 \\ 0.2845 \\ \end{matrix} \right]$ (8.41)
 \displaystyle\begin{align} & \left[ \begin{matrix} \begin{matrix} 2.0292 \\ 1.5831 \\ 1.2098 \\ \end{matrix} \\ 0.8797 \\ 0.575 \\ 0.2845 \\ \end{matrix} \right]+\left[ \begin{matrix} \begin{matrix} 4 \\ 4 \\ 4 \\ \end{matrix} \\ 4 \\ 4 \\ 4 \\ \end{matrix} \right]=\left[ \begin{matrix} \begin{matrix} 6.0292 \\ 5.5831 \\ 5.2098 \\ \end{matrix} \\ 4.8797 \\ 4.575 \\ 4.2845 \\ \end{matrix} \right] & {{d}_{exact}}=\left[ \begin{matrix} \begin{matrix} 6.0331 \\ 5.5853 \\ 5.2109 \\ \end{matrix} \\ 4.8802 \\ 4.5753 \\ 4.2847 \\ \end{matrix} \right] \end{align} (8.42)
 $\displaystyle 6 elements$

#### Using 8 elements

Shape functions for 8 elements

 \displaystyle\begin{align} & N^{(1)}=\frac{1}{0.125}[\begin{matrix} 0.125-x & x \end{matrix}]\\ &N^{(2)}=\frac{1}{0.125}[\begin{matrix} 0.25-x & x-0.125 \end{matrix}]\\ & N^{(3)}=\frac{1}{0.125}[\begin{matrix} 0.375-x & x-0.25 \end{matrix}]\\ &N^{(4)}=\frac{1}{0.125}[\begin{matrix} 0.5-x & x-0.375 \end{matrix}]\\ & N^{(5)}=\frac{1}{0.125}[\begin{matrix} 0.625-x & x-0.5 \end{matrix}]\\ & N^{(6)}=\frac{1}{0.125}[\begin{matrix} 0.75-x & x-0.625 \end{matrix}]\\ & N^{(7)}=\frac{1}{0.125}[\begin{matrix} 0.875-x & x-0.75 \end{matrix}]\\ & N^{(8)}=\frac{1}{0.125}[\begin{matrix} 1-x & x-0.875 \end{matrix}]\\ \end{align} (8.43)
 $\displaystyle \left[ \begin{matrix} 17.5 & -17.5 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ -17.5 & 38 & -20.5 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & -20.5 & 44 & -23.5 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & -23.5 & 50 & -26.5 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & -26.5 & 56 & -29.5 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & -29.5 & 62 & -32.5 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & -32.5 & 68 & -35.5 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & -35.5 & 74 & -38.5 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & -38.5 & 38.5 \\ \end{matrix} \right]\left[ \begin{matrix} \begin{matrix} \begin{matrix} {{d}_{1}} \\ {{d}_{2}} \\ {{d}_{3}} \\ \end{matrix} \\ {{d}_{4}} \\ {{d}_{5}} \\ {{d}_{6}} \\ \end{matrix} \\ {{d}_{7}} \\ {{d}_{8}} \\ {{d}_{9}} \\ \end{matrix} \right]=\left[ \begin{matrix} \begin{matrix} \begin{matrix} \begin{matrix} 6.01302084 \\ 0.078125 \\ \end{matrix} \\ 0.15625 \\ \end{matrix} \\ 0.234375 \\ 0.31250004 \\ 0.390625 \\ \end{matrix} \\ 0.46874373 \\ 0.546875 \\ 0.29947916 \\ \end{matrix} \right]+\left[ \begin{matrix} \begin{matrix} \begin{matrix} 0 \\ 0 \\ 0 \\ \end{matrix} \\ 0 \\ 0 \\ 0 \\ \end{matrix} \\ 0 \\ 0 \\ {{r}_{1}} \\ \end{matrix} \right]$ (8.44)
 $\displaystyle \left[ \begin{matrix} \begin{matrix} \begin{matrix} \begin{matrix} {{d}_{1}} \\ {{d}_{2}} \\ \end{matrix} \\ {{d}_{3}} \\ \end{matrix} \\ {{d}_{4}} \\ {{d}_{5}} \\ \end{matrix} \\ {{d}_{6}} \\ {{d}_{7}} \\ {{d}_{8}} \\ \end{matrix} \right]=\left[ \begin{matrix} \begin{matrix} \begin{matrix} 2.0312 \\ 1.6876 \\ \end{matrix} \\ 1.3904 \\ 1.1246 \\ 0.8800 \\ \end{matrix} \\ 0.6497 \\ 0.4286 \\ 0.213 \\ \end{matrix} \right]$ (8.45)
 \displaystyle\begin{align} &\left[ \begin{matrix} \begin{matrix} \begin{matrix} 2.0312 \\ 1.6876 \\ \end{matrix} \\ 1.3904 \\ 1.1246 \\ 0.8800 \\ \end{matrix} \\ 0.6497 \\ 0.4286 \\ 0.213 \\ \end{matrix} \right]+\left[ \begin{matrix} \begin{matrix} \begin{matrix} 4 \\ 4 \\ \end{matrix} \\ 4 \\ 4 \\ 4 \\ \end{matrix} \\ 4 \\ 4 \\ 4 \\ \end{matrix} \right]=\left[ \begin{matrix} \begin{matrix} \begin{matrix} 6.0312 \\ 5.6876 \\ \end{matrix} \\ 5.3904 \\ 5.1246 \\ 4.88 \\ \end{matrix} \\ 4.6497 \\ 4.4286 \\ 4.213 \\ \end{matrix} \right] &{{d}_{exact}}=\left[ \begin{matrix} \begin{matrix} \begin{matrix} 6.0331 \\ 5.6866 \\ \end{matrix} \\ 5.3911 \\ 5.1249 \\ 4.8802 \\ \end{matrix} \\ 4.6498 \\ 4.4286 \\ 4.213 \\ \end{matrix} \right] \end{align} (8.46)
 $\displaystyle 8 elements$

# References

1. Fish and Belytschko FEA Anlaysis

# Contributing Members & Referenced Lecture

 Team Contribution Table Problem Number Lecture Assigned To Solved By Typed By Proofread By 5.1 Lecture 26-2 Chris,Tim Chris,Tim Chris,Tim XXX 5.2 Lecture 26-2 Chris,Tim Chris,Tim Chris,Tim XXX 5.3 Lecture 29-6 Danny,Phil XXX XXX XXX 5.4 Lecture 29-6 Danny,Phil XXX XXX XXX 5.5 Lecture 29-7 Chris,Tim Oztekin Oztekin XXX 5.6 Lecture 30-6 Danny,Phil XXX XXX xxx 5.7 Lecture 30-6 Erman,Taj Shahtaj Shahtaj Oztekin 5.8 Lecture 30-7 Erman,Taj Oztekin Oztekin XXX


text:
### Mathematics Class IV

Numbers
The Four Operations
Multiples and Factors
HCF and LCM
Fractions
Decimals
Metric Measures
Time
Money, Profit And Loss
Geometry
Symmetry and Pattern
Perimeter and Area
Data Handling

# Writing a Decimal Fraction as Common Fraction

Lesson Progress
0% Complete

Step 1:
Write the number given in the decimal fraction without decimal point as numerator.

Step 2:
Write the denominator with .

Step 3:
Add number of zeroes equal to the number of digits after the decimal point on the right side of .

Example 1.

Write as a fraction.
Write the number given in the decimal fraction without decimal point as numerator.

Write the denominator with .

Add number of zeroes equal to the number of digits after the decimal point on the right side of .

Example 2.

Write as a fraction.
Write the number given in the decimal fraction without decimal point as numerator.

Write the denominator with .

Add number of zeroes equal to the number of digits after the decimal point on the right side of .

Scroll to Top

format: plain


text:
HOSTING A TOTAL OF 318 FORMULAS WITH CALCULATORS

## Radiation pressure (absorption)

Radiation pressure is the pressure exerted upon any surface exposed to electromagnetic radiation. The term can also refer to pressures generated by other forms of wave motion, e.g. acoustic radiation pressure. According to special relativity, the momentum of a photon is equal to its energy divided by the speed of light.

if absorbed, the pressure is the power flux density divided by the speed of light

## $f\frac{a}{c}\mathrm{COS\Theta }$

Here,c=speed of light,a=power density,f=force.

ENTER THE VARIABLES TO BE USED IN THE FORMULA

Similar formulas which you may find interesting.


text:
# Queries¶

## Outline¶

1. First we’ll make a connection to the Gaia server,

2. We will explore information about the database and the tables it contains,

3. We will write a query and send it to the server, and finally

## Query Language¶

In order to select data from a database, you have to compose a query, which is a program written in a “query language”. The query language we’ll use is ADQL, which stands for “Astronomical Data Query Language”.

ADQL is a dialect of SQL (Structured Query Language), which is by far the most commonly used query language. Almost everything you will learn about ADQL also works in SQL.

The reference manual for ADQL is here. But you might find it easier to learn from this ADQL Cookbook.

## Using Jupyter¶

If you have not worked with Jupyter notebooks before, you might start with the tutorial on from Jupyter.org called “Try Classic Notebook”, or this tutorial from DataQuest.

There are two environments you can use to write and run notebooks:

• “Jupyter Notebook” is the original, and

• “Jupyter Lab” is a newer environment with more features.

For these lessons, you can use either one.

If you are too impatient for the tutorials, here’s are the most important things to know:

1. Notebooks are made up of code cells and text cells (and a few other less common kinds). Code cells contain code; text cells, like this one, contain explanatory text written in Markdown.

2. To run a code cell, click the cell to select it and press Shift-Enter. The output of the code should appear below the cell.

3. In general, notebooks only run correctly if you run every code cell in order from top to bottom. If you run cells out of order, you are likely to get errors.

4. You can modify existing cells, but then you have to run them again to see the effect.

5. You can add new cells, but again, you might have to be careful about the order you run them in.

6. If you have added or modified cells and the behavior of the notebook seems strange, you can restart the “kernel”, which clears all of the variables and functions you have defined, and run the cells again from the beginning.

• If you are using Jupyter notebook, open the Kernel menu and select “Restart and Run All”.

• In Jupyter Lab…

• In Colab, open the Runtime menu and select “Restart and run all”

Before you go on, you might want to explore the other menus and the toolbar to see what else you can do.

## Connecting to Gaia¶

The library we’ll use to get Gaia data is Astroquery. Astroquery provides Gaia, which is an object that represents a connection to the Gaia database.

We can connect to the Gaia database like this:

from astroquery.gaia import Gaia

Created TAP+ (v1.2.1) - Connection:
Host: gea.esac.esa.int
Use HTTPS: True
Port: 443
SSL Port: 443
Created TAP+ (v1.2.1) - Connection:
Use HTTPS: True
Port: 443
SSL Port: 443

Running this import statement has the effect of creating a TAP+ connection; TAP stands for “Table Access Protocol”, which is a network protocol for sending queries to the database and getting back the results.

## Databases and Tables¶

What is a database, anyway? Most generally, it can be any collection of data, but when we are talking about ADQL or SQL:

• A database is a collection of one or more named tables.

• Each table is a 2-D array with one or more named columns of data.

We can use Gaia.load_tables to get the names of the tables in the Gaia database. With the option only_names=True, it loads information about the tables, called “metadata”, not the data itself.

tables = Gaia.load_tables(only_names=True)

INFO: Retrieving tables... [astroquery.utils.tap.core]
INFO: Parsing tables... [astroquery.utils.tap.core]
INFO: Done. [astroquery.utils.tap.core]

for table in tables:
print(table.name)

external.apassdr9
external.galex_ais
external.ravedr5_com
external.ravedr5_dr5
external.ravedr5_gra
external.ravedr5_on
external.sdssdr13_photoprimary
external.skymapperdr1_master
external.skymapperdr2_master
external.tmass_xsc
public.hipparcos
public.hipparcos_newreduction
public.hubble_sc
public.igsl_source
public.igsl_source_catalog_ids
public.tycho2
public.dual
tap_config.coord_sys
tap_config.properties
tap_schema.columns
tap_schema.key_columns
tap_schema.keys
tap_schema.schemas
tap_schema.tables
gaiaedr3.gaia_source
gaiaedr3.agn_cross_id
gaiaedr3.commanded_scan_law
gaiaedr3.dr2_neighbourhood
gaiaedr3.frame_rotator_source
gaiaedr3.hipparcos2_best_neighbour
gaiaedr3.hipparcos2_neighbourhood
gaiaedr3.panstarrs1_best_neighbour
gaiaedr3.panstarrs1_join
gaiaedr3.panstarrs1_neighbourhood
gaiaedr3.sdssdr13_best_neighbour
gaiaedr3.sdssdr13_join
gaiaedr3.sdssdr13_neighbourhood
gaiaedr3.skymapperdr2_best_neighbour
gaiaedr3.skymapperdr2_join
gaiaedr3.skymapperdr2_neighbourhood
gaiaedr3.tycho2tdsc_merge_best_neighbour
gaiaedr3.tycho2tdsc_merge_neighbourhood
gaiaedr3.urat1_best_neighbour
gaiaedr3.urat1_neighbourhood
gaiaedr3.gaia_source_simulation
gaiaedr3.gaia_universe_model
gaiaedr3.tycho2tdsc_merge

So that’s a lot of tables. The ones we’ll use are:

• gaiadr2.gaia_source, which contains Gaia data from data release 2,

• gaiadr2.panstarrs1_original_valid, which contains the photometry data we’ll use from PanSTARRS, and

• gaiadr2.panstarrs1_best_neighbour, which we’ll use to cross-match each star observed by Gaia with the same star observed by PanSTARRS.

We can use load_table (not load_tables) to get the metadata for a single table. The name of this function is misleading, because it only downloads metadata, not the contents of the table.

meta = Gaia.load_table('gaiadr2.gaia_source')
meta

Retrieving table 'gaiadr2.gaia_source'
Done.

<astroquery.utils.tap.model.taptable.TapTableMeta at 0x7f2e23f089d0>

Jupyter shows that the result is an object of type TapTableMeta, but it does not display the contents.

To see the metadata, we have to print the object.

print(meta)

TAP Table name: gaiadr2.gaiadr2.gaia_source
Description: This table has an entry for every Gaia observed source as listed in the
Main Database accumulating catalogue version from which the catalogue
release has been generated. It contains the basic source parameters,
that is only final data (no epoch data) and no spectra (neither final
nor epoch).
Num. columns: 96

## Columns¶

The following loop prints the names of the columns in the table.

for column in meta.columns:
print(column.name)

solution_id
designation
source_id
random_index
ref_epoch
ra
ra_error
dec
dec_error
parallax
parallax_error
parallax_over_error
pmra
pmra_error
pmdec
pmdec_error
ra_dec_corr
ra_parallax_corr
ra_pmra_corr
ra_pmdec_corr
dec_parallax_corr
dec_pmra_corr
dec_pmdec_corr
parallax_pmra_corr
parallax_pmdec_corr
pmra_pmdec_corr
astrometric_n_obs_al
astrometric_n_obs_ac
astrometric_n_good_obs_al
astrometric_gof_al
astrometric_chi2_al
astrometric_excess_noise
astrometric_excess_noise_sig
astrometric_params_solved
astrometric_primary_flag
astrometric_weight_al
astrometric_pseudo_colour
astrometric_pseudo_colour_error
mean_varpi_factor_al
astrometric_matched_observations
visibility_periods_used
astrometric_sigma5d_max
frame_rotator_object_type
matched_observations
duplicated_source
phot_g_n_obs
phot_g_mean_flux
phot_g_mean_flux_error
phot_g_mean_flux_over_error
phot_g_mean_mag
phot_bp_n_obs
phot_bp_mean_flux
phot_bp_mean_flux_error
phot_bp_mean_flux_over_error
phot_bp_mean_mag
phot_rp_n_obs
phot_rp_mean_flux
phot_rp_mean_flux_error
phot_rp_mean_flux_over_error
phot_rp_mean_mag
phot_bp_rp_excess_factor
phot_proc_mode
bp_rp
bp_g
g_rp
rv_nb_transits
rv_template_teff
rv_template_logg
rv_template_fe_h
phot_variable_flag
l
b
ecl_lon
ecl_lat
priam_flags
teff_val
teff_percentile_lower
teff_percentile_upper
a_g_val
a_g_percentile_lower
a_g_percentile_upper
e_bp_min_rp_val
e_bp_min_rp_percentile_lower
e_bp_min_rp_percentile_upper
flame_flags
lum_val
lum_percentile_lower
lum_percentile_upper
epoch_photometry_url

You can probably guess what many of these columns are by looking at the names, but you should resist the temptation to guess. To find out what the columns mean, read the documentation.

If you want to know what can go wrong when you don’t read the documentation, you might like this article.

### Exercise¶

One of the other tables we’ll use is gaiadr2.panstarrs1_original_valid. Use load_table to get the metadata for this table. How many columns are there and what are their names?

# Solution

print(meta2)

for column in meta2.columns:
print(column.name)

Retrieving table 'gaiadr2.panstarrs1_original_valid'
Done.
Description: The Panoramic Survey Telescope and Rapid Response System (Pan-STARRS) is
a system for wide-field astronomical imaging developed and operated by
the Institute for Astronomy at the University of Hawaii. Pan-STARRS1
(PS1) is the first part of Pan-STARRS to be completed and is the basis
for Data Release 1 (DR1). The PS1 survey used a 1.8 meter telescope and
its 1.4 Gigapixel camera to image the sky in five broadband filters (g,
r, i, z, y).

The current table contains a filtered subsample of the 10 723 304 629
entries listed in the original ObjectThin table.
We used only ObjectThin and MeanObject tables to extract
panstarrs1OriginalValid table, this means that objects detected only in
stack images are not included here. The main reason for us to avoid the
use of objects detected in stack images is that their astrometry is not
as good as the mean objects astrometry: “The stack positions (raStack,
decStack) have considerably larger systematic astrometric errors than
the mean epoch positions (raMean, decMean).” The astrometry for the
MeanObject positions uses Gaia DR1 as a reference catalog, while the
stack positions use 2MASS as a reference catalog.

In details, we filtered out all objects where:

-   nDetections = 1

-   no good quality data in Pan-STARRS, objInfoFlag 33554432 not set

-   mean astrometry could not be measured, objInfoFlag 524288 set

-   stack position used for mean astrometry, objInfoFlag 1048576 set

-   error on all magnitudes equal to 0 or to -999;

-   all magnitudes set to -999;

-   error on RA or DEC greater than 1 arcsec.

The number of objects in panstarrs1OriginalValid is 2 264 263 282.

The panstarrs1OriginalValid table contains only a subset of the columns
available in the combined ObjectThin and MeanObject tables. A
description of the original ObjectThin and MeanObjects tables can be
found at:
https://outerspace.stsci.edu/display/PANSTARRS/PS1+Database+object+and+detection+tables

http://mastweb.stsci.edu/ps1casjobs/home.aspx
Documentation:
https://outerspace.stsci.edu/display/PANSTARRS
http://pswww.ifa.hawaii.edu/pswww/
References:
The Pan-STARRS1 Surveys, Chambers, K.C., et al. 2016, arXiv:1612.05560
Pan-STARRS Data Processing System, Magnier, E. A., et al. 2016,
arXiv:1612.05240
Pan-STARRS Pixel Processing: Detrending, Warping, Stacking, Waters, C.
Z., et al. 2016, arXiv:1612.05245
Pan-STARRS Pixel Analysis: Source Detection and Characterization,
Magnier, E. A., et al. 2016, arXiv:1612.05244
Pan-STARRS Photometric and Astrometric Calibration, Magnier, E. A., et
al. 2016, arXiv:1612.05242
The Pan-STARRS1 Database and Data Products, Flewelling, H. A., et al.
2016, arXiv:1612.05243

Catalogue curator:
SSDC - ASI Space Science Data Center
https://www.ssdc.asi.it/
Num. columns: 26
obj_name
obj_id
ra
dec
ra_error
dec_error
epoch_mean
g_mean_psf_mag
g_mean_psf_mag_error
g_flags
r_mean_psf_mag
r_mean_psf_mag_error
r_flags
i_mean_psf_mag
i_mean_psf_mag_error
i_flags
z_mean_psf_mag
z_mean_psf_mag_error
z_flags
y_mean_psf_mag
y_mean_psf_mag_error
y_flags
n_detections
zone_id
obj_info_flag
quality_flag

## Writing queries¶

By now you might be wondering how we download the actual data. With tables this big, you generally don’t. Instead, you use queries to select only the data you want.

A query is a string written in a query language like SQL; for the Gaia database, the query language is a dialect of SQL called ADQL.

Here’s an example of an ADQL query.

query1 = """SELECT
TOP 10
source_id, ref_epoch, ra, dec, parallax

Python note: We use a triple-quoted string here so we can include line breaks in the query, which makes it easier to read.

The words in uppercase are ADQL keywords:

• SELECT indicates that we are selecting data (as opposed to adding or modifying data).

• TOP indicates that we only want the first 10 rows of the table, which is useful for testing a query before asking for all of the data.

• FROM specifies which table we want data from.

The third line is a list of column names, indicating which columns we want.

In this example, the keywords are capitalized and the column names are lowercase. This is a common style, but it is not required. ADQL and SQL are not case-sensitive.

To run this query, we use the Gaia object, which represents our connection to the Gaia database, and invoke launch_job:

job1 = Gaia.launch_job(query1)
job1

<astroquery.utils.tap.model.job.Job at 0x7f2e23f2afa0>

The result is an object that represents the job running on a Gaia server.

If you print it, it displays metadata for the forthcoming table.

print(job1)

--------- ------- ---- ------------------------------------------------------------------ -----
source_id   int64      Unique source identifier (unique within a particular Data Release)     0
ref_epoch float64   yr                                                    Reference epoch     0
ra float64  deg                                                    Right ascension     0
dec float64  deg                                                        Declination     0
parallax float64  mas                                                           Parallax     2
Jobid: None
Phase: COMPLETED
Owner: None
Output file: sync_20201229114647.xml.gz
Results: None

Don’t worry about Results: None. That does not actually mean there are no results.

However, Phase: COMPLETED indicates that the job is complete, so we can get the results like this:

results1 = job1.get_results()
type(results1)

astropy.table.table.Table

Optional detail: Why is table repeated three times? The first is the name of the module, the second is the name of the submodule, and the third is the name of the class. Most of the time we only care about the last one. It’s like the Linnean name for gorilla, which is Gorilla gorilla gorilla.

The result is an Astropy Table, which is similar to a table in an SQL database except:

• SQL databases are stored on disk drives, so they are persistent; that is, they “survive” even if you turn off the computer. An Astropy Table is stored in memory; it disappears when you turn off the computer (or shut down this Jupyter notebook).

• SQL databases are designed to process queries. An Astropy Table can perform some query-like operations, like selecting columns and rows. But these operations use Python syntax, not SQL.

Jupyter knows how to display the contents of a Table.

results1

Table length=10
yrdegdegmas
int64float64float64float64float64
67585097575941414402015.5290.4899010727352-30.343172184207830.48023816159705535
67585086924379769602015.5290.64176055082714-30.289720131422252.2625971293368154
67585270362508492802015.5290.35703708993327-30.193657826181596-0.2763960334229464
67585644582982429442015.5290.77376294142846-29.7653684392252380.5907906528352993
67585586128421557762015.5290.67550630386245-29.921133960781690.2858563565989917
67585563794599219202015.5290.5454006212404-29.900709054816964-1.0012355835832834
67585375202603851522015.5290.7341683169994-30.158181298418626--
67585063559800106242015.5290.586950869878-30.371376421675690.3769870991981157
67585616064339681282015.5290.82697136098506-29.75247697212053--
67585641834144084482015.5290.7133096669958-29.7817436736796070.9376387942856869

Each column has a name, units, and a data type.

For example, the units of ra and dec are degrees, and their data type is float64, which is a 64-bit floating-point number, used to store measurements with a fraction part.

This information comes from the Gaia database, and has been stored in the Astropy Table by Astroquery.

### Exercise¶

Read the documentation of this table and choose a column that looks interesting to you. Add the column name to the query and run it again. What are the units of the column you selected? What is its data type?

# Solution

## Asynchronous queries¶

launch_job asks the server to run the job “synchronously”, which normally means it runs immediately. But synchronous jobs are limited to 2000 rows. For queries that return more rows, you should run “asynchronously”, which mean they might take longer to get started.

If you are not sure how many rows a query will return, you can use the SQL command COUNT to find out how many rows are in the result without actually returning them. We’ll see an example of this later.

The results of an asynchronous query are stored in a file on the server, so you can start a query and come back later to get the results.

For anonymous users, files are kept for three days.

As an example, let’s try a query that’s similar to query1, with two changes:

• It selects the first 3000 rows, so it is bigger than we should run synchronously.

• It uses a new keyword, WHERE.

query2 = """SELECT TOP 3000
source_id, ref_epoch, ra, dec, parallax
WHERE parallax < 1
"""

A WHERE clause indicates which rows we want; in this case, the query selects only rows “where” parallax is less than 1. This has the effect of selecting stars with relatively low parallax, which are farther away. We’ll use this clause to exclude nearby stars that are unlikely to be part of GD-1.

WHERE is one of the most common clauses in ADQL/SQL, and one of the most useful, because it allows us to select only the rows we need from the database.

We use launch_job_async to submit an asynchronous query.

job2 = Gaia.launch_job_async(query2)
print(job2)

INFO: Query finished. [astroquery.utils.tap.core]

name    dtype  unit                            description
--------- ------- ---- ------------------------------------------------------------------
source_id   int64      Unique source identifier (unique within a particular Data Release)
ref_epoch float64   yr                                                    Reference epoch
ra float64  deg                                                    Right ascension
dec float64  deg                                                        Declination
parallax float64  mas                                                           Parallax
Jobid: 1609260407863O
Phase: COMPLETED
Owner: None
Output file: async_20201229114648.vot
Results: None

And here are the results.

results2 = job2.get_results()
results2

Table length=3000
yrdegdegmas
int64float64float64float64float64
67585097575941414402015.5290.4899010727352-30.343172184207830.48023816159705535
67585270362508492802015.5290.35703708993327-30.193657826181596-0.2763960334229464
67585644582982429442015.5290.77376294142846-29.7653684392252380.5907906528352993
67585586128421557762015.5290.67550630386245-29.921133960781690.2858563565989917
67585563794599219202015.5290.5454006212404-29.900709054816964-1.0012355835832834
67585063559800106242015.5290.586950869878-30.371376421675690.3769870991981157
67585641834144084482015.5290.7133096669958-29.7817436736796070.9376387942856869
67585297893226531842015.5290.35449465190385-30.10231481465004-1.3305133038319952
67585074855553360642015.5290.6745961953205-30.359651480610280.057302711920868686
...............
46604648940382954242015.579.74509945051724-66.496560127377590.1631257821260955
46604635196504801282015.580.06150680482719-66.42891601595986-0.5822669736733328
46604876745540710402015.580.03536662790846-66.33136792885134-1.4978065466579966
46604947183008867842015.579.6181997343967-66.174723353182120.1488652133861382
46604764732741593602015.579.12508291180023-66.356448447572750.10671540523101529
46604710100939402242015.579.1908378165463-66.489168402940960.8152150079365807
46604983948071677442015.580.23431532382547-66.08688681815092-0.14060596426163252
46604987383917228802015.580.20205898946918-66.043923132951260.1894642263668475
46604740680910969602015.578.93686107060341-66.457955767458960.4665513634059366
46604373031626977282015.580.63872120225027-66.46144560233674-1.1319163320535979

You might notice that some values of parallax are negative. As this FAQ explains, “Negative parallaxes are caused by errors in the observations.” Negative parallaxes have “no physical meaning,” but they can be a “useful diagnostic on the quality of the astrometric solution.”

### Exercise¶

The clauses in a query have to be in the right order. Go back and change the order of the clauses in query2 and run it again.

The query should fail, but notice that you don’t get much useful debugging information.

For this reason, developing and debugging ADQL queries can be really hard. A few suggestions that might help:

• Whenever possible, start with a working query, either an example you find online or a query you have used in the past.

• Make small changes and test each change before you continue.

• While you are debugging, use TOP to limit the number of rows in the result. That will make each test run faster, which reduces your development time.

• Launching test queries synchronously might make them start faster, too.

# Solution

## Operators¶

In a WHERE clause, you can use any of the SQL comparison operators; here are the most common ones:

Symbol

Operation

>

greater than

<

less than

>=

greater than or equal

<=

less than or equal

=

equal

!= or <>

not equal

Most of these are the same as Python, but some are not. In particular, notice that the equality operator is =, not ==. Be careful to keep your Python out of your ADQL!

You can combine comparisons using the logical operators:

• AND: true if both comparisons are true

• OR: true if either or both comparisons are true

Finally, you can use NOT to invert the result of a comparison.

### Exercise¶

Read about SQL operators here and then modify the previous query to select rows where bp_rp is between -0.75 and 2.

# Solution

# Here's a solution using > and < operators

query = """SELECT TOP 10
source_id, ref_epoch, ra, dec, parallax
WHERE parallax < 1
AND bp_rp > -0.75 AND bp_rp < 2
"""

# And here's a solution using the BETWEEN operator

query = """SELECT TOP 10
source_id, ref_epoch, ra, dec, parallax
WHERE parallax < 1
AND bp_rp BETWEEN -0.75 AND 2
"""

This Hertzsprung-Russell diagram shows the BP-RP color and luminosity of stars in the Gaia catalog (Copyright: ESA/Gaia/DPAC, CC BY-SA 3.0 IGO).

Selecting stars with bp-rp less than 2 excludes many class M dwarf stars, which are low temperature, low luminosity. A star like that at GD-1’s distance would be hard to detect, so if it is detected, it it more likely to be in the foreground.

## Cleaning up¶

Asynchronous jobs have a jobid.

job1.jobid, job2.jobid

(None, '1609260407863O')

Which you can use to remove the job from the server.

Gaia.remove_jobs([job2.jobid])

Removed jobs: '['1609260407863O']'.

If you don’t remove it job from the server, it will be removed eventually, so don’t feel too bad if you don’t clean up after yourself.

## Formatting queries¶

So far the queries have been string “literals”, meaning that the entire string is part of the program. But writing queries yourself can be slow, repetitive, and error-prone.

It is often a good idea to write Python code that assembles a query for you. One useful tool for that is the string format method.

As an example, we’ll divide the previous query into two parts; a list of column names and a “base” for the query that contains everything except the column names.

Here’s the list of columns we’ll select.

columns = 'source_id, ra, dec, pmra, pmdec, parallax, radial_velocity'

And here’s the base; it’s a string that contains at least one format specifier in curly brackets (braces).

query3_base = """SELECT TOP 10
{columns}
WHERE parallax < 1
AND bp_rp BETWEEN -0.75 AND 2
"""

This base query contains one format specifier, {columns}, which is a placeholder for the list of column names we will provide.

To assemble the query, we invoke format on the base string and provide a keyword argument that assigns a value to columns.

query3 = query3_base.format(columns=columns)

The result is a string with line breaks. If you display it, the line breaks appear as \n.

query3

'SELECT TOP 10 \nsource_id, ra, dec, pmra, pmdec, parallax, radial_velocity\nFROM gaiadr2.gaia_source\nWHERE parallax < 1\n  AND bp_rp BETWEEN -0.75 AND 2\n'

But if you print it, the line breaks appear as… line breaks.

print(query3)

SELECT TOP 10
source_id, ra, dec, pmra, pmdec, parallax, radial_velocity
WHERE parallax < 1
AND bp_rp BETWEEN -0.75 AND 2

Notice that the format specifier has been replaced with the value of columns.

Let’s run it and see if it works:

job3 = Gaia.launch_job(query3)
print(job3)

--------------- ------- -------- ------------------------------------------------------------------ -----
source_id   int64          Unique source identifier (unique within a particular Data Release)     0
ra float64      deg                                                    Right ascension     0
dec float64      deg                                                        Declination     0
pmra float64 mas / yr                         Proper motion in right ascension direction     0
pmdec float64 mas / yr                             Proper motion in declination direction     0
parallax float64      mas                                                           Parallax     0
Jobid: None
Phase: COMPLETED
Owner: None
Output file: sync_20201229114653.xml.gz
Results: None

results3 = job3.get_results()
results3

Table length=10
degdegmas / yrmas / yrmaskm / s
int64float64float64float64float64float64float64
466046637150777484879.49100199261952-66.414117776726682.6073927760139983-1.1968907801968522-0.13534455558687877--
466049870403198400080.19346436358076-66.049076752328561.644062563874402-1.954975321006136-0.1774586376397--
466045894982481702479.68666778641992-66.522209596625571.47776972758805810.3328961930362448-0.030149510331454386--
466045197477392652879.89972073493868-66.750378580711912.4635048563021065-0.7474574729840501-0.005219591141134416--
466047403373135705678.93214036467484-66.464836445916671.2369118132455594-1.71871430346386770.6018032937392243--
466045410507787468879.91045649235578-66.681068461359711.47478554075336770.022464361420165786-0.19949109843083218--
466049485575374054479.63101881952036-66.150633158135361.62278632632221130.204908932619051520.0013030724292998944--
466043953657130816080.24817961417894-66.463032706640921.8652105429518493-0.5618584442373111-0.10578442360531497--
466047018544046656079.25942306061864-66.511207775496053.3735611720039890.25326243152876704-0.6932631338638376--
466044146069036556880.38862155241748-66.344044941350461.6588478675146066-0.09646951910977163-0.2100859303087445--

Good so far.

### Exercise¶

This query always selects sources with parallax less than 1. But suppose you want to take that upper bound as an input.

Modify query3_base to replace 1 with a format specifier like {max_parallax}. Now, when you call format, add a keyword argument that assigns a value to max_parallax, and confirm that the format specifier gets replaced with the value you provide.

# Solution

query4_base = """SELECT TOP 10
{columns}
WHERE parallax < {max_parallax} AND
bp_rp BETWEEN -0.75 AND 2
"""

query4 = query4_base.format(columns=columns,
max_parallax=0.5)
print(query)

SELECT TOP 10
source_id, ref_epoch, ra, dec, parallax
WHERE parallax < 1
AND bp_rp BETWEEN -0.75 AND 2

## Summary¶

This notebook demonstrates the following steps:

1. Making a connection to the Gaia server,

2. Exploring information about the database and the tables it contains,

3. Writing a query and sending it to the server, and finally

4. Downloading the response from the server as an Astropy Table.

In the next lesson we will extend these queries to select a particular region of the sky.

## Best practices¶

• If you can’t download an entire dataset (or it’s not practical) use queries to select the data you need.

• Read the metadata and the documentation to make sure you understand the tables, their columns, and what they mean.

• Develop queries incrementally: start with something simple, test it, and add a little bit at a time.

• Use ADQL features like TOP and COUNT to test before you run a query that might return a lot of data.

• If you know your query will return fewer than 3000 rows, you can run it synchronously, which might complete faster (but it doesn’t seem to make much difference). If it might return more than 3000 rows, you should run it asynchronously.

• ADQL and SQL are not case-sensitive, so you don’t have to capitalize the keywords, but you should.

• ADQL and SQL don’t require you to break a query into multiple lines, but you should.

Jupyter notebooks can be good for developing and testing code, but they have some drawbacks. In particular, if you run the cells out of order, you might find that variables don’t have the values you expect.

There are a few things you can do to mitigate these problems:

• Make each section of the notebook self-contained. Try not to use the same variable name in more than one section.

• Keep notebooks short. Look for places where you can break your analysis into phases with one notebook per phase.

One of the other tables we’ll use is gaiadr2.panstarrs1_original_valid. Use load_table to get the metadata for this table. How many columns are there and what are their names?


text:
# Math Help - [SOLVED] Tignometry Help

1. ## [SOLVED] Tignometry Help

Question:
Use the symmetric and periodic properties of the sine, cosine and tangent functions to establish the following results.

$(a) \sin (90-\theta)^o = \cos \theta^o$
$(b) \sin (270 + \theta)^o = -\cos \theta^o$
$(c) \tan (\theta - 180)^o = \tan \theta^o$

Attempt:
I don't know how to start

2. Hello,

For the first one, there are two ways imho :
- sketch a trigonometrical circle and you will see that this relation is true
- use the formula sin(a-b)=sin(a)cos(b)-sin(b)cos(a)

For the second one, use the periodicity :
- the sine function is 360-periodic. Hence $sin(270+\theta)=sin((270+\theta)-360)=sin(\theta-90)$

- As the sine function is an odd function, sin(-x)=-sin(x)
- $sin(\theta-90)=-sin(90-\theta)$ and you retrieve the first relation


text:
Version:

This documentation is not for the latest stable Salvus version.

This tutorial is presented as Python code running inside a Jupyter Notebook, the recommended way to use Salvus. To run it yourself you can copy/type each individual cell or directly download the full notebook, including all required files.

This notebook teaches how to compute gradients with respect to a mesh's material parameters using only SalvusFlow and not SalvusProject. SalvusProject internally handles many of the complexities but the approach presented here grants a bit more flexibility.

Copy
# Variable used in the notebook to determine which site
# is used to run the simulations.
import os
SALVUS_FLOW_SITE_NAME = os.environ.get("SITE_NAME", "local")
%matplotlib inline

import pathlib

import numpy as np
import matplotlib.pyplot as plt
import xarray as xr

import salvus.namespace as sn

## Step 1: Generate target and initial model

This example is purely synthetic so we'll generate a target model with a Gaussian blob perturbation in $v_s$ in the center of the domain. This will be used to generate the "observed data" in this case. The initial model (the synthetics for cause of this tutorial) is a homogeneous model without the blob.

# Use xarray to define the blob in vs.
def get_target_model():
x = np.linspace(-500, 500, 200)
y = np.linspace(-500, 500, 200)
xx, yy = np.meshgrid(x, y, indexing="ij")

# Simple gaussian blob around the center of the domain.
d = (xx ** 2 + yy ** 2) ** 0.5
vs = np.exp(-0.5 * (d / 100.0) ** 2) * 250.0 + 1000.0

ds = xr.Dataset(
data_vars={"vs": (["x", "y"], vs)}, coords={"x": x, "y": y}
)

return ds

target_model = get_target_model()
target_model.vs.plot()
<matplotlib.collections.QuadMesh at 0x7f7b17381b90>

Create a mesh - this is directly used as the initial model.

# Build mesh for the initial model.
initial_model_mesh = sn.simple_mesh.CartesianHomogeneousIsotropicElastic2D(
vp=1500.0,
vs=1000.0,
rho=2000.0,
x_max=float(target_model.x.data.ptp()),
y_max=float(target_model.y.data.ptp()),
max_frequency=10.0,
tensor_order=1,
).create_mesh()

# Shift.
initial_model_mesh.points[:, 0] += target_model.x.data.min()
initial_model_mesh.points[:, 1] += target_model.y.data.min()

Now interpolate the blob onto it thus creating the target mesh. Have a look at the perturbation in $v_s$.

target_model_mesh = initial_model_mesh.copy()

target_model_mesh = sn.toolbox.interpolate_cartesian(
data=target_model, mesh=target_model_mesh, extrapolate=True
)
target_model_mesh
/miniconda/envs/salvus/lib/python3.7/_collections_abc.py:702: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
return len(self._mapping)
/miniconda/envs/salvus/lib/python3.7/_collections_abc.py:720: MatplotlibDeprecationWarning: The global colormaps dictionary is no longer considered public API.
yield from self._mapping

When run interactively in a Jupyter Notebook, this output cell contains a widget. Click the button below to load a preview of it. Please note that most interactive functionality does not work without a running Python kernel.
<salvus.mesh.unstructured_mesh.UnstructuredMesh at 0x7f7b1088d890>

This is also the starting point for setting up a project for the inversion. We can directly create the project from the model defined above. We will call this model the true_model.

## Step 2: Setup the simulation with a single source and a ring of receivers

# Explosive source in one corner of the mesh.
src = sn.simple_config.source.cartesian.MomentTensorPoint2D(
mxx=10000.0,
myy=10000.0,
mxy=0.0,
x=300.0,
y=300.0,
source_time_function=sn.simple_config.stf.Ricker(center_frequency=5.0),
)

# A centered ring of receivers.
)

# Create two simulation objects.
# (1) The target will serve as observed data.
w_target = sn.simple_config.simulation.Waveform(
)
w_target.physics.wave_equation.end_time_in_seconds = 1.0
w_target

# (2) The initial model will be used to compute synthetics.
w_initial = sn.simple_config.simulation.Waveform(
)
w_initial.physics.wave_equation.end_time_in_seconds = 1.0

# Must store the checkpoints for the subsequent adjoint run.
w_initial.output.volume_data.format = "hdf5"
w_initial.output.volume_data.filename = "output.h5"
w_initial.output.volume_data.sampling_interval_in_time_steps = (
"auto-for-checkpointing"
)

## Step 3: Run Simulations

For the target/observed model and the initial/current model.

# Create the "observed data".
j_target = sn.api.run(
input_file=w_target,
site_name=SALVUS_FLOW_SITE_NAME,
output_folder="output_target",
overwrite=True,
)

# Forward run for the synthetics.
j_initial = sn.api.run(
input_file=w_initial,
site_name=SALVUS_FLOW_SITE_NAME,
output_folder="output_initial",
overwrite=True,
# Don't delete the remote files as they are needed
delete_remote_files=False,
)
Job job_2009171308302123_07f14f4a4e running on local with 4 rank(s).
Site information:
* Salvus version: 0.11.16
* Floating point size: 32
-> Current Task: Time loop complete* Downloaded 341.8 KB of results to output_target.
* Total run time: 1.33 seconds.
* Pure simulation time: 0.78 seconds.
Job job_2009171308495663_bdd42f15ac running on local with 4 rank(s).
Site information:
* Salvus version: 0.11.16
* Floating point size: 32
-> Current Task: Time loop complete* Downloaded 342.3 KB of results to output_initial.
* Total run time: 3.72 seconds.
* Pure simulation time: 1.81 seconds.

# Load the data from both runs, and visualize.
#
# Note that this is a vectorial receiver and thus the component
# has to be specified.
target_event = j_target.get_as_event()
initial_event = j_initial.get_as_event()

initial_event.plot(receiver_field="displacement", component="X")