User:Eml5526.s11.team6/hwk5

= Problem 5.1 Solution of the PDE using weak form.=

Strong Form

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

$$ \frac{\partial}{\partial x}\left [ (2+3x)\frac{\partial u}{\partial x} \right ]+5x=0 \qquad \forall  x   \epsilon \left ]0,1   \right [$$ (5.1.1)
 * }
 * }

Boundary Conditions

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

$$ u(0)=4, \qquad -\frac{\partial u}{\partial x}_{(x=1)}=6 $$
 *  (5.1.2)
 * }
 * }

Solve for $$ u(x) $$ and $$ u^h(x) $$
To solve the above problem for $$ u(x) $$ and $$ u^h(x) $$ using appropriate basis functions of the following types until convergence of u^h(0.5) -- 1. Fourier basis function 2. Polynomial basis function 3. Exponential basis function

Introducing the arbitrary weighted function
Weak form can be found as follows, with weighted function

Now, multiply eq. 5.1.1 by $$ \displaystyle w $$ and integrate between 0 to 1.

Use integration by parts for the first term in eq. 5.1.4

Using eq. (5.1.2) and eq. (5.1.3), we get

Put eq. (5.1.6) in eq. (5.1.5), we get the weak form as follows

Using the polynomial basis function
For getting appropriate constraint breaking solution, we choose Also, Take $$ \displaystyle \beta=0$$ for the constraint breaking solution. Therefore, From eq. 5.1.2 and eq. 5.1.10, we can say that Therefore, eq. 5.1.10 becomes Similarly, the weighted function can be written as Eq. 5.1.3 gives, Therefore eq. 5.1.13 becomes Hence, from eq. 5.2.12 and 5.2.15 Put this in eq. 5.1.7, we get Using eq. 5.1.15 we can find $$ \displaystyle w(1) $$ and put in eq. 5.1.18, we get By integrating and verifying with| Wolfram Alpha, we get The following figure shows comparison of analytical and numerical solution for different values of n. As expected, the numerical solution approaches analytical solution as n increases.

Approximate solution converges to analytical solution at n=4. To minimize the error,higher values of n need to be considered.

Using the exponential basis function
For getting appropriate constraint breaking solution, we choose Also, Take $$ \displaystyle \beta=0$$ for the constraint breaking solution. Therefore, From eq. 5.1.2 and eq. 5.1.25, we can say that Therefore, eq. 5.1.25 becomes Similarly, the weighted function can be written as Eq. 5.1.3 gives, Therefore eq. 5.1.28 becomes Hence, from eq. 5.2.27 and 5.2.30 Put this in eq. 5.1.7, we get Using eq. 5.1.30 we can find $$ \displaystyle w(1) $$ and put in eq. 5.1.33, we get By integrating and verifying with| Wolfram Alpha, we get

Following Matlab code is used to find coefficients a_j by changing values of i and j in the above equation After getting the coefficients of numerical solution, we compare numerical solution with analytical solution for different values of n.

Then graph of error that is difference between analytical solution and numerical solution is plotted at x=0.5



As expected, with increase in number of n,numerical solution approaches analytical solution.

From the above graph, error is almost zero for n=6 and above. But to get the error in the order of less than $$10^{-6}$$,we solve it further for higher values of n.

Using the Fourier basis function
The trial solution then becomes the following for n=1 to convergence:

The essential boundary condition,$$u(x=0)=4$$ must be satisfied, and the basis functions must be selected so $$ w(n=0) = 0 $$. Solving the above basis function for u(x=1)=4 yielded the following: $$ {d}_{1}=4, i=1,2,\cdots $$, $$ j=1,2,\cdots $$

The remaining displacement coefficients must be solved for $$d_{j}, \;\;j=1,2,\cdots $$, by constructing the following stiffness and force matrices using the switching Fourier basis function.

Matlab was used to generate the following stiffness and force matrices shown below.

Since Kd=F, d can be solved for using Matlab to perform the matrix algebra. The resulting d matrix is shown below.

The trial solution of the governing equation then becomes the following.

Now we plot for Exact solution u(x) and Approximate solution uh(x) and also Convergence of error vs. number of degrees of freedom (n)

As seen from the Figure 5.1.1 below, it is evident that the approximate solution uh(x) very closely imitates the exact solution u(x) Also it can be seen from Figure 5.1.2. that the error converges to almost zero value when more than 8 basis functions are used





=Problem 5.2 - Solutions using weak form and appropriate families of basis functions=

Strong Form
The strong form as referred from General one dimensional Model 1.0/Data set 1b on page 2 of Mtg 9:-

Boundary Conditions
The Essential Boundary condition as applied on the boundary $$\Gamma_g = [1]$$, is -- $$ u(1)=4 $$, The Natural Boundary condition as applied on the boundary $$\Gamma_h = [0]$$, is -- $$ -\frac{\partial u(x=0)}{\partial x}=6 $$

Solve for $$ u(x) $$ and $$ u^h(x) $$
To solve the above problem for $$ u(x) $$ and $$ u^h(x) $$ using appropriate basis functions of the following types until convergence of $$u^h(0.5)$$ to the order of $$10^{-6}$$ -- 1. Fourier basis function 2. Polynomial basis function 3. Exponential basis function

Solution
Using the technique of integration by parts, we have derived the weak form as follows --

We make following choice for the arbitrary weighted function for our convenience -

Now, multiply Eq. 5.2.2 by $$ \displaystyle w $$ and integrate between 0 to 1.

Use integration by parts for the first term in Eq. 5.2.3

Using Eq. (5.2.3) and Eq. (Eq. 5.2.4), we get

From Eq. (5.2.4) in Eq. (5.2.5), we get the weak form as follows -

Now, using the technique from Mtg 18 and Mtg 19 we can write the discretized weak form as -

We choose the basis functions such that $$ w(n=0) = 1 $$, to eliminate the unknown flux at this boundary. Using technique from page -1 Mtg 22, we can write - $$ {{d}_{0}}=(\frac {g}{b_0}) = (\frac {4}{1})= 4, $$ In order to plot the exact solution u(x), solving the above ordinary differential $$ Eq. 5.2.1 $$, we get the exact solution of governing equation as follows -

Now we select legitimate families of basis functions as follows -

Family of Fourier basis function
Now we can write the Approximated trial solution $$ u^h (x) $$ as follows -

Here we have -  [ j = 0,1,2,3,4...]

Now we construct the stiffness matrix as follows:

Note- Here 1st equation is valid when 'i' and 'j' both are even, 2nd expression is valid when 'i' and 'j' both are is odd, 3rd expression is valid for remaining combination of values of 'i' and 'j'

Similarly, the force matrix can be shown as below-

Note- Here 1st expression is valid when 'i' and 'j' both are even, 2nd expression is valid when 'i' and 'j' both are odd

The stiffness matrix and Force matrix are used to solve for the remaining coefficients $${{d}_{j}},_ – ^ – j=1,2,\cdots $$

Also now we need to have convergence of the solutions to the order of 10^(-6), hence we use a Matlab code demonstrated below to build the above matrices which are obtained then as follows-

By the equation that$$Kd=F$$, we can solve the remaining of $$d$$ as

Hence, the trial or approximated solution of the governing equation from $$ Eq. 5.2.1 $$ is-

Now we plot for Exact solution $$ u (x)$$ and Approximate solution $$ u^h (x) $$ and also Convergence of error vs. number of degrees of freedom (n)

As seen from the Figure 5.2.1 below, it is evident that the approximate solution $$ u^h(x) $$ very closely imitates the exact solution $$ u(x) $$ Also it can be seen from Figure 5.2.2. that the error very rapidly converges to almost zero value when 5 or 6 basis functions are used





Introducing the arbitrary weighted function
Weak form can be found as follows, with weighted function

Now, multiply eq. 5.2.1 by $$ \displaystyle w $$ and integrate between 0 to 1.

Use integration by parts for the first term in eq. 5.2.14

Using eq. (5.2.2) and eq. (5.2.13), we get

Put eq. (5.2.16) in eq. (5.2.15), we get the weak form as follows

Using the polynomial basis function

For getting appropriate constraint breaking solution, we choose Also, Take $$ \displaystyle \beta=1$$ for the constraint breaking solution. Therefore, From eq. 5.2.2 and eq. 5.2.20, we can say that Therefore, eq. 5.2.20 becomes Similarly, the weighted function can be written as Eq. 5.2.13 gives, Therefore eq. 5.2.23 becomes Hence, from eq. 5.2.22 and 5.2.25 Put this in eq. 5.2.17, we get Using eq. 5.2.25 we can find $$ \displaystyle w(0) $$ and put in eq. 5.2.28, we get By integrating and verifying with| Wolfram Alpha, we get Following Matlab code is used to find coefficients a_j by changing values of i and j in the above equation The figure shows comparison of analytical and numerical solution for polynomial function. At n=5,numerical solution approaches analytical solution.Error decreases below $$ 10^{-6}$$ for n=10

Exponential basis function
Using the exponential basis function

For getting appropriate constraint breaking solution, we choose Also, Take $$ \displaystyle \beta=1$$ for the constraint breaking solution. Therefore, From eq. 5.2.13 and eq. 5.2.35, we can say that Therefore, eq. 5.2.35 becomes Similarly, the weighted function can be written as Eq. 5.1.3 gives, Therefore eq. 5.2.38 becomes Hence, from eq. 5.2.37 and 5.2.40 Put this in eq. 5.2.17, we get Using eq. 5.2.40 we can find $$ \displaystyle w(1) $$ and put in eq. 5.1.18, we get By integrating and verifying with| Wolfram Alpha1 ,| Wolfram Alpha2,| Wolfram Alpha3 we get

Following Matlab code is used to find coefficients a_j by changing values of i and j in the above equation After getting the coefficients of numerical solution, we compare numerical solution with analytical solution for different values of n.

Then graph of error that is difference between analytical solution and numerical solution is plotted at x=0.5



As expected, with increase in number of n,numerical solution approaches analytical solution.

From the above graph, error is almost zero for n=6 and above. But to get the error in the order of less than $$10^{-6}$$,we solve it further for higher values of n.

=Problem 5.3 Using Lagrange Interpolation Basis Functions to Approximate U(x)=

Given
Refer to lecture slide [[media:Fe1.s11.mtg29.djvu|29-5]] for problem assignment. Given the following Strong Form of a PDE:

Find
1) Solve (Eq. 5.3.1) using the weak form with Lagrange Interpolation Basis Functions (LIBF), equidistant nodes, and m = 4, 6, 8… until convergence of $$ u^{h} \left (0.5 \right ) $$ to $$ O \left ( 10^{-6} \right ) $$.

2) Explain how LIBF are used as Constraint Breaking Solutions.

3) Plot all LIBF used.

4) Plot $$ u \left ( x \right ) $$, $$ u(x)^{h} $$, and convergence (error vs. m)

Solving the PDE using the Exact Solution
First we must obtain the exact solution of (Eq. 5.3.1) in order to determine the error of our approximated function. To determine the exact solution, we take the integral of (Eq. 5.3.1.a) to obtain the following equation

Next, we rearrange (Eq. 5.3.2) and apply the natural boundary condition (Eq. 5.3.1.b) in order to solve for the constant A_1.

Therefore,

Substitute $$ A_1 $$ back into (Eq. 5.3.2), rearrange to isolate the term $$ \frac{du}{dx} $$ and then integrate in order to determine $$ u_{} $$.

Now we substitute the essential boundary condition (Eq. 5.3.1.c) into (Eq. 5.3.5) to solve for the constant $$ A_2 $$.

Solving for $$ A_2 $$ yields

Finally, we substitute (Eq. 5.3.7) into (Eq. 5.3.5) to find the exact solution for $$ \displaystyle u $$.

Determining the Weak Form
Now that we have our exact solution, we can begin to derive our approximated function. The approximated function is found by first obtaining the weak form of (Eq. 5.3.1).

First we multiply the governing partial differential equation (Eq. 5.3.1.a) and the natural boundary condition (Eq. 5.3.1.b) by an arbitrary weight function and integrate over the x domain.

Additionally, the weight function must disappear at the point where the essential boundary condition is applied, therefore,

Using integration by parts on the first term in (Eq. 5.3.9) yields

Substitute (Eq. 5.3.11) into (Eq. 5.3.9) and apply (Eq. 5.3.10) to obtain

Finally, we substitute in the natural boundary condition (Eq. 5.3.1.b) into (Eq. 5.3.12) and rearrange the terms to obtain the weak form of the PDE given in (Eq. 5.3.1)

Approximating U and W using LIBF
Our solution $$ u \left ( x \right ) $$ and our weight functions can be approximated by the following equations:

where $$ c_{i} $$ and $$ d_{j} $$ are constants with respect to x and $$ b_{i} $$ and $$ b_{j} $$ are the LIBF. To approximate our solution u, we divide our x domain into m subdivisions. Additionally, we will use m LIBF to approximate our solution. These LIBF will act as a shape function so that each LIBF will equal 1 at one specific point x, and will equal 0 at all other points, i.e.,

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

Additionally, the weight function must disappear at the essential boundary condition, therefore,

Also, the essential boundary condition (Eq. 5.3.1.c) can be substituted into (Eq. 5.3.14.a) to obtain

Next, we expand (Eq. 5.3.14) into 4 separate terms, i.e.,

Substitute the LIBF into (Eq. 5.3.22) and then take the derivative with respect to x to obtain the following:

Finally, we substitute (Eq. 5.3.20) through (Eq. 5.3.23) into (Eq. 5.3.13). The result can be separated into separate terms multiplied by $$ c_{i} $$.

Equation $$\displaystyle (Eq. 5.3.24)$$ must hold for any arbitrary weight function $$w(x)$$, so long as it satisfies the homogeneous essential boundary condition. It then also follows that it must also hold for any arbitrary $$c_2$$, $$c_3$$, and $$c_4$$, leading to the following algebraic equation:

Pre-multiplying both sides of $$\displaystyle (Eq. 5.3.25)$$ by the inverse of the coefficient matrix on the left side yields:

(Eq. 5.3.28) is combined with (Eq. 5.3.21) and the LIBF to create our approximation for the solution to the PDE, $$ u \left ( x \right ) $$, i.e.,

The same process can be repeated for m=6 and m=8. Results of these further iterations are shown in the sections below.

Explain how LIBF are used as Constraint Breaking Solutions
A Constraint Breaking Solution uses basis functions that abide by the following property,

$$ b_{1}(\beta)\neq 0 $$ and $$ b_{i}(\beta)=0 \quad\text{for}\quad \; i=2,3,...,n $$

where $$ \beta $$ is the location where the essential boundary condition is dictated ( $$ \beta = 0 $$ in our instance).

Upon closer analysis of LIBF, refer to (Eq. 5.3.16-19) for example, it becomes clear that these basis functions abide by this rule. For instance, for basis function $$ b_{1}$$ taken at the point where the essential boundary condition is defined, $$ b_{1} = L_{1,4} \left ( x=0 \right ) = 1 $$. Furthermore, all other LIBF $$ \left ( L_{2,4}, L_{3,4}, \quad\text{and}\quad L_{4,4} \right ) $$ are equal to zero at $$ x=0 $$ because $$ x_{1}=0 $$.

Therefore, substituting physical values of x into our LIBF illustrates that LIBF can be used as Constraint Breaking Solutions.

Plot all LIBF used






Plot u(x), u(x)^{h}, and convergence (error vs. m)


Notice that the approximations and the exact solution are very similar. The differences in the equations are too small for this figure to display.



The graph of the error vs m helps us understand the convergence of our approximated functions as m increases in value. The more exact values of error are:

m=4: 0.0158

m=6: -5.7975e-004

m=8: 2.1571e-005

Thus, we have convergence of our solution on the order of 10^-5 using m=8.

Matlab Code
=Problem 5.4 Using LIBF with equidistant nodes to solve G1DM1.0/D1=

GIVEN
One dimensional Model 1.0/Data set 1.0 on page 2 of Mtg 9:

where $${a_2}\left( x \right){\text{ =  2 + 3x}}$$,$$f\left( {x,t} \right) = 5x $$

1-D LIBF on page 3 of Mtg 26:

FIND
A) Explain how LIBF are used as CBS

B) Plot all LIBF used

C) Use matlab to integrate and solve the weak form

D) Plot $$ u_m^h\ vs\ u$$ and $$ u_m^h\left(0.5 \right)-u\left(0.5 \right)\ vs\ m$$

SOLUTION
A) Take m=4 for example to explain how LIBF works as CBS.

When m=4, we can write the basis function according to Eq. 5.4.4 (N was denoted as basis function from now on):

According to Eq.5.4.2, we should have CBS when x=1, which means that the last basis function should not equal zero when x=1.

At the same time, when x=1, we have the following values for LIBF based on its own properties:

Since we have the weight function as $$ w\left( x \right) = {c_i}N\left( x \right) $$ and from Eq. 5.4.9 we can choose arbitrary value for c from i=1~3 and make c equal zero for the last term, which absolutely satisfy CBS for which have non-zero value at the essential boundary for the last basis function and equal zero for all the others.

B)-D)Integrate Eq. 5.4.1 by parts over the whole domain, we can have

$$ \left[ {w(x)\left( {(2+3x)\frac} \right)} \right]_{}^{x = 1} - \left[ {w(x)\left( {(2+3x)\frac} \right)} \right]_{x = 0}^{} - \int\limits_0^1 {\frac\left( {(2+3x)\frac} \right)} dx + \int\limits_0^1 {w(x).(2+3x)} dx = 0 $$

Since LIBF satisfies CBS, we have w=0 when x=1, and then the above equation can give us the weak form:

If we use LIBF both as trial solution and weight function, the following matlab codes can solve the weak form for m nodes







The above three images show that whatever the value m equals, LIBF show its property which can be used as CBS, since its values vanish all other nodes except just one node where the natural boundary condition happens.



It is not easy to find out the difference from the approximation solution and exact solution, but the gaps do exist since the error is non-zero.



We can find out when the solution convergence and what is the error at x=0.5. The error is so small which indicates that LIBF, as a global basis function, has advantages solving weak form and obviously can satisfy CBS wherever essential boundary happens.

=Problem 5.5 Calculix Guide (Continued)=

Given
Refer to lecture slide [[media:Fe1.s11.mtg29.djvu|29-7]] for problem assignment. The Calculix Guide was originally started in Homework Problem 4.7.

Find
1) For the disk problem in Homework Problem 4.7, extract the node information (node numbers and coordinates) and the element information (element numbers and element nodes).

2) Generate 3 meshes of the same disk with triangular elements (increase the number of elements).

3) Install ccx, run examples, and write a report for "dummies" (explain commands, screenshots, etc...).

1.Extracting Node and Element Information from Calculix
For this tutorial, we will continue with the source code for a disk given in Homework 4.7.

First, we pre-process this code using the keystroke command F10. The disk is then graphically generated.

Now, we can obtain node information by performing the following tasks:

1. Plot the nodes and their global node number labels by entering the command "plot na all" and pressing enter. This is shown in Figure 5.5.1.



2. Next, we change our cursor into a box that can capture data within it. This is done by defining a set using the command "qadd". Type "qadd" and then press enter. Notice that your cursor changes to a miniature box. This is shown in Figure 5.5.2.



3. Now, we redefine the shape of our data capturing box to a more convenient size by typing "r". This will define the upper left corner of our new data capturing box. Move your cursor to where you would like the lower right corner of your data capturing box and type "r" again. Notice that your data capturing box has changed its size to fit the corners defined by your "r" commands. This is shown in Figure 5.5.3.



4. Finally, to capture the requested node information (the node number and its coordinates) simply mouse over a node of interest. Once the node is within your data capturing box, simply press "n". This command will output the node information to your output screen. This is shown in Figure 5.5.4. Once you are done capturing the node information, press "q" to quit using the data capturing box.



We can obtain element information (element number and element nodes) by performing the following tasks:

1. Plot the elements of our disk by typing in the command "plot ea all" and pressing the enter key. Notice that the disk looks significantly different than when the nodes were plotted. This view highlights the elements used in our model. See Figure 5.5.5.



2. Next you set up your data capturing box. This method is the same method that was utilized when finding the node information (outlined in steps 2 and 3 of the node information section).

3. Finally, to capture the requested element information simply mouse over an element of interest. Once the element is within your data capturing box, simply press "e". This command will output the element information to your output screen. This is shown in Figure 5.5.6. Once you are done capturing the element information, press "q" to quit using the data capturing box.



2.Generate three meshes of same disc using triangular elements
Meshes are redefined using commands:

'elty3'

'elty6'

and also by redefining the points







3.Installation procedure and report for dummies
To install CalculiX for Windows, start by downloading the executable file available at http://www.bconverged.com/download.php.

Before installing CalculiX, we recommend that you create a folder that you will want to use as your workspace for CalculiX files.

This executable file installs both-CGX and CCX on your computer.

The procedure to install can be found at Team 6-Homework 4.7 []

Introduction

The basic procedures and commands can be found at CCX-help menu as shown below:



You can also refer to ccx manual below for further details:


 * 1) CCX Manual []

=Problem 5.6 Quadratic Lagrangian Element Basis Functions=

Find
Plot three figures similar to those of the linear Lagrangian element basis functions..

Solution
The generic equation for the quadratic Lagrangian element basis functions (QLEBF) is presented below in Equation 5.6.1;


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

$$ L_{j,n}=\prod_{k=1}^{n}\frac{(x-x_{k})}{(x_{i}-x_{k})} $$
 * $$\displaystyle (Eq.5.6.1) $$
 * }
 * }

Equation 5.6.1 can be expanded to the following equation to fit each of our cases as shown in Equations 5.6.2-4 below with x1=1, x2=2, and x3=3.


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

$$ L_{1,3}=\frac{(x-x_{2})(x-x_{3})}{(x_{1}-x_{2})(x_{1}-x_{3})} $$
 * $$\displaystyle (Eq.5.6.2) $$
 * }
 * }


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

$$ L_{2,3}=\frac{(x-x_{3})(x-x_{1})}{(x_{2}-x_{3})(x_{2}-x_{1})} $$
 * $$\displaystyle (Eq.5.6.3) $$
 * }
 * }


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

$$ L_{3,3}=\frac{(x-x_{1})(x-x_{2})}{(x_{3}-x_{1})(x_{3}-x_{2})} $$
 * $$\displaystyle (Eq.5.6.4) $$
 * }
 * }

Since the quadratic Lagrangian element basis functions have three nodes, three QLEBFs were used to generate the following plots 5.6.1-4 in MATLAB. The code used to generate these plots is presented below.









Comments
Each of the plotted quadratic Lagrangian element basis functions has three nodes as opposed to the two nodes of the linear basis functions due to the n=3 condition. It is observed from the plots that the sum of the basis functions at each node must be equal to 1 regardless of the number of nodes. This also holds true for the linear Lagrangian element basis functions as well.

=Problem 5.7:Solving a strong form with Linear Lagrangian Element Basis Function (LLEBF)=

Given
The given ordinary differential equation is:

Boundary conditions
The Essential Boundary condition is -- $$ u\left(0\right) = 4$$ The Natural Boundary condition is -- $$\frac{du}{dx}\left(1\right) = -6 $$

Descritization
Here we want to use Uniform Descritization such that Element nodes are equidistant with number of elements as $$ nel = 4, 6, 8, ..$$

LLEBF as Constraint Breaking Solution (CBS)
Consider a basis function $$b_{i}$$ such that we have, $$b_{0}\left(\beta\right) = 1$$, and $$b_{i}\left(\beta\right) = 0; i = 1, 2, ...$$. In such a case as per the discussion in Mtg 20, $$b_{i}$$ becomes a Constraint Breaking Solution (CBS)

LLEBFs (N1 and N2 below) for a 2 node linear element can be used to represent the unknown variable u(x)in a given domain which is written as -


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

$$u\left ( x \right )=u_{1}\frac{x-x_{2}}{x_{1} - x_{2}}+u_{2}\frac{x-x_{1}}{x_{2} - x_{1}} =u_{1}N_{1}+u_{2}N_{2} $$
 * $$\displaystyle (5.7.2) $$
 * }
 * }

Here $$u_{1}, u_{2}$$ are the nodal values of the unknown u(x) at the nodes 1 and 2 respectively.

The Lagrange Polynomials follow the property of Partition of Unity i.e. they have the value of 1 at the respective node and 0 at the others nodes. It means the basis function $$ b_{i} $$ carries value as 1 at the i th node and has a zero value at any other nodes in the domain.

Plot all LLEBF used for $$ nel =4 $$
The straight line plots (common for all elements used with $$ nel =4 $$) irrespective of the elemental sizes for the 2 LLEBFs given in Eq. 5.7.2, following the Partition of Unity model are shown in the Figure 5.7.1. Note- In the following Figure 5.7.1., N[i] and N[i+1] are the nodal basis functions and i's denote the nodes



Solve the system using (nel = 4) and increase the number of elements or $$ nel =4,6,8..$$ until the error at x = 0.5 converges to the order of 10^-6
The weighted residual form with weighting function $$ w (x) $$ is given as:

We develop the corresponding weak form using the technique of integration by parts, with $$ w(0) = 0$$ owing to the fact that the flux is unknown at $$ x = 0 $$. The expression for the weak form is as below -

To use the element stiffness matrices to the global stiffness matrix $$ K $$,

The stiffness matrix is obtained as below -

The force matrix is obtained as below -

We assemble the system using a Matlab code given below-

Plot for Exact solution $$ u (x)$$ and Approximate solution $$ u^h (x) $$ along x
It can be observed from the Figure 5.7.2 (Using 4 nodes) and Figure 5.7.3 (Using 6 nodes) that the exact solution is nearly linear and can be very well approximated by the LLEBFs considered with least 4 number of nodes and obviously also when the number of nodes is increased to 6. This is as the approximated solution $$ u^h (x) $$ very closely imitates the curve for the exact solution $$ u (x)$$.



Convergence - Plot of error vs. number of nodes (n)
At x = 0.5, to plot the error ( $$ u (x)$$ - $$ u^h (x) $$ ), we use Eq. 5.7.2 with the interpolation functions.

Even with (n = 4), sufficiently low number of elements, the system gives considerable convergence to an order of to 10^-3. We need more elements for it to further converge to the order of of 10^-6. Figure 5.7.4 depicts the error as a function of the number of elements.



=Problem 5.8: Solving Data Set G1DM1.0/D1 Using Linear Lagrangian Element Basis Functions and Uniform Discretization=

Boundary Conditions
The Essential Boundary condition as applied on the boundary $$\Gamma_g = [1]$$, is

The Natural Boundary condition as applied on the boundary $$\Gamma_h = [0]$$, is

Descritization
Here we want to use Uniform Descritization such that Element nodes are equidistant with number of elements as $$ nel = 4, 6, 8, ..$$

Find
1. Explain how Linear Lagrangian Element Basis Function (LLEBF) are used as Costraint Breaking Solution. 2. Plot all LLEBF used for $$ nel =4 $$. 3. Use matlab to numerically integrate and find the solution for $$n_{el}$$ = 2,4,6,8... until the error evaluated at x=0.5 converges to the order of $$10^{-6}$$ 4. Plot for Exact solution $$ u (x)$$ and Approximate solution $$ u^h (x) $$. 5. Convergence - Plot of error vs. number of nodes (n).

Linear Lagrangian Element Basis Functions as a Constraint Breaking Solution
Linear Lagrangian element basis functions assume the value of the Kronecker Delta at nodes, i.e,

For this case, where $$\displaystyle \Gamma_g = 1$$, define the terminating node, $$\displaystyle x_n = 1$$ so that all except the last of the LLEBF will be equal to zero when evaluated at this point while simultaneously the final basis function, $$ b_n$$ will be equal to unity, i.e.

This breaks the constraints of the weighting coefficients $$\displaystyle \{c_i,\;\;i = 0,2,...,n-1\}$$ and as such any arbitrary set may be chosen, as simply selecting $$\displaystyle c_n = 0$$ will satisfy the homogeneous essential boundary condition.

Solution of the System
In the element sense the basis functions are composed of linear Lagrange equations of the form

Now creating the element shape function matrix for this linear 2-node element:

Now, substituting $$\displaystyle l^e = x_2^e-x_1^e\;$$ into $$ \displaystyle (Eq.5.8.8)$$

The element shape function derivative matrix for these elements is

Now constructing the element stiffness matrix

The element force matrix is

Finally the global stiffness and force matrices are assembled via equations (5.13) and (5.14) (F&B p. 96)

where $$\displaystyle L^e$$ is the element dependent scatter operator.

Plot for Exact solution $$ u (x)$$ and Approximate solution $$ u^h (x) $$ along x
It can be observed from the Figure 5.8.2 (Using 4 nodes) and Figure 5.8.3 (Using 6 nodes) that the exact solution is nearly linear and can be very well approximated by the LLEBFs considered with least 4 number of nodes and obviously also when the number of nodes is increased to 6. This is as the approximated solution $$ u^h (x) $$ very closely imitates the curve for the exact solution $$ u (x)$$.



Convergence - Plot of error vs. number of nodes (n)
Even with (n = 4), sufficiently low number of elements, the system gives considerable convergence to an order of to 10^-3. We need more elements for it to further converge to the order of of 10^-6. Figure 5.7.4 depicts the error as a function of the number of elements.



=Contributing Members=

Eml5526.s11.team6.joglekar 16:38, 23 March 2011 (UTC)

eml5526.s11.team6.lee 23:45, 23 March 2011 (UTC)

Eml5526.s11.team6.tupsakhare 01:30, 24 March 2011 (UTC)

Eml5526.s11.team6.kurth 02:02, 24 March 2011 (UTC)

Eml5526.s11.team6.deshpande 02:12, 24 March 2011 (UTC)

Eml5526.s11.team6.vork 02:37, 24 March 2011 (UTC)

Eml5526.s11.team6.gravois 03:14, 24 March 2011 (UTC)