User:Eml5526.s11.team5/sub5

= 5.1 =

P.D.E

 * {| style="width:100%" border="0"

$$ \frac{\partial}{\partial x} \left[ (2+3x) \frac{\partial u}{\partial x} \right] + 5x = 0 \ \ \forall x \in \left] 0,1 \right[ $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *  $$ \displaystyle                                                                          (Eq. 5.1.1)
 * }

Boundary Conditions

 * {| style="width:100%" border="0"

$$u \ at \ {x=0} \ is \ 4 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *  $$ \displaystyle                                                                          (Eq. 5.1.2)
 * }


 * {| style="width:100%" border="0"

$$  \frac{\partial u}{\partial x} \  at \ x = 1 \ is \ \frac{12}{5} $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *  $$ \displaystyle                                                                          (Eq. 5.1.3)
 * }

Find: Approximate solution using weak form and its comparison to exact solution

 * Consider basis function satisfying constraint basis solution $$ \boldsymbol{F_{I}}, \ where \ I = p, f , e $$


 * Compute $$ u_n^h (x = 0.5) \quad $$ for $$ n = 2, 4, 6..... \quad $$ until convergence of $$ u^h (0.5)\ to \ 0(10^{-6}) $$


 * error


 * {| style="width:100%" border="0"

$$ e_n(0.5) = u_n(0.5) - u_n^h(0.5) \quad $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *  $$ \displaystyle                                                                       (Eq. 5.1.4)
 * }


 * Plot $$ e_n(0.5) \quad $$ vs. $$ n \quad $$

Solution
Strong form


 * {| style="width:100%" border="0"

$$ \frac{\partial}{\partial x} \left[ (2+3x) \frac{\partial u}{\partial x} \right] + 5x = 0 \ \ \forall x \in \left] 0,1 \right[ $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *  $$ \displaystyle                                                                       (Eq. 5.1.5)
 * }

Boundary conditions


 * {| style="width:100%" border="0"

$$u \ at \ {x=0} \ is \ 4 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *  $$ \displaystyle                                                                          (Eq. 5.1.6)
 * }


 * {| style="width:100%" border="0"

$$  \frac{\partial u}{\partial x} \  at \ x = 1 \ is \ \frac{12}{5} $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *  $$ \displaystyle                                                                          (Eq. 5.1.7)
 * }

Finding of the Exact Solution
The exact solution of the displacement can be found as follows.

Integrating eqn 5.2.5, we get,


 * {| style="width:100%" border="0"

$$ (2+3x)\frac{\partial u}{\partial x} + \frac{5x^2}{2} = c_1 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 *  $$ \displaystyle                                                                          (Eq. 5.1.8)
 * }

Substituting the natural B.C in the above equation, we get,

<span id="(1)">
 * {| style="width:100%" border="0"

$$ c_1 = 14.5 \quad$$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 5.1.9)
 * }

Therefore eqn 5.1.8 becomes,

<span id="(1)">
 * {| style="width:100%" border="0"

$$ (2+3x)\frac{\partial u}{\partial x} + \frac{5x^2}{2} = 14.5 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 5.1.10)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$\frac{\partial u}{\partial x} = \frac{14.5-\frac{5x^2}{2}}{(2+3x)} $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 5.1.11)
 * }

Integrating the above equation, we get,Reference : www.wolframalpha.com

<span id="(1)">
 * {| style="width:100%" border="0"

$$u = \frac{1}{108}(-45x^2 + 60x + 482 log (2+3x) +60 ) + c_2 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 5.1.12)
 * }

Substituting the Essential B.C in the above equation, we get,

<span id="(1)">
 * {| style="width:100%" border="0"

$$c_2 = 0.3509542497 \quad$$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 5.1.13)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$u = \frac{1}{108}(-45x^2 + 60x + 482 log (2+3x) +60 ) + 0.3509542497 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 5.1.14)
 * }

Finding of the Approximate Solution using Weak form
Weak Form

<span id="(1)">
 * {| style="width:100%" border="0"

$$  \sum_{i} c_i \left [ \sum_{j} \left \{ \tilde{K_{ij}} \tilde{d_j} \right \} - \tilde{F_i} \right ] = 0 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 5.1.15)
 * }

where,

<span id="(1)">
 * {| style="width:100%" border="0"

$$\tilde{K_{ij}} = \int_{0}^{1} b_i' \ a_2 \ b_j ' \ dx $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 5.1.16)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ \tilde{F_i} = b_i (\alpha )h + \int_{0}^{1} b_i f dx $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 5.1.17)
 * }

Finding of the Approximate Solution using Polynomial Basis Function satisfying Constraint Breaking Solution in Weak form
<span id="(1)">
 * {| style="width:100%" border="0"

$$b_j = x^j \ where \ j = 0,1,2,... $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 5.1.18)
 * }

Let us try for n = 2;

<span id="(1)">
 * {| style="width:100%" border="0"

$$ b_j = [ 1 \quad (x) \quad (x)^2]$$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 5.1.19)
 * }

We know that,

<span id="(1)">
 * {| style="width:100%" border="0"

$$ u^h (x) = \sum_{0}^{2} d_j b_j(x) = d_0 + (x)d_1 + (x)^2 d_2 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 5.1.20)
 * }

Substituting the essential B.C;

<span id="(1)">
 * {| style="width:100%" border="0"

$$ d_0 = 4 \quad $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 5.1.21)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ b_j' = [ 0 \quad 1 \quad 2(x)] $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 5.1.22)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$K_{01} = K_{10} = K_{02} = K_{20} = 0 \quad $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 5.1.23)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ K_{11} = \int_{0}^{1}(1)(2+3x)(1)dx = 3.5 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 5.1.24)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ K_{12} = K_{21} = \int_{0}^{1}(1)(2+3x)2(x)(1)dx = 2\int_{0}^{1}(3x^2-x-2)dx = 4 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 5.1.25)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$K_{22} = \int_{0}^{1}2(x)(2+3x)2(x)dx = 4\int_{0}^{1}(2+3x)(x)dx = 5.6667 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 5.1.26)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ F_0 = b_0(\alpha)(12) + \int_{0}^{1}b_0 (5x)dx = (1)(12) + \int_{0}^{1} 5xdx = 14.5 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 5.1.27)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$F_1 = b_1(\alpha)(12) + \int_{0}^{1}b_1 (5x)dx = (1)(12) + \int_{0}^{1} 5x(x)dx = 13.6667 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 5.1.28)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ F_2 = b_2(\alpha)(12) + \int_{0}^{1}b_2 (5x)dx = (12) + \int_{0}^{1} 5x(x)^2dx = 13.25 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 5.1.29)
 * }

Hence we obtain 3 linearly independent equations, (i.e) one from the essential boundary condition and the other two from the weak form. These three equations can be written in matrix form as follows.

<span id="(1)">
 * {| style="width:100%" border="0"

$$ K.d = F \quad $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ \begin{bmatrix} 1 & 0 & 0\\ 0& 3.5 & 4\\  0&  4 & 5.6667  \end{bmatrix} \begin{Bmatrix} d_0\\ d_1 \\ d_2 \end{Bmatrix} = \begin{Bmatrix} 4\\ 13.6667 \\ 13.25 \end{Bmatrix}$$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 5.1.30)
 * }

Solving for $$ d_0 \quad $$, $$ d_1 \quad$$ and $$ d_2 \quad$$ , we get,

<span id="(1)">
 * {| style="width:100%" border="0"

$$ d_0 = 4; \qquad d_1 = 6.3769; \qquad d_2 = -2.1631$$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 5.1.31)
 * }

Therefore the displacement equation is given by ,

<span id="(1)">
 * {| style="width:100%" border="0"

$$ u^h(x) = 4+6.3769(x) -2.1631(x)^2 \quad $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 5.1.32)
 * }

At $$ x = 0.5 \quad $$ ,

<span id="(1)">
 * {| style="width:100%" border="0"

$$ u^h(0.5) = 4+6.3769(0.5) -2.1631(0.5)^2 = 6.647675 \quad $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 5.1.33)
 * }

The exact solution can be found from the eqn. 5.1.14. It is given by ,

<span id="(1)">
 * {| style="width:100%" border="0"

$$ u(0.5) = 6.671155 \quad$$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

Hence the error we get for $$ n = 2 $$ is ,

<span id="(1)">
 * {| style="width:100%" border="0"

$$u^h(0.5) - u(0.5) = -0.023480 \quad$$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 5.1.34)
 * }

MATLAB code with output
The K matrix and F matrix which we get for the convergence upto 10^-6 is as follows:

K =

1        0         0         0         0         0         0         0         0         0    3.5000    4.0000    4.2500    4.4000    4.5000    4.5714    4.6250    4.6667         0    4.0000    5.6667    6.6000    7.2000    7.6190    7.9286    8.1667    8.3556         0    4.2500    6.6000    8.1000    9.1429    9.9107   10.5000   10.9667   11.3455         0    4.4000    7.2000    9.1429   10.5714   11.6667   12.5333   13.2364   13.8182         0    4.5000    7.6190    9.9107   11.6667   13.0556   14.1818   15.1136   15.8974         0    4.5714    7.9286   10.5000   12.5333   14.1818   15.5455   16.6923   17.6703         0    4.6250    8.1667   10.9667   13.2364   15.1136   16.6923   18.0385   19.2000         0    4.6667    8.3556   11.3455   13.8182   15.8974   17.6703   19.2000   20.5333

F =

[4  13.6667   13.2500   13.0000   12.8333   12.7143   12.6250   12.5556   12.5000]^T

d =

[4.0000   7.2498   -5.4295    4.9208   -5.0302    4.5192   -3.0193    1.2522   -0.2348]^T

The outputs are as follow:

@x=0.5 Exact U=6.671156e+000

U(weak form) for n=2 = 6.647645e+000

Error at n=2 is 2.351072e-002

@x=0.5 Exact U=6.671156e+000

U(weak form) for n=3 = 6.678544e+000

Error at n=3 is 7.388478e-003

@x=0.5 Exact U=6.671156e+000

U(weak form) for n=4 = 6.671260e+000

Error at n=4 is 1.047460e-004

@x=0.5 Exact U=6.671156e+000

U(weak form) for n=5 = 6.670884e+000

Error at n=5 is 2.712998e-004

@x=0.5 Exact U=6.671156e+000

U(weak form) for n=6 = 6.671152e+000

Error at n=6 is 3.583888e-006

@x=0.5 Exact U=6.671156e+000

U(weak form) for n=7 = 6.671166e+000

Error at n=7 is 1.009645e-005

@x=0.5 Exact U=6.671156e+000

U(weak form) for n=8 = 6.671156e+000

Error at n=8 is 9.634049e-008



Finding of the Approximate Solution using Fourier Basis Function satisfying Constraint Breaking Solution in Weak form
Let us take $$ n = 1 \quad $$ Referencewww.wolframalpha.com

Referencewww.wolframalpha.com

Referencewww.wolframalpha.com

Referencewww.wolframalpha.com

Referencewww.wolframalpha.com

MATLAB code with output
The K matrix and F matrix which we get for the convergence upto 10^-6 is as follows:

K =

1        0         0         0         0         0         0         0         0         0    1.1444   -1.3612    3.2564   -0.5408    3.7008    2.3208    1.2335    4.7249         0   -1.3612    2.3556   -4.1866    2.4001   -5.8462   -0.0674   -4.7299   -2.9460         0    3.2564   -4.1866    9.5121   -2.3503   11.6194    5.4649    5.7892   12.8497         0   -0.5408    2.4001   -2.3503    4.4879   -5.4725    5.3501   -8.9290    3.8999         0    3.7008   -5.8462   11.6194   -5.4725   16.8127    2.2054   14.2443   12.2128         0    2.3208   -0.0674    5.4649    5.3501    2.2054   14.6873   -9.6620   19.6948         0    1.2335   -4.7299    5.7892   -8.9290   14.2443   -9.6620   23.4828   -3.0985         0    4.7249   -2.9460   12.8497    3.8999   12.2128   19.6948   -3.0985   32.5172

F =

[4  -6.1075   11.6035  -18.9907   13.0886  -27.2503    3.4218  -23.8065   -8.5011]^T

d =

[4.0000   47.4091   48.5755   -3.4993  -30.6212   -4.1185    7.2053    0.9062   -0.4248]^T

The outputs are as follow:

@x=0.5 Exact U=6.671156e+000

U(weak form) for n=1 = 6.619837e+000

Error at n=1 is 5.131875e-002

@x=0.5 Exact U=6.671156e+000

U(weak form) for n=2 = 6.671392e+000

Error at n=2 is 2.367707e-004

@x=0.5 Exact U=6.671156e+000

U(weak form) for n=3 = 6.671145e+000

Error at n=3 is 1.055939e-005

@x=0.5 Exact U=6.671156e+000

U(weak form) for n=4 = 6.671156e+000

Error at n=4 is 3.459119e-007



Finding of the Approximate Solution using Exponential Basis Function satisfying Constraint Breaking Solution in Weak form
For $$ n = 2 \quad $$ Referencewww.wolframalpha.com

Referencewww.wolframalpha.com

Referencewww.wolframalpha.com

Referencewww.wolframalpha.com

Referencewww.wolframalpha.com

MATLAB code with output
The K matrix and F matrix which we get for the convergence upto 10^-6 is as follows:

K =

1.0e+008 *

1        0         0         0         0         0         0         0         0         0    0.0000    0.0000    0.0000    0.0000    0.0000    0.0000    0.0001    0.0003         0    0.0000    0.0000    0.0000    0.0000    0.0001    0.0002    0.0006    0.0017         0    0.0000    0.0000    0.0000    0.0001    0.0003    0.0008    0.0022    0.0062         0    0.0000    0.0000    0.0001    0.0003    0.0008    0.0025    0.0072    0.0206         0    0.0000    0.0001    0.0003    0.0008    0.0026    0.0077    0.0225    0.0649         0    0.0000    0.0002    0.0008    0.0025    0.0077    0.0232    0.0682    0.1973         0    0.0001    0.0006    0.0022    0.0072    0.0225    0.0682    0.2014    0.5858         0    0.0003    0.0017    0.0062    0.0206    0.0649    0.1973    0.5858    1.7106

F =

1.0e+004 * [4  0.0023    0.0085    0.0249    0.0692    0.1885    0.5107    1.3817    3.7387]^T

d =

[4.0000  171.4273 -268.6688  254.1458 -153.5759   59.9076  -14.6354    2.0398   -0.1239]^T

The outputs are as follow:

@x=0.5 Exact U=6.671156e+000

U(weak form) for n=2 = 6.328507e+000

Error at n=2 is 3.426491e-001

@x=0.5 Exact U=6.671156e+000

U(weak form) for n=3 = 6.661733e+000

Error at n=3 is 9.422741e-003

@x=0.5 Exact U=6.671156e+000

U(weak form) for n=4 = 6.681880e+000

Error at n=4 is 1.072431e-002

@x=0.5 Exact U=6.671156e+000

U(weak form) for n=5 = 6.670738e+000

Error at n=5 is 4.179779e-004

@x=0.5 Exact U=6.671156e+000

U(weak form) for n=6 = 6.669725e+000

Error at n=6 is 1.431086e-003

@x=0.5 Exact U=6.671156e+000

U(weak form) for n=7 = 6.670941e+000

Error at n=7 is 2.143177e-004

@x=0.5 Exact U=6.671156e+000

U(weak form) for n=8 = 6.671171e+000

Error at n=8 is 1.566386e-005

>>

Convergence Plot


--Eml5526.s11.team5.srv 02:08, 24 March 2011 (UTC)

= 5.2 =

P.D.E
<span id="(1)">
 * {| style="width:100%" border="0"

$$ \frac{\partial}{\partial x} \left[ (2+3x) \frac{\partial u}{\partial x} \right] + 5x = 0 \ \ \forall x \in \left] 0,1 \right[ $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 5.2.1)
 * }

Boundary Conditions
<span id="(1)">
 * {| style="width:100%" border="0"

$$u \ at \ {x=1} \ is \ 4 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 5.2.2)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$  \frac{\partial u}{\partial x} \  at \ x = 0 \ is \ -6 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 5.2.3)
 * }

Find: Approximate solution using weak form and its comparison to exact solution
 Find approximate solution using weak form and compare it to exact one


 * Consider basis function satisfying constraint basis solution $$ \boldsymbol{F_{I}}, \ where \ I = p, f , e $$


 * Compute $$ u_n^h (x = 0.5) \quad $$ for $$ n = 2, 4, 6..... \quad $$ until convergence of $$ u^h (0.5)\ to \ 0(10^{-6}) $$


 * error

<span id="(1)">
 * {| style="width:100%" border="0"

$$ e_n(0.5) = u_n(0.5) - u_n^h(0.5) \quad $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                       (Eq. 5.2.4)
 * }


 * Plot $$ e_n(0.5) \quad $$ vs. $$ n \quad $$

Solution
Strong form

<span id="(1)">
 * {| style="width:100%" border="0"

$$ \frac{\partial}{\partial x} \left[ (2+3x) \frac{\partial u}{\partial x} \right] + 5x = 0 \ \ \forall x \in \left] 0,1 \right[ $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                       (Eq. 5.2.5)
 * }

Boundary conditions

<span id="(1)">
 * {| style="width:100%" border="0"

$$u \ at \ {x=1} \ is \ 4 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 5.2.6)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$  \frac{\partial u}{\partial x} \  at \ x = 0 \ is \ -6 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 5.2.7)
 * }

Finding of the Exact Solution
The exact solution of the displacement can be found as follows.

Differentiating eqn 5.2.5, we get,

<span id="(1)">
 * {| style="width:100%" border="0"

$$ (2+3x)\frac{\partial u}{\partial x} + \frac{5x^2}{2} = c_1 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 5.2.8)
 * }

Substituting the natural B.C in the above equation, we get,

<span id="(1)">
 * {| style="width:100%" border="0"

$$ c_1 = -12 \quad$$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 5.2.9)
 * }

Therefore eqn 5.2.8 becomes,

<span id="(1)">
 * {| style="width:100%" border="0"

$$ (2+3x)\frac{\partial u}{\partial x} + \frac{5x^2}{2} = -12 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 5.2.10)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$\frac{\partial u}{\partial x} = \frac{-12-\frac{5x^2}{2}}{(2+3x)} $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 5.2.11)
 * }

Integrating the above equation, we get,

<span id="(1)">
 * {| style="width:100%" border="0"

$$u = \frac{1}{108}(-45x^2 + 60x - 472 log (2+3x) +60 ) + c_2 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 5.2.12)
 * }

Substituting the Essential B.C in the above equation, we get,

<span id="(1)">
 * {| style="width:100%" border="0"

$$c_2 = 10.3394 \quad$$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 5.2.13)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$u = \frac{1}{108}(-45x^2 + 60x - 472 log (2+3x) +60 ) + 10.3393953 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 5.2.14)
 * }

Finding of the Approximate Solution using Weak form
Weak Form

<span id="(1)">
 * {| style="width:100%" border="0"

$$  \sum_{i} c_i \left [ \sum_{j} \left \{ \tilde{K_{ij}} \tilde{d_j} \right \} - \tilde{F_i} \right ] = 0 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 5.2.15)
 * }

where,

<span id="(1)">
 * {| style="width:100%" border="0"

$$\tilde{K_{ij}} = \int_{0}^{1} b_i' \ a_2 \ b_j ' \ dx $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 5.2.16)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ \tilde{F_i} = b_i (\alpha )h + \int_{0}^{1} b_i f dx $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 5.2.17)
 * }

Finding of the Approximate Solution using Polynomial Basis Function satisfying Constraint Breaking Solution in Weak form
<span id="(1)">
 * {| style="width:100%" border="0"

$$b_j = (x-1)^j \ where \ j = 0,1,2,... $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 5.2.18)
 * }

Let us try for n = 2;

<span id="(1)">
 * {| style="width:100%" border="0"

$$ b_j = [ 1 \quad (x-1) \quad (x-1)^2]$$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 5.2.19)
 * }

We know that,

<span id="(1)">
 * {| style="width:100%" border="0"

$$ u^h (x) = \sum_{0}^{2} d_j b_j(x) = d_0 + (x-1)d_1 + (x-1)^2 d_2 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 5.2.20)
 * }

Substituting the essential B.C;

<span id="(1)">
 * {| style="width:100%" border="0"

$$ d_0 = 4 \quad $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 5.2.21)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ b_j' = [ 0 \quad 1 \quad 2(x-1)] $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 5.2.22)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$K_{01} = K_{10} = K_{02} = K_{20} = 0 \quad $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 5.2.23)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ K_{11} = \int_{0}^{1}(1)(2+3x)(1)dx = 3.5 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 5.2.24)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ K_{12} = K_{21} = \int_{0}^{1}(1)(2+3x)2(x-1)(1)dx = 2\int_{0}^{1}(3x^2-x-2)dx = -3 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 5.2.25)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$K_{22} = \int_{0}^{1}2(x-1)(2+3x)2(x-1)dx = 4\int_{0}^{1}(2+3x)(x-1)dx = 3.6667 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 5.2.26)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ F_0 = b_0(\alpha)(12) + \int_{0}^{1}b_0 (5x)dx = (1)(12) + \int_{0}^{1} 5xdx = 14.5 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 5.2.27)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$F_1 = b_1(\alpha)(12) + \int_{0}^{1}b_1 (5x)dx = (-1)(12) + \int_{0}^{1} 5x(x-1)dx = -12.833 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 5.2.28)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ F_2 = b_2(\alpha)(12) + \int_{0}^{1}b_2 (5x)dx = (12) + \int_{0}^{1} 5x(x-1)^2dx = 12.4166 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 5.2.29)
 * }

Hence we obtain 3 linearly independent equations, (i.e) one from the essential boundary condition and the other two from the weak form. These three equations can be written in matrix form as follows.

<span id="(1)">
 * {| style="width:100%" border="0"

$$ K.d = F \quad $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ \begin{bmatrix} 1 & 0 & 0\\ 0& 3.5 & -3\\  0&  -3 & 3.6667  \end{bmatrix} \begin{Bmatrix} d_0\\ d_1 \\ d_2 \end{Bmatrix} = \begin{Bmatrix} 4\\ -12.833 \\ 12.4166 \end{Bmatrix}$$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 5.2.30)
 * }

Solving for $$ d_0 \quad $$, $$ d_1 \quad$$ and $$ d_2 \quad$$ , we get,

<span id="(1)">
 * {| style="width:100%" border="0"

$$ d_0 = 4; \qquad d_1 = -2.55768; \qquad d_2 = 1.2937$$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 5.2.31)
 * }

Therefore the displacement equation is given by ,

<span id="(1)">
 * {| style="width:100%" border="0"

$$ u^h(x) = 4-2.55768(x-1) + 1.2937(x-1)^2 \quad $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 5.2.32)
 * }

At $$ x = 0.5 \quad $$ ,

<span id="(1)">
 * {| style="width:100%" border="0"

$$ u^h(0.5) = 4-2.55768(0.5-1) + 1.2937(0.5-1)^2 = 5.602265 \quad $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 5.2.33)
 * }

The exact solution can be found from the eqn. 5.2.14. It is given by ,

<span id="(1)">
 * {| style="width:100%" border="0"

$$ u(0.5) = 5.59352829 \quad$$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

Hence the error we get for $$ n = 2 $$ is ,

<span id="(1)">
 * {| style="width:100%" border="0"

$$u^h(0.5) - u(0.5) = 8.74117 \times 10^{-3} \quad$$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 5.2.34)
 * }

MATLAB code with output
The K,d and F matrix are as follow:

K =

1        0         0         0         0         0         0         0         0         0    3.5000   -3.0000    2.7500   -2.6000    2.5000   -2.4286    2.3750   -2.3333         0   -3.0000    3.6667   -3.9000    4.0000   -4.0476    4.0714   -4.0833    4.0889         0    2.7500   -3.9000    4.5000   -4.8571    5.0893   -5.2500    5.3667   -5.4545         0   -2.6000    4.0000   -4.8571    5.4286   -5.8333    6.1333   -6.3636    6.5455         0    2.5000   -4.0476    5.0893   -5.8333    6.3889   -6.8182    7.1591   -7.4359         0   -2.4286    4.0714   -5.2500    6.1333   -6.8182    7.3636   -7.8077    8.1758         0    2.3750   -4.0833    5.3667   -6.3636    7.1591   -7.8077    8.3462   -8.8000         0   -2.3333    4.0889   -5.4545    6.5455   -7.4359    8.1758   -8.8000    9.3333

F =

[4 -12.8333   12.4167  -12.2500   12.1667  -12.1190   12.0893  -12.0694   12.0556]^T

d =

[4.0000  -2.8999    0.3731   -0.2799    0.3240    0.4387    0.8107    0.6131    0.2299]^T



The outputs are as follow:

@x=0.5 Exact U=5.593524e+000

U(weak form) for n=2 = 5.602355e+000

Error at n=2 is 8.831243e-003

@x=0.5 Exact U=5.593524e+000

U(weak form) for n=3 = 5.585545e+000

Error at n=3 is 7.978831e-003

@x=0.5 Exact U=5.593524e+000

U(weak form) for n=4 = 5.593383e+000

Error at n=4 is 1.409517e-004

@x=0.5 Exact U=5.593524e+000

U(weak form) for n=5 = 5.593788e+000

Error at n=5 is 2.637038e-004

@x=0.5 Exact U=5.593524e+000

U(weak form) for n=6 = 5.593527e+000

Error at n=6 is 3.409184e-006

@x=0.5 Exact U=5.593524e+000

U(weak form) for n=7 = 5.593514e+000

Error at n=7 is 9.890659e-006

@x=0.5 Exact U=5.593524e+000

U(weak form) for n=8 = 5.593524e+000

Error at n=8 is 1.058885e-007



u,u^h vs x Plot




Finding of the Approximate Solution using Fourier Basis Function satisfying Constraint Breaking Solution in Weak form
For $$ n = 1 \quad $$ www.wolframalpha.com

www.wolframalpha.com

www.wolframalpha.com

www.wolframalpha.com

www.wolframalpha.com

MATLAB code with output
The K,d and F matrix are as follows:

K =

1        0         0         0         0         0         0         0         0         0    0.7643    1.1171    2.3046    0.8846    3.0596   -0.7732    2.1100   -2.5560         0    1.1171    2.7357    3.6746    3.8194    5.9293    2.8546    6.5625    0.9196         0    2.3046    3.6746    7.1367    3.4374   10.0789   -1.1801    8.2449   -6.7905         0    0.8846    3.8194    3.4374    6.8633    7.1898    8.2933   11.0829    7.5264         0    3.0596    5.9293   10.0789    7.1898   16.1542    2.4145   17.1555   -5.6179         0   -0.7732    2.8546   -1.1801    8.2933    2.4145   15.3458   11.1219   19.5889         0    2.1100    6.5625    8.2449   11.0829   17.1555   11.1219   25.5917    4.9200         0   -2.5560    0.9196   -6.7905    7.5264   -5.6179   19.5889    4.9200   30.4083

F =

[4  -5.7179  -10.8903  -17.7236  -12.2749  -25.2744   -3.2817  -21.8270    7.5951]^T

d =

[4.0000 -62.5673   12.4807   25.5617  -15.3711   -4.9628    6.3734    0.2645   -0.9396]^T

The outputs are as follows:

@x=0.5 Exact U=5.593524e+000

U(weak form) for n=1 = 5.605964e+000

Error at n=1 is 1.244042e-002

@x=0.5 Exact U=5.593524e+000

U(weak form) for n=2 = 5.593171e+000

Error at n=2 is 3.531293e-004

@x=0.5 Exact U=5.593524e+000

U(weak form) for n=3 = 5.593534e+000

Error at n=3 is 9.958579e-006

@x=0.5 Exact U=5.593524e+000

U(weak form) for n=4 = 5.593523e+000

Error at n=4 is 3.393677e-007

Finding of the Approximate Solution using Exponential Basis Function satisfying Constraint Breaking Solution in Weak form
For $$ n = 2 \quad $$ www.wolframalpha.com

www.wolframalpha.com

www.wolframalpha.com

www.wolframalpha.com

www.wolframalpha.com

MATLAB code with output
The K,d and F matrix are as follows:

K =

1        0         0         0         0         0         0    1.7162    2.6335    3.1703    3.5125    3.7469         0    2.6335    4.2271    5.2687    5.9950    6.5286         0    3.1703    5.2687    6.7444    7.8343    8.6709         0    3.5125    5.9950    7.8343    9.2489   10.3699         0    3.7469    6.5286    8.6709   10.3699   11.7498

F =

[4  -8.2460  -11.4568  -12.7638  -13.3370  -13.6178]^T

d =

[4.0000 -101.6013 227.7000 -274.9216  169.6985  -42.1655]^T

The outputs areas follows:

@x=0.5 Exact U=5.593524e+000

U(weak form) for n=2 = 5.708257e+000

Error at n=2 is 1.147334e-001

@x=0.5 Exact U=5.593524e+000

U(weak form) for n=3 = 5.577020e+000

Error at n=3 is 1.650408e-002

@x=0.5 Exact U=5.593524e+000

U(weak form) for n=4 = 5.580883e+000

Error at n=4 is 1.264123e-002

@x=0.5 Exact U=5.593524e+000

U(weak form) for n=5 = 5.593442e+000

Error at n=5 is 8.216957e-005

Convergence Plot


= 5.3 =

Given: The Partial Differential Equation of G1MD1.0/D1b with Boundary Conditions
<span id="(1)">
 * {| style="width:100%" border="0"

Using the Lagrange Interpolant Basis Function (LIBF) with uniform discretization (equidistant nodes) where m = 4,6,8..., solve the following General 1-Dimensional Model 1.0/ Data Set 1.0b where $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$\frac{d }\left[ {{a_2}\left( x \right)\frac} \right] + f\left( {x,t} \right) = 0 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle (Eq. 5.3.1)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

and where $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$${a_2}\left( x \right){\text{ =  2 + 3x}} f\left( {x,t} \right) = 5x $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ u\left( {0} \right) = 4 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle (Eq. 5.3.2)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ -2\frac\left( 1 \right) = 12 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle (Eq. 5.3.3)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

with the essential boundary condition at 0 and the natural boundary condition at 1, as specified in, until convergence of $$ u_{m}^{h}(0.5) $$ is within 0(10^-6). $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

Find: Approximate Solutions using Lagrange Interpolant Basis Functions (LIBF)
<span id="(1)">
 * {| style="width:100%" border="0"

(1) Explain how the LIBF's are used as a "Constraint Breaking Solution" (CBS) $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

(2) Plot all LIBF's used. $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

(3) Using MATLAB's "quad" function, Wolframalpha, etc, integrate the LIBF's. $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

(4) Plot $$ u_m^h\ vs\ u$$ and $$ u_m^h\left(0.5 \right)-u\left(0.5 \right)\ vs\ m$$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

Solution: Deriving the Approximate Solution for Convergence with the Exact Form
<span id="(1)">
 * {| style="width:100%" border="0"

For the PDE and boundary conditions given above in Equation 5.3.1 through 5.3.3, the exact solution u(x) has been previously derived, illustrated in Equations 5.1.8 through 5.1.13. The exact solution itself may be seen as Equation 5.1.14. $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

In finding an approximate solution by using the LIBF, a series of Lagrangian equations are constructed for each node required; where the number of equations within a series equals the value of the particular node. $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

As the equation of the LIBF is seen below, $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ {L_{i,m}} = \prod\limits_{k = 1,k \ne i}^m {\frac} $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle (Eq. 5.3.4)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

a four (4) node system (m=4) would therefore be constructed of four (4) Lagrangian equations, as the examples below illustrate. $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ {L_{1,4}} = \frac $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle (Eq. 5.3.5)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ {L_{2,4}} = \frac $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle (Eq. 5.3.6)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ {L_{3,4}} = \frac $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle (Eq. 5.3.7)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ {L_{4,4}} = \frac $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle (Eq. 5.3.8)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

A six (6) node system would yield similar results; in that case, six Lagrangian equations. $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

So in order to find the displacement or approximate equation to compare to the exact solution, Equations 5.1.16 and 5.1.17 must be implemented. This requires the integration of the derived Langrangian basis functions, where from our example-equations seen above as Equations 5.3.5 through 5.3.8, $$ b_{1}^{'} $$ is defined to be $$ L_{1,4} $$, $$ b_{2}^{'} $$ is defined to be $$ L_{2,4} $$, and so on for m=4. Performing the necessary integrations will result in the tabulation of the stiffness $$ \tilde{K_{ij}} $$ matrix as well as the force matrix, $$ \tilde{F_i} $$. It should be noted that in implementing Equation 5.1.17, h is the natural boundary condition, evaluated at $$ \alpha =x $$. f equates to (5x), given by G1MD1.0/D1b $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

In taking the first row of the evaluated stiffness matrix $$ \tilde{K_{ij}} $$, and evaluating each term by the essential boundary condition, and by setting the first term of the force matrix $$ \tilde{F_i} $$ equal to 4, each value of d (approximate solution terms) may be solved. This procedure is similar to Equations 5.1.30 and 5.1.31 seen from Problem 5.1 above. $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

Lastly, similar to the results seen in Equation 5.1.32, the derived displacement or approximate equation may be derived from the d values. $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

With the approximated solution derived, the error that exists between this solution and the exact solution (for each node value) may be derived by evaluating each equation at 0.5 as specified in the problem statement and noting the difference. This procedure is further illustrated below. $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

The above procedure described above to derive the displacement or approximate solution may be repeated for each node value in ascending order, until an acceptable error tolerance is reached. For this problem, the MATLAB script seen below runs in iteration until the error tolerance reaches 0(10^-6) as specified and consequently, a particular node value will yield the desired results; the complete convergence of the approximate solution with the exact solution. $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

The following MATLAB script illustrates such iteration until convergence occures. The node value escalates, starting at 4 and ending at 10, or m=10. Each node value within the MATLAB script generates its own series of Lagrangian functions, each in-turn being integrated to solve for $$ \tilde{K_{ij}} $$ and $$ \tilde{F_i} $$, and therefore ultimately more accurate approximate solutions with each iteration and increasing node value. $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

MATLAB code with output
The K, d and F matrix are as follows:

K = 1.0e+003 *

0.0010        0         0         0         0         0        0         0         0         0   -0.0604    0.2192   -0.4084    0.5959   -0.6994    0.6118   -0.3952    0.1893   -0.0661    0.0131    0.0842   -0.4084    0.9467   -1.5042    1.8445   -1.7338    1.2243   -0.6545    0.2570   -0.0558   -0.1254    0.5959   -1.5042    2.6761   -3.5430    3.5444   -2.6918    1.5524   -0.6524    0.1479    0.1419   -0.6994    1.8445   -3.5430    5.1007   -5.4289    4.3217   -2.6144    1.1416   -0.2647   -0.1177    0.6118   -1.7338    3.5444   -5.4289    6.1450   -5.1089    3.1731   -1.4183    0.3332    0.0705   -0.3952    1.2243   -2.6918    4.3217   -5.1089    4.4584   -2.8582    1.2823   -0.3031   -0.0301    0.1893   -0.6545    1.5524   -2.6144    3.1731   -2.8582    1.9428   -0.9079    0.2074    0.0090   -0.0661    0.2570   -0.6524    1.1416   -1.4183    1.2823   -0.9079    0.5021   -0.1473   -0.0015    0.0131   -0.0558    0.1479   -0.2647    0.3332   -0.3031    0.2074   -0.1473    0.0707

F =

[4.0000 0.0318  0.2767 -0.2545  1.0648 -0.7424  1.3340  -0.2164  0.8466 12.1521]^T

d =

[4.0000 4.7446  5.3868  5.9485  6.4444  6.8852  7.2787  7.6308  7.9461  8.2283]^T

The outputs are as follows:

@x=0.5 Exact U=6.671156e+000

U(weak form) for n=4 = 6.678544e+000

Error at n=4 is 7.388478e-003

@x=0.5 Exact U=6.671156e+000

U(weak form) for n=6 = 6.670884e+000

Error at n=6 is 2.712999e-004

@x=0.5 Exact U=6.671156e+000

U(weak form) for n=8 = 6.671166e+000

Error at n=8 is 1.009478e-005

@x=0.5 Exact U=6.671156e+000

U(weak form) for n=10 = 6.671155e+000

Error at n=10 is 4.060753e-007

Solution to (1): LIBF used as a CBS
A Constraint Breaking Solution (CBS) utilizes an appropriate basis function, $$ {b_{i}(x), i= 1,2,3,...n} \quad $$ to eliminate constraints such that

$$ b_{0}(\beta)\neq 0 $$ and $$ b_{i}(\beta)=0 \quad\text{for}\quad \; i=1,2,3,...,n \quad $$ where $$ \beta \quad $$ is the location where the essential boundary condition is held.

To see how LIBF is used, consider Equation 5.3.4 and a 4 node system (m=4) with Equations 5.3.5 through 5.3.8.

<P> The value at all nodes will be zero except at one through use of the boundary conidition. Use $$ {L_{1,4}} \quad $$ to illustrate this concept. $$ b_0 ={L_{1,4}} \quad $$ and $$ {L_{1,4}}(x=0)=1 \quad $$, and $$( {L_{2,4}}, {L_{3,4}}, {L_{4,4}})=0 \quad $$ at x=0. This, demonstrates how the constraints were broken and thus how the LIBF can be used as a CBS.

Solution to (2): Plots of all LIBF's Used
After generating the aforementioned MATLAB script, 4 iterations (where ultimately m=10) are required before the approximated solution converges with the exact solution within the error tolerance 0(10-6).









Solution to (3): Using MATLAB to Integrate
As seen in the above MATLAB script, the MATLAB function "int" is used for integration to obtain both $$ \tilde{K_{ij}} \quad $$ and $$  \tilde{F_i} \quad $$ from equations 5.1.16 and 5.1.17

Solution to (4): Convergence and Error Plot
After running 4 iterations, when m=10, and the error tolerance is met; the convergence of $$  {u^h}_m \quad $$ (the approximated solution) with u (the exact solution) is illustrated in the below figure, with both equations being plotted on the same figure versus x.





The above Figure 5.3.6, illustrates the error deviation occurring between the approximate solution and exact solution, both evaluated at x=0.5, decreasing with respect to the increasing node value. As can be seen, the error tolerance is met when m=10.

--Eml5526.s11.team5.smith 02:58, 24 March 2011 (UTC)

--Eml5526.s11.team5.mcdaniel 03:02, 24 March 2011 (UTC)

= 5.4 =

Given: The Partial Differential Equation of G1MD1.0/D1 with Boundary Conditions
<span id="(1)">
 * {| style="width:100%" border="0"

Using the Lagrange Interpolant Basis Function (LIBF) with uniform discretization (equidistant nodes) where m = 4,6,8..., solve the following General 1-Dimensional Model 1.0/ Data Set 1.0 where $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$\frac{d }\left[ {{a_2}\left( x \right)\frac} \right] + f\left( {x,t} \right) = 0 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle (Eq. 5.4.1)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

and where $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$${a_2}\left( x \right){\text{ =  2 + 3x}} f\left( {x,t} \right) = 5x $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ u\left( {1} \right) = 4 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle (Eq. 5.4.2)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ -2\frac\left( 0 \right) = 12 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle (Eq. 5.4.3)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

with the essential boundary condition at 1 and the natural boundary condition at 0, as specified in, until convergence of $$ u_{m}^{h}(0.5) $$ is within 0(10^-6). $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

Find: Approximate Solutions using Lagrange Interpolant Basis Functions (LIBF)
<span id="(1)">
 * {| style="width:100%" border="0"

(1) Explain how the LIBF's are used as a "Constraint Breaking Solution" (CBS) $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

(2) Plot all LIBF's used. $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

(3) Using MATLAB's "quad" function, Wolframalpha, etc, integrate the LIBF's. $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

(4) Plot $$ u_m^h\ vs\ u$$ and $$ u_m^h\left(0.5 \right)-u\left(0.5 \right)\ vs\ m$$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

Solution: Deriving the Approximate Solution for Convergence with the Exact Form
<span id="(1)">
 * {| style="width:100%" border="0"

For the PDE and boundary conditions given above in Equation 5.4.1 through 5.4.3, the exact solution u(x) has been previously derived, illustrated in Equations 5.2.8 through 5.2.13. The exact solution itself may be seen as Equation 5.2.14. $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

In finding an approximate solution by using the LIBF, the procedure explained in the "Solution" section of problem 5.3 above is utilized once again. $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

The MATLAB script below again illustrates the convergence that occurs. The node value escalates as before, starting at 4 and ending at 10, or m=10. More approximate solutions with each iteration are derived. $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

MATLAB code with output
The K, d and F matrix are as follows:

K = 1.0e+003 *

0.0010        0         0         0         0         0         0         0        0         0   -0.1473    0.5021   -0.9079    1.2823   -1.4183    1.1416   -0.6524    0.2570   -0.0661    0.0090    0.2074   -0.9079    1.9428   -2.8582    3.1731   -2.6144    1.5524   -0.6545    0.1893   -0.0301   -0.3031    1.2823   -2.8582    4.4584   -5.1089    4.3217   -2.6918    1.2243   -0.3952    0.0705    0.3332   -1.4183    3.1731   -5.1089    6.1450   -5.4289    3.5444   -1.7338    0.6118   -0.1177   -0.2647    1.1416   -2.6144    4.3217   -5.4289    5.1007   -3.5430    1.8445   -0.6994    0.1419    0.1479   -0.6524    1.5524   -2.6918    3.5444   -3.5430    2.6761   -1.5042    0.5959   -0.1254   -0.0558    0.2570   -0.6545    1.2243   -1.7338    1.8445   -1.5042    0.9467   -0.4084    0.0842    0.0131   -0.0661    0.1893   -0.3952    0.6118   -0.6994    0.5959   -0.4084    0.2192   -0.0604   -0.0015    0.0090   -0.0301    0.0705   -0.1177    0.1419   -0.1254    0.0842   -0.0604    0.0295

F =

[4.0000 0.8466 -0.2164  1.3340 -0.7424  1.0648 -0.2545  0.2767  0.0318 12.0073]^T

d =

[4.0000 4.3272  4.6666  5.0215  5.3966  5.7978  6.2325  6.7112  7.2485  7.8656]^T

The outputs are as follows:

@x=0.5 Exact U=5.593524e+000

U(weak form) for n=4 = 5.585545e+000

Error at n=4 is 7.978810e-003

@x=0.5 Exact U=5.593524e+000

U(weak form) for n=6 = 5.593788e+000

Error at n=6 is 2.637247e-004

@x=0.5 Exact U=5.593524e+000

U(weak form) for n=8 = 5.593514e+000

Error at n=8 is 9.869443e-006

@x=0.5 Exact U=5.593524e+000

U(weak form) for n=10 = 5.593524e+000

Error at n=10 is 4.186644e-007

Solution to (1): LIBF used as a CBS
A Constraint Breaking Solution (CBS) utilizes an appropriate basis function, $$ {b_{i}(x), i= 1,2,3,...n} \quad $$ to eliminate contraints such that

$$ b_{0}(\beta)\neq 0 $$ and $$ b_{i}(\beta)=0 \quad\text{for}\quad \; i=1,2,3,...,n \quad $$ where $$ \beta \quad $$ is the location where the essential boundary condition is held.

To see how LIBF is used, consider Equation 5.3.4 and a 4 node system (m=4) with Equations 5.3.5 through 5.3.8.

<P> The value at all nodes will be zero except at one through use of the boundary condition. Use $$ {L_{4,4}} \quad $$ to illustrate this concept. $$ b_0 ={L_{4,4}} \quad $$ and $$ {L_{4,4}}(x=1)=1 \quad $$, and $$ ({L_{1,4}}, {L_{2,4}}, {L_{3,4}})=0 \quad $$ at x=1. This, demonstrates how the constraints were broken and thus how the LIBF can be used as a CBS.

Solution to (2): Plots of all LIBF's Used
After generating the aforementioned MATLAB script, 4 iterations (where ultimately m=10) are required before the approximated solution converges with the exact solution within the error tolerance 0(10-6).









Solution to (3): Using MATLAB to Integrate
As seen in the above MATLAB script, the MATLAB function "int" is used for integration to obtain both $$ \tilde{K_{ij}} \quad $$ and $$  \tilde{F_i} \quad $$ from equations 5.1.16 and 5.1.17

Solution to (4): Convergence and Error Plot
After running 4 iterations, when m=10, and the error tolerance is met; the convergence of $$  {u^h}_m \quad $$ (the approximated solution) with u (the exact solution) is illustrated in the below figure, with both equations being plotted on the same figure versus x.





The above Figure 5.4.6, illustrates the error deviation occurring between the approximate solution and exact solution, both evaluated at x=0.5, decreasing with respect to the increasing node value. As can be seen, the error tolerance is met when m=10.

--Eml5526.s11.team5.smith 03:00, 24 March 2011 (UTC)

--Eml5526.s11.team5.mcdaniel 03:04, 24 March 2011 (UTC)

=5.5 Calculix CCX examples=

Given: open source Finite Element code (Calulix)
<span id="(1)">
 * {| style="width:100%" border="0"

Calculix is an open source Finite Element code (with ABAQUS like input) available at http://www.dhondt.de/
 * style="width:95%" |
 * style="width:95%" |
 * }

Find: installation and run Calculix Crunchix(CCX)
<span id="(1)">
 * {| style="width:100%" border="0"

1) For the Disk problem, find
 * style="width:95%" |
 * style="width:95%" |
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

a)Node info
 * style="width:95%" |
 * style="width:95%" |
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

b)Element info
 * style="width:95%" |
 * style="width:95%" |
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

2) Generate 3 meshes of same disk with triangular elements. (increase the number of elements)
 * style="width:95%" |
 * style="width:95%" |
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

3) Install CCX, run examples and write a report explaining the commands used.
 * style="width:95%" |
 * style="width:95%" |
 * }

1) Node, Element Info
[| CGX Tutorial]

<span id="(1)">
 * {| style="width:100%" border="0"

The Node Info and Element Info can be extracted by including a command SEND ALL ABQ, after the mesh has been created.
 * style="width:95%" |
 * style="width:95%" |
 * }



<span id="(1)">
 * {| style="width:100%" border="0"

The node info, element info is shown below:
 * style="width:95%" |
 * style="width:95%" |
 * }

2) Triangular mesh
<span id="(1)">
 * {| style="width:100%" border="0"

Triangular mesh can be generated if we set the element type to tr3 or tr6.
 * style="width:95%" |
 * style="width:95%" |
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

Example: 'ELTY all tr3' or ' ELTY all tr6'
 * style="width:95%" |
 * style="width:95%" |
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

We have used element type 'tr3' to generate the triangular mesh.
 * style="width:95%" |
 * style="width:95%" |
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

The mesh size, i.e. the number of elements can be increased by increasing the number of divisions on the lines used to generate the geometry.
 * style="width:95%" |
 * style="width:95%" |
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

The code used to generate the meshes, the meshed geometries are shown below:
 * style="width:95%" |
 * style="width:95%" |
 * }

CCX Installation
<span id="(1)">
 * {| style="width:100%" border="0"

Windows executable downloaded from [| this link] comes along with CCX modules. For more details on installation please refer:[| CCX Installation Tutorial]
 * style="width:95%" |
 * style="width:95%" |
 * }

CCX Examples
<span id="(1)">
 * {| style="width:100%" border="0"

Basic examples demonstrating various features in CCX can be downloaded from:[| CCX Examples]
 * style="width:95%" |
 * style="width:95%" |
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

We were able to successfully run the below examples:
 * style="width:95%" |
 * style="width:95%" |
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

After the code is run, we need to post process the beam for displacements, forces. To post process the results, click on 'Tools>PostProcess' or hit 'SHIFT+F10'. The below screen shots demonstrate how the displacements and forces can be plotted.
 * style="width:95%" |
 * style="width:95%" |
 * }

<span id="(1)">
 * {| style="width:100%" border="0"


 * style="width:95%" |
 * Displacement Plot:
 * Displacement Plot:


 * }







<span id="(1)">
 * {| style="width:100%" border="0"


 * style="width:95%" |
 * Force Plot Plot:
 * Force Plot Plot:


 * }







<span id="(1)">
 * {| style="width:100%" border="0"


 * style="width:95%" |
 * Plate Mesh
 * Plate Mesh


 * }



<span id="(1)">
 * {| style="width:100%" border="0"


 * style="width:95%" |
 * Stress Plot
 * Stress Plot


 * }



<span id="(1)">
 * {| style="width:100%" border="0"


 * style="width:95%" |
 * Displacement Plot
 * Displacement Plot


 * }



<span id="(1)">
 * {| style="width:100%" border="0"


 * style="width:95%" |
 * Meshed Geometry
 * Meshed Geometry


 * }



<span id="(1)">
 * {| style="width:100%" border="0"


 * style="width:95%" |
 * Stress Plot
 * Stress Plot


 * }



<span id="(1)">
 * {| style="width:100%" border="0"


 * style="width:95%" |
 * Displacement Plot
 * Displacement Plot


 * }



CCX Report
<span id="(1)">
 * {| style="width:100%" border="0"

Commands Used to generate the .inp file are explained below:
 * style="width:95%" |
 * style="width:95%" |
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

Note: Few commands used for preprocessing were explained earlier at: [Preprocessing commands explained]
 * style="width:95%" |
 * style="width:95%" |
 * }

NSET
Keyword type: model definition

This option is used to assign nodes to a node set. The parameter NSET containing the name of the set is required (maximum 80 characters), whereas the parameter GENERATE (without value) is optional. If present, nodal ranges can be expressed by their initial value, their final value, and an increment. If a set with the same name already exists, it is reopened and complemented. The name of a set is case insensitive. Internally, it is modified into upper case and a 'N' is appended to denote it as node set.

Example:


 * NSET,NSET=N1

1,8,831,208


 * NSET,NSET=N2

100,N1

assigns the nodes with number 1, 8, 831 and 208 to (node) set N1 and the nodes with numbers 1, 8, 831, 208 (= set N1) and 100 to set N2

BOUNDARY
Keyword type: step or model definition

This option is used to prescribe boundary conditions. This includes:

temperature, displacements and rotations for structures total temperature, mass flow and total pressure for gas networks temperature, mass flow and static pressure for liquid networks temperature, mass flow and fluid depth for channels static temperature, velocity and static pressure for 3D-fluids. For liquids and structures the total and static temperature virtually coincide, therefore both are represented by the term temperature.

The following degrees of freedom are being used:

for structures:

1: translation in the local x-direction

2: translation in the local y-direction

3: translation in the local z-direction

4: rotation about the local x-axis

5: rotation about the local y-axis

6: rotation about the local z-axis

11: temperature

for gas networks:

1: mass flow

2: total pressure

11: total temperature

for liquid networks:

1: mass flow

2: static pressure

11: temperature

for liquid channels:

1: mass flow

2: fluid depth

11: temperature

for 3D-fluids:

1: velocity in the local x-direction

2: velocity in the local y-direction

3: velocity in the local z-direction

8: static pressure

11: static temperature

ELEMENT
Keyword type: model definition

With this option elements are defined. There is one required parameter, TYPE and one optional parameter, ELSET. The parameter TYPE defines the kind of element which is being defined. The following types can be selected:

General 3D solids

C3D4 (4-node linear tetrahedral element)

C3D6 (6-node linear triangular prism element)

C3D8 (3D 8-node linear isoparametric element)

C3D8R (the C3D8 element with reduced integration)

C3D10 (10-node quadratic tetrahedral element)

C3D15 (15-node quadratic triangular prism element)

C3D20 (3D 20-node quadratic isoparametric element)

C3D20R (the C3D20 element with reduced integration)

ABAQUS 3D solids for heat transfer (names are provided for compatibility)

DC3D4: identical to C3D4

DC3D6: identical to C3D6

DC3D8: identical to C3D8

DC3D10: identical to C3D10

DC3D15: identical to C3D15

DC3D20: identical to C3D20

Shell elements

S6 (6-node triangular shell element)

S8 (8-node quadratic shell element)

S8R (the S8 element with reduced integration)

Plane stress elements

CPS6 (6-node triangular plane stress element)

CPS8 (8-node quadratic plane stress element)

CPS8R (the CPS8 element with reduced integration)

Plane strain elements

CPE6 (6-node triangular plane strain element)

CPE8 (8-node quadratic plane strain element)

CPE8R (the CPS8 element with reduced integration)

Axisymmetric elements

CAX6 (6-node triangular axisymmetric element)

CAX8 (8-node quadratic axisymmetric element)

CAX8R (the CAX8 element with reduced integration)

Beam elements

B32 (3-node beam element)

B32R (the B32 element with reduced integration)

Special elements

D (3-node fluid element)

GAPUNI (2-node unidirectional gap element)

Example:


 * ELEMENT,ELSET=Eall,TYPE=C3D20R

1,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,

16,17,18,19,20

defines one 20-node element with reduced integration and stores it in set Eall.

EQUATION
Keyword type: model definition

With this option, a linear equation constraint between arbitrary displacement components at any nodes where these components are active can be imposed. The equation is assumed to be homogeneous, and all variables are to be written on the left hand side of the equation. The first variable is considered to be the dependent one, and is subsequently eliminated from the equation, i.e. the corresponding degree of freedom does not show up in the stiffness matrix. This reduces the size of the matrix. A node can only be used once as the dependent node in an equation or in a SPC.

First line:

Following lines, in a set: First line of set:
 * EQUATION

Number of terms in the equation. Following lines of set (maximum 12 entries per line):

Node number of the first variable.

Degree of freedom at above node for the first variable.

Value of the coefficient of the first variable.

Node number of the second variable.

Degree of freedom at above node for the second variable.

Value of the coefficient of the second variable.

Etc..

Continue, giving node number, degree of freedom, value of the coefficient, etc. Repeat the above line as often as needed if there are more than four terms in the *EQUATION. Specify exactly four terms per line for each constraint, except for the last line which may have less than four terms.

Example:

3 3,2,2.3,28,1,4.05,17,1,-8.22
 * EQUATION

defines an equation of the form, where u, v and w are the displacement for degree of freedom one, two and three, respectively

MATERIAL
Keyword type: model definition

This option is used to indicate the start of a material definition. A material data block is defined by the options between a *MATERIAL line and either another *MATERIAL line or a keyword line that does not define material properties. All material options within a data block will be assumed to define the same material. If a property is defined more than once for a material, the last definition is used. There is one required parameter, NAME, defining the name of the material with which it can be referenced in element property options (e.g. *SOLID SECTION). The name can contain up to 80 characters.

Material data requests outside the defined ranges are extrapolated in a constant way and a warning is generated. Be aware that this occasionally occurs due to rounding errors.

First line:

Enter the NAME parameter and its value.
 * MATERIAL

Example:


 * MATERIAL,NAME=EL

starts a material block with name EL.

ELASTIC
Keyword type: model definition, material

This option is used to define the elastic properties of a material. There is one optional parameter TYPE. Default is TYPE=ISO, other values are TYPE=ORTHO and TYPE=ENGINEERING CONSTANTS for orthotropic materials and TYPE=ANISO for anisotropic materials. All constants may be temperature dependent. For orthotropic and fully anisotropic materials, the coefficients satisfy the equation:

$$ S_{IJ}=D_{IJKL}E_{KL} (Eq 149) $$

where is the second Piola-Kirchhoff stress and  is the Lagrange deformation tensor. For linear calculations, these reduce to the generic stress and strain tensors.

First line:


 * ELASTIC

Enter the TYPE parameter and its value, if needed

Following line for TYPE=ISO:

Young's modulus.

Poisson's ratio.

Temperature.

Repeat this line if needed to define complete temperature dependence.

Example:


 * ELASTIC,TYPE=ORTHO

500000.,157200.,400000.,157200.,157200.,300000.,126200.,126200., 126200.,294.

defines an orthotropic material for temperature T=294. Since the definition includes values for only one temperature, they are valid for all temperatures.

SOLID SECTION
Keyword type: model definition

This option is used to assign material properties to 3D, plane stress, plane strain and axisymmetric element sets. The parameters ELSET and MATERIAL are required, the parameter ORIENTATION is optional. The parameter ELSET defines the element set to which the material specified by the parameter MATERIAL applies. The parameter ORIENTATION allows to assign local axes to the element set. If activated, the material properties are applied to the local axis. This is only relevant for non isotropic material behavior. For plane stress and plane strain elements the thickness can be specified on the second line. Default is 1.

First line:


 * SOLID SECTION

Enter any needed parameters.

Second line (only relevant for plane stress, plane strain and axisymmetric elements; can be omitted for 3D elements):

thickness (plane stress and plane strain elements)

Example:


 * SOLID SECTION,MATERIAL=EL,ELSET=Eall,ORIENTATION=OR1

assigns material EL with orientation OR1 to all elements in (element) set Eall.

ELSET
Keyword type: model definition

This option is used to assign elements to an element set. The parameter ELSET containing the name of the set is required (maximum 80 characters), whereas the parameter GENERATE (without value) is optional. If present, element ranges can be expressed by their initial value, their final value, and an increment. If a set with the same name already exists, it is reopened and complemented. The name of a set is case insensitive. Internally, it is modified into upper case and a 'E' is appended to denote it as element set.

First line:


 * ELSET

Enter any needed parameters and their values.

Following line if the GENERATE parameter is omitted:

List of elements and/or sets of elements previously defined to be assigned to this element set (maximum 16 entries per line).

Example:


 * ELSET,ELSET=E1,GENERATE

20,25


 * ELSET,ELSET=E2

E1,50,51

assigns the elements with numbers 20, 21, 22, 23, 24 and 25 to element set E1 and the elements with numbers 20, 21, 22, 23, 24, 25 (= set E1), 50 and 51 to element set E2.

SURFACE
Keyword type: model definition

This option is used to define surfaces made up of nodes or surfaces made up of element faces. A mixture of nodes and element faces belonging to one and the same surface is not possible. There are two parameters: NAME and TYPE. The parameter NAME containing the name of the surface is required. The TYPE parameter takes the value NODE for nodal surfaces and ELEMENT for element face surfaces. Default is TYPE=ELEMENT

Example:


 * SURFACE,NAME=left,TYPE=NODE

part,

1,

8

assigns the nodes with number 1, and 8 and the nodes belonging to node set part to a surface with name left.

Example:


 * SURFACE,NAME=new

38,S6

assigns the face 6 of element 38 to a surface with name new

CONTACT PAIR
Keyword type: model definition

This option is used to express that two surfaces can make contact. There is one required parameter: INTERACTION, and two optional parameters: SMALL SLIDING and ADJUST. The dependent surface is called the slave surface, the independent surface is the master surface. Surfaces are defined using the *SURFACE. The dependent surface can be defined as a nodal surface (option TYPE=NODE on the *SURFACE keyword) or as an element face surface (default for the *SURFACE card), whereas the independent surface has to be defined as an element face surface.

The INTERACTION parameter takes the name of the surface interaction (keyword *SURFACE INTERACTION) which applies to the contact pair. The surface interaction defines the nature of the contact (hard versus soft contact..)

First line:


 * CONTACT PAIR

enter the required parameter INTERACTION and its parameter.

Following line:

Name of the slave surface (can be nodal or element face based).

Name of the master surface (must be based on element faces).

Example:


 * CONTACT PAIR,INTERACTION=IN1,ADJUST=0.01

dep,ind

defines a contact pair consisting of the surface dep as dependent surface and the element face surface ind als independent surface. The name of the surface interaction is IN1. All slave nodes for which the clearance is smaller than or equal to 0.01 will moved onto the master surface.

SURFACE INTERACTION
Keyword type: model definition

This option is used to indicate the start of a surface interaction definition. A surface interaction data block is defined by the options between a *SURFACE INTERACTION line and either another *SURFACE INTERACTION line or a keyword line that does not define surface interaction properties. All surface interaction options within a data block will be assumed to define the same surface interaction. If a property is defined more than once for a surface interaction, the last definition is used. There is one required parameter, NAME, defining the name of the surface interaction with which it can be referenced in surface interaction property options (e.g. *SURFACE BEHAVIOR). The name can contain up to 80 characters.

Surface interaction data requests outside the defined ranges are extrapolated in a constant way. Be aware that this occasionally occurs due to rounding errors.

First line:


 * SURFACE INTERACTION

Enter the NAME parameter and its value.

Example:


 * SURFACE INTERACTION,NAME=SI1

starts a material block with name SI1.

BEAM SECTION
Keyword type: model definition

This option is used to assign material properties to beam element sets. The parameters ELSET, MATERIAL and SECTION are required, the parameters ORIENTATION, OFFSET1 and OFFSET2 are optional. The parameter ELSET defines the shell element set to which the material specified by the parameter MATERIAL applies. The parameter ORIENTATION allows to assign local axes to the element set. If activated, the material properties are applied to the local axis. This is only relevant for non isotropic material behavior.

The parameter SECTION defines the cross section of the beam and can take the value RECT for a rectangular cross section and CIRC for an elliptical cross section. A rectangular cross section is defined by its thickness in two perpendicular directions, an elliptical cross section is defined by the length of its principal axes. These directions are defined by specifying direction 1 on the third line of the present keyword card

First line:


 * BEAM SECTION

Enter any needed parameters.

Second line:

thickness in 1-direction

thickness in 2-direction

Third line:

global x-coordinate of a unit vector in 1-direction (default:0)

global y-coordinate of a unit vector in 1-direction (default:0)

global z-coordinate of a unit vector in 1-direction (default:-1)

Example:


 * BEAM SECTION,MATERIAL=EL,ELSET=Eall,OFFSET1=-0.5,SECTION=RECT

3.,1.

1.,0.,0.

assigns material EL to all elements in (element) set Eall. The reference line is in 1-direction on the back surface, in 2-direction on the central surface. The thickness in 1-direction is 3 unit lengths, in 2-direction 1 unit length. The 1-direction is the global x-axis.

STEP
Keyword type: step

This card describes the start of a new STEP. PERTURBATION, NLGEOM, INC, INCF and TURBULENCE MODEL are the optional parameters.

First and only line:


 * STEP

Enter any needed parameters and their values

Example:


 * STEP,INC=1000,INCF=20000,TURBULENCE MODEL=SST

starts a step and increases the maximum number of thermomechanical increments to complete the step to 1000. The maximum number of 3D fluid increments is set to 20000 and for the turbulence model the SST model was chosen.

STATIC
Keyword type: step

This procedure is used to perform a static analysis. The load consists of the sum of the load of the last *STATIC step and the load specified in the present step with replacement of redefined loads.

First line:


 * STATIC

Enter any needed parameters and their values. Second line (only relevant for nonlinear analyses; for linear analyses, the step length is always 1) Initial time increment. This value will be modified due to automatic incrementation, unless the parameter DIRECT was specified (default 1.). Time period of the step (default 1.). Minimum time increment allowed. Only active if DIRECT is not specified. Default is the initial time increment or 1.e-5 times the time period of the step, whichever is smaller. Maximum time increment allowed. Only active if DIRECT is not specified. Default is 1.e+30.

Example:

.1,1.
 * STATIC,DIRECT

defines a static step and selects the SPOOLES solver as linear equation solver in the step (default). If the step is a linear one, the other parameters are of no importance. If the step is nonlinear, the second line indicates that the initial time increment is .1 and the total step time is 1. Furthermore, the parameter DIRECT leads to a fixed time increment. Thus, if successful, the calculation consists of 10 increments of length 0.1.

CLOAD
Keyword type: step

This option allows concentrated forces to be applied to any node in the model which is not fixed by a single or multiple point constraint. Optional parameters are OP, AMPLITUDE, TIME DELAY, LOAD CASE and SECTOR. OP can take the value NEW or MOD. OP=MOD is default and implies that the concentrated loads applied to different nodes are kept over all steps starting from the last perturbation step. Specifying a force in a node for which a force was defined in a previous step replaces this value. OP=NEW implies that all previously applied concentrated loads are removed. If multiple *CLOAD cards are present in a step this parameter takes effect for the first *CLOAD card only.

First line:


 * CLOAD

Enter any needed parameters and their value.

Following line:

Node number or node set label.

Degree of freedom.

Magnitude of the load

Repeat this line if needed.

Example:

1000,3,10.3
 * CLOAD,OP=NEW,AMPLITUDE=A1,TIME DELAY=20.

removes all previous point load forces and applies a force with magnitude 10.3 and amplitude A1 (shifted in positive time direction by 20 time units) for degree of freedom three (global if no transformation was defined for node 1000, else local) of node 1000.

NODE PRINT
Keyword type: step

This option is used to print selected nodal variables in file jobname.dat. The following variables can be selected:

Displacements (key=U)

Structural temperatures, total temperatures in networks and static temperatures in 3D fluids (key=NT)

Pressures in networks (key=PN). These are the total pressures for gases, static pressures for liquids and liquid depth for channels. The fluid section types dictate the kind of network.

Static pressures in 3D fluids (key=PS)

Velocities in 3D fluids (key=V)

Mass flows in networks (key=MF)

External forces (key=RF)

External concentrated heat sources (key=RFL)

First line:


 * NODE PRINT

Enter the parameter NSET and its value.

Second line:

Identifying keys for the variables to be printed, separated by kommas.

Example:


 * NODE PRINT,NSET=N1

RF

requests the storage of the reaction forces in the nodes belonging to (node) set N1 in the .dat file.

EL PRINT
Keyword type: step

This option is used to print selected element variables in an ASCII file with the name jobname.dat. Some of the element variables are printed in the integration points, some are whole element variables. The following variables can be selected:

Integration point variables

true (Cauchy) stress (key=S)

strain (key=E). This is the Lagrangian strain for (hyper)elastic materials and incremental plasticity and the Eulerian strain for deformation plasticity.

equivalent plastic strain (key=PEEQ)

equivalent creep strain (key=CEEQ; is converted internally into PEEQ since the viscoplastic theory does not distinguish between the two; consequently, the user will find PEEQ in the dat file, not CEEQ)

the energy density (key=ENER)

the internal state variables (key=SDV)

heat flux (key=HFL)

Whole element variables

the internal energy (key=ELSE)

the kinetic energy (key=ELKE)

the volume (key=EVOL)

First line:


 * EL PRINT

Enter the parameter ELSET and its value.

Second line:

Identifying keys for the variables to be printed, separated by kommas.

Example:


 * EL PRINT,ELSET=Copper

E

requests to store the strains at the integration points in the elements of set Copper in the .dat file.

NODE FILE
Keyword type: step

This option is used to print selected nodal variables in file jobname.frd for subsequent viewing by CalculiX GraphiX. The following variables can be selected:

Displacements (key=U)

Displacements: magnitude and phase (key=PU, only for *STEADY STATE DYNAMICS calculations) and *FREQUENCY calculations with cyclic symmetry).

Maximum displacements orthogonal to a given vector at all times for *FREQUENCY calculations with cyclic symmetry. The components of the vector are the coordinates of a node stored in a node set with the name RAY. This node and node set must have been defined by the user (key=MAXU).

Temperatures (key=NT). This includes both structural temperatures and total fluid temperatures in a network.

Temperatures: magnitude and phase (key=PNT, only for *STEADY STATE DYNAMICS calculations).

External forces (key=RF)

External concentrated heat sources (key=RFL)

Static temperatures in networks and 3D fluids (key=TS)

Total temperatures in networks and 3D fluids (key=TT)

Mass flows in networks (key=MF). The mass flow through a fluid element is stored in the middle node of the element. In the end nodes the mass flow is not unique, since more than two element can be connected to the node. For end nodes the sum of the mass flow leaving the node is stored. Notice that at nodes where mass flow is leaving the network the value will be wrong if no proper exit element (with node number 0) is attached to that node.

Velocities in 3D fluids (key=V)

Mach numbers in 3D compressible fluids (key=MACH)

Static pressures in liquid networks and 3D fluids (key=PS)

Total pressures in gas networks and 3D fluids (key=PT)

Fluid depth in channel networks (key=DEPT)

Critical depth in channel networks (key=HCRI)

Pressure coefficient in 3D compressible fluids (key=CP)

First line:


 * NODE FILE

Enter any needed parameters and their values.

Second line:

Identifying keys for the variables to be printed, separated by kommas.

Example:


 * NODE FILE,FREQUENCY=2,TIME POINTS=T1

RF,NT

requests the storage of reaction forces and temperatures in the .frd file every second increment. In addition, output will be stored for all time points defined by the T1 time points sequence.

EL FILE
Keyword type: step

This option is used to save selected element variables averaged at the nodal points in a frd file (extension .frd) for subsequent viewing by CalculiX GraphiX. The following element variables can be selected:

true (Cauchy) stress (key=S). For beam elements this tensor is replaced by the section forces if SECTION FORCES is selected.

stress: magnitude and phase (key=PHS, only for *STEADY STATE DYNAMICS calculations and *FREQUENCY calculations with cyclic symmetry). maximum worst principal stress at all times for *FREQUENCY calculations with cyclic symmetry. It is stored for nodes belonging to the node set with name STRESSDOMAIN. This node set must have been defined by the user with the *NSET command. The worst principal stress is the maximum of the absolute value of the principal stresses (key=MAXS).

strain (key=E). This is the Lagrangian strain for (hyper)elastic materials and incremental plasticity and the Eulerian strain for deformation plasticity.

equivalent plastic strain (key=PEEQ)

equivalent creep strain (key=CEEQ; is converted internally into PEEQ since the viscoplastic theory does not distinguish between the two; consequently, the user will find PEEQ in the frd file, not CEEQ)

the energy density (key=ENER)

the internal state variables (key=SDV)

heat flux (key=HFL)

First line:


 * EL FILE

Enter any needed parameters and their values.

Second line:

Identifying keys for the variables to be printed, separated by kommas.

Example:


 * EL FILE

S,PEEQ

requests that the (Cauchy) stresses and the equivalent plastic strain is stored in .frd format for subsequent viewing with CalculiX GraphiX.

END STEP
Keyword type: step

This option concludes the definition of a step.

First and only line:


 * END STEP

Example:


 * END STEP

concludes a step. Each *STEP card must at some point be followed by an *END STEP card.

--Eml5526.s11.team5.vijay 01:48, 24 March 2011 (UTC)

=5.6 Quad Lagrangian Element Basis Function=

Given: Quad Lagrangian Basis Function
Lagrangian element with 3 nodes per element.

Find
Plot the basis functions for each node

Solution
Lagrange basis functions for an element of order m can be written as follow:

$$ L_{i,m}(x)=\prod_{k=1}^{m} \frac{x-x_{k}}{x_{i}-x_{k}} \forall k\ne i$$

For a Quad lagrangian Basis function there are 3 nodes per element (m=3). We can write the Quad Lagrangian Basis Function as follows:

$$ L_{1,3}(x)= (\frac{x-x_{2}}{x_{1}-x_{2}})(\frac{x-x_{3}}{x_{1}-x_{3}}) $$

$$ L_{2,3}(x)= (\frac{x-x_{1}}{x_{2}-x_{1}})(\frac{x-x_{3}}{x_{2}-x_{3}}) $$

$$ L_{3,3}(x)= (\frac{x-x_{1}}{x_{3}-x_{1}})(\frac{x-x_{2}}{x_{3}-x_{2}}) $$

Let us consider the element shown below:



In the above element $$X_{1}=0, X_{2}=0.5, X_{3}=1 $$

Substituting the above values, the basis functions at each node are as follows:

--Eml5526.s11.team5.vijay 01:49, 24 March 2011 (UTC)

= 5.7 =

P.D.E
<span id="(1)">
 * {| style="width:100%" border="0"

$$ \frac{\partial}{\partial x} \left[ (2+3x) \frac{\partial u}{\partial x} \right] + 5x = 0 \ \ \forall x \in \left] 0,1 \right[ $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 5.7.1)
 * }

Boundary Conditions
<span id="(1)">
 * {| style="width:100%" border="0"

$$u \ at \ {x=0} \ is \ 4 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 5.7.2)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$  \frac{\partial u}{\partial x} \  at \ x = 1 \ is \ \frac{12}{5} $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 5.7.3)
 * }

Find: Explain how LLEBF is used as CBS, get approximate solution by using LLEBF and compare to exact solution

 * 1) Explain how Linear Lagrangian Element Basis Functions (LLEBF) are used as a Constraint Breaking Solution (CBS).


 * 2) Plot all LLEBF used for 4 number of elements (nel=4).


 * 3) Use Matlab to integrate and solve with nel=4,6,8,... until error at x=0.5 is less then 10^-6.


 * 4) Plot the exact vs the approximated solution and the error at 0.5 vs number of elements.

Solution
1)

<span id="(1)">
 * {| style="width:100%" border="0"

A function $$\left \{ b_i(x), i=0,1,2,... \right \} $$ is a CBS if $$ b_0(\Beta) \ne 0$$   and   $$ b_i(\Beta)=0 \quad$$  for  $$i=1,2,...,n \quad$$
 * style="width:95%" |
 * style="width:95%" |

In this problem we are using a linear Lagrange polynomial to interpolate $$u(x) \quad$$ :


 * $$u(x)= u(x_1)\frac{x-x_{2}}{x_{1} - x_{2}}+u(x_2) \frac{x-x_{1}}{x_{2} - x_{1}}\quad$$

As you can see when $$x=x_1 \quad$$ the value contributed from $$u(x_2)\quad$$  is 0 and the value contributed from $$u(x_1)\quad$$ is 1, which is not equal to 0. This makes it a CBS.
 * }

2) <span id="(1)">
 * {| style="width:100%" border="0"

The Plot of all LLBEF for element 1 is given in figure 5.7.1. Since each element has the same length and same number of nodes thier plot of LLBEF's will also be identical. The plot of all 4 elements together is seen in figure 5.7.2.
 * style="width:95%" |
 * style="width:95%" |


 * }

Figure 5.7.1. Single element representation by LLEBS function in discritized truss Figure 5.7.2.Global view of a truss with LLEBS bases shown for each element

3)

Given on 30-5

<span id="(1)">
 * {| style="width:100%" border="0"

$$ K_{ij}^e = \int_{w^e} {b^e_i}'(x)a_2(x){b^e_j}'(x) dx \ $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

, where $$ {b^e_i} \ $$ is linear Lagrange function

<span id="(1)">
 * {| style="width:100%" border="0"

$$ {b^e_i} = L_{i,n}= \frac{x-x_j}{x_i - x_j} \ $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

in this specific case of linear Lagrange interpolation

<span id="(1)">
 * {| style="width:100%" border="0"

$$ {b^e_i} = L_{1,2}= \frac{x-x_2}{x_1 - x_2} = \frac{x_2 - x}{x_2 - x_1}\ $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ {b^e_i} = L_{2,2}= \frac{x-x_1}{x_2 - x_1} \ $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

then by taking derivative

<span id="(1)">
 * {| style="width:100%" border="0"

$$ {b^e_1}' = L_{1,2}'= \frac{-1}{x_2 - x_1} = \frac{-1}{l_e} \ $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ {b^e_2}' = L_{2,2}'= \frac{1}{x_2 - x_1} = \frac{1}{l_e} \ $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

Then we construct Gram matrix

<span id="(1)">
 * {| style="width:100%" border="0"

$$ K_{ij}^e = \int_{w^e} {b^e_i}'(x)a_2(x){b^e_j}'(x) dx =  \int_{w^e} \Gamma a_2(x) dx \ $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ \Gamma =
 * style="width:95%" |
 * style="width:95%" |

\begin{bmatrix} b_1b_1     & b_1b_2     \\ b_2b_1    & b_2b_2 \end{bmatrix}

= \frac{1}{l_e^2}

\begin{bmatrix} 1    & -1    \\ 1     & -1 \end{bmatrix}

$$ $$
 * <p style="text-align:right"> $$ \displaystyle
 * }

remembering that $$ a2 = 2+3x \ $$, we get

<span id="(1)">
 * {| style="width:100%" border="0"

$$ \int_{w^e}
 * style="width:95%" |
 * style="width:95%" |

\frac{1}{l_e^2}

\begin{bmatrix} 1    & -1    \\ 1     & -1 \end{bmatrix}

2+3x dx \ $$ $$
 * <p style="text-align:right"> $$ \displaystyle
 * }

after integration we get some (2x2) matrix, and noting the element nodes, we construct the global matrix

For the force matrix we do the same. We assemble the local element vector "N"

<span id="(1)">
 * {| style="width:100%" border="0"

$$ N^e = [ b_i \ b_j] \ $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

and then we keep assembling the elements

<span id="(1)">
 * {| style="width:100%" border="0"

$$ F = \int_{w^e} N^{eT} f(x) dx \ $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

After we get stiffness and force matrices, we apply boundary conditions

<span id="(1)">
 * {| style="width:100%" border="0"

$$u \ at \ {x=0} \ is \ 4 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 5.7.2)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$  \frac{\partial u}{\partial x} \  at \ x = 1 \ is -6 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 5.7.3)
 * }

to get final stiffness matrix $$ K \ $$ and $$ F \ $$

We generated matlab code that computes stiffness matrix $$ K \ $$ and force matrix $$ F \ $$, puts necessary B.C.s and then plots the results

, where errorMatirix, gives number of elements that were used and corresponding error

In order to get to desired accuracy $$ O(10^{-6} \ $$, we subsequently increased the number of elements till computational time was no longer plausible. The number of elements that we stopped was 80. The solution was so close to the desired value that it was converging past 30 elements at an extremely slow rate, while computational time with each added element increased significantly. Below in drop-down table are shown numerical results of number of elements used and respective absolute error associated with elements at x = 0.5.

4)

Plots below show the error convergence towards 0 as well as convergence rate. Also note the power on "y" axes, and relative convergence of error passed 30 elements. It is really slow at that point.

Figure 5.7.3. Analytical solution of truss segmented into 80 elements vs. analytical solution

Figure 5.7.4.Absolute error change as a function of number of segmentation (elements) of truss in Data Set 1 (D1)

= 5.8 =

P.D.E
<span id="(1)">
 * {| style="width:100%" border="0"

$$ \frac{\partial}{\partial x} \left[ (2+3x) \frac{\partial u}{\partial x} \right] + 5x = 0 \ \ \forall x \in \left] 0,1 \right[ $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 5.8.1)
 * }

Boundary Conditions
<span id="(1)">
 * {| style="width:100%" border="0"

$$u \ at \ {x=1} \ is \ 4 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 5.8.2)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$  \frac{\partial u}{\partial x} \  at \ x = 0 \ is -6 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 5.8.3)
 * }

Find: Data set 1 D1 for generic truss

 * 1) Explain how Linear Lagrangian Element Basis Functions (LLEBF) are used as a Constraint Breaking Solution (CBS).


 * 2) Plot all LLEBF used for 4 number of elements (nel=4).


 * 3) Use Matlab to integrate and solve with nel=4,6,8,... until error at x=0.5 is less then 10^-6.


 * 4) Plot the exact vs the approximated solution and the error at 0.5 vs number of elements.

Solution
1)

<span id="(1)">
 * {| style="width:100%" border="0"

A function $$\left \{ b_i(x), i=0,1,2,... \right \} $$ is a CBS if $$ b_0(\Beta) \ne 0$$   and   $$ b_i(\Beta)=0 \quad$$  for  $$i=1,2,...,n \quad$$
 * style="width:95%" |
 * style="width:95%" |

In this problem we are using a linear Lagrange polynomial to interpolate $$u(x) \quad$$ :


 * $$u(x)= u(x_1)\frac{x-x_{2}}{x_{1} - x_{2}}+u(x_2) \frac{x-x_{1}}{x_{2} - x_{1}}\quad$$

As you can see when $$x=x_1 \quad$$ the value contributed from $$u(x_2)\quad$$  is 0 and the value contributed from $$u(x_1)\quad$$ is 1, which is not equal to 0. This makes it a CBS.
 * }

2) <span id="(1)">
 * {| style="width:100%" border="0"

The Plot of all LLBEF for element 1 is given in figure 5.7.1. Since each element has the same length and same number of nodes thier plot of LLBEF's will also be identical. The plot of all 4 elements together is seen in figure 5.7.2.
 * style="width:95%" |
 * style="width:95%" |


 * }

Single element representation by LLEBS function in discritized truss

Global view of a truss with LLEBS bases shown for each element

3)

Given on 30-5

<span id="(1)">
 * {| style="width:100%" border="0"

$$ K_{ij}^e = \int_{w^e} {b^e_i}'(x)a_2(x){b^e_j}'(x) dx \ $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

, where $$ {b^e_i} \ $$ is linear Lagrange function

<span id="(1)">
 * {| style="width:100%" border="0"

$$ {b^e_i} = L_{i,n}= \frac{x-x_j}{x_i - x_j} \ $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

in this specific case of linear Lagrange interpolation

<span id="(1)">
 * {| style="width:100%" border="0"

$$ {b^e_i} = L_{1,2}= \frac{x-x_2}{x_1 - x_2} = \frac{x_2 - x}{x_2 - x_1}\ $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ {b^e_i} = L_{2,2}= \frac{x-x_1}{x_2 - x_1} \ $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

then by taking derivative

<span id="(1)">
 * {| style="width:100%" border="0"

$$ {b^e_1}' = L_{1,2}'= \frac{-1}{x_2 - x_1} = \frac{-1}{l_e} \ $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ {b^e_2}' = L_{2,2}'= \frac{1}{x_2 - x_1} = \frac{1}{l_e} \ $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

Then we construct Gram matrix

<span id="(1)">
 * {| style="width:100%" border="0"

$$ K_{ij}^e = \int_{w^e} {b^e_i}'(x)a_2(x){b^e_j}'(x) dx =  \int_{w^e} \Gamma a_2(x) dx \ $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$ \Gamma =
 * style="width:95%" |
 * style="width:95%" |

\begin{bmatrix} b_1b_1     & b_1b_2     \\ b_2b_1    & b_2b_2 \end{bmatrix}

= \frac{1}{l_e^2}

\begin{bmatrix} 1    & -1    \\ 1     & -1 \end{bmatrix}

$$ $$
 * <p style="text-align:right"> $$ \displaystyle
 * }

remembering that $$ a2 = 2+3x \ $$, we get

<span id="(1)">
 * {| style="width:100%" border="0"

$$ \int_{w^e}
 * style="width:95%" |
 * style="width:95%" |

\frac{1}{l_e^2}

\begin{bmatrix} 1    & -1    \\ 1     & -1 \end{bmatrix}

2+3x dx \ $$ $$
 * <p style="text-align:right"> $$ \displaystyle
 * }

after integration we get some (2x2) matrix, and noting the element nodes, we construct the global matrix

For the force matrix we do the same. We assemble the local element vector "N"

<span id="(1)">
 * {| style="width:100%" border="0"

$$ N^e = [ b_i \ b_j] \ $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

and then we keep assembling the elements

<span id="(1)">
 * {| style="width:100%" border="0"

$$ F = \int_{w^e} N^{eT} f(x) dx \ $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle
 * }

After we get stiffness and force matrices, we apply boundary conditions

<span id="(1)">
 * {| style="width:100%" border="0"

$$u \ at \ {x=1} \ is \ 4 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 5.8.4)
 * }

<span id="(1)">
 * {| style="width:100%" border="0"

$$  \frac{\partial u}{\partial x} \  at \ x = 0 \ is -6 $$ $$
 * style="width:95%" |
 * style="width:95%" |
 * <p style="text-align:right"> $$ \displaystyle                                                                          (Eq. 5.8.5)
 * }

to get final stiffness matrix $$ K \ $$ and $$ F \ $$

We generated matlab code that computes stiffness matrix $$ K \ $$ and force matrix $$ F \ $$, puts necessary B.C.s and then plots the results

, where errorMatirix, gives number of elements that were used and corresponding error

In order to get to desired accuracy $$ O(10^{-6} \ $$, we subsequently increased the number of elements till computational time was no longer plausible. The number of elements that we stopped was 80. The solution was so close to the desired value that it was converging past 30 elements at an extremely slow rate, while computational time with each added element increased significantly. Below in drop-down table are shown numerical results of number of elements used and respective absolute error associated with elements at x = 0.5.

4)

Plots below show the error convergence towards 0 as well as convergence rate. Also note the power on "y" axes, and relative convergence of error passed 30 elements. It is really slow at that point.

Figure 5.7.3. Analytical solution of truss segmented into 80 elements vs. analytical solution

Figure 5.7.4.Absolute error change as a function of number of segmentation (elements) of truss in Data Set 1 (D1)

=Contributing Members =

= References =