User:Eml5526.s11.team01/Hwk5

 HOMEWORK 5

=Problem 5.1: Weak form with appropriate functions such as Polynomial, Fourier & Exponential until convergence.=

Given
General 1-D Model 1 Data set 1. The Partial Differential equation for this data set is,

Find
Solve the Weak Form of the above general model 1 data set 1 for the Functions such as

1)Polynomial Basis Function,

2)Exponential Basis Function &

3)Fourier Basis Function

to find out the convergence point.

1) Polynomial Basis Function
Solve the Differential Equation

in the interval (0,1) with the boundary conditions

The Weak Form of the given Differential Equation will be obtained as follows

So we need to integrate the above equation using integration by parts, after which we will get the following equation,

We now have to select the appropriate Polynomial Basis Function which will follow the constraint breaking solution. so selecting

Clearly,

We need, So Choosing,

We now have, Using the essential Boundary condition we have, Therefore, In the identical way we can also write,

With the condition,

Therefore,

Now we need to calculate the derivative of $$U$$ & $$w$$, Therefore,

Substituting this into the weak form of $$\left(5.1.9\right)$$, we get,

Cancelling $$c_{i}$$ and Organizing, we get, Simplifying the last equation we get,

We can obtain the value of $$\displaystyle a_{j}$$ from the above equation.

The analytical exact solution on the other hand is obtained as,

therefore,

The Values of p & r are obtained by using the following boundary conditions,

To solve this system of equations we developed a MATLAB code. It was found out that with $$n=2$$ a good approximation is obtained. However with$$n=8$$ the error was smaller than $$10^{-6}$$. The MATLAB code is as follows.

Appendix
 Matlab Code: 

The analytical and the numerical solution along with the Error between the two are plotted as below.





2) Fourier Basis Function
To follow the constraint Breaking Solution we choose,

We need to satisfy the essential boundary condition So we choose the basis Function such that So we can write, Now the trial solution will be as follows,

Now we need to solve for the remaining coefficients By constructing a Stiffness & Force Matrices as follows.

Here the first equation gets satisfied for even values of i & j,

Whereas 2nd equation is satisfied for odd values of i & j,

further 3rd equation is satisfied for remaining combinations of i & j.

Here the first equation gets satisfied for even values of i & j,

Whereas 2nd equation is satisfied for odd values of i & j,

We developed a MATLAB code to generate the Force and the stiffness Matrix and We get the following values.

We can now use the relationship which is, $$\displaystyle Kd=F$$ to solve for $$\displaystyle d$$.

Ultimately the trial solution of the equation is,

We have developed at MATLAB code to solve these system of equations & to deliver the Error between analytical and the numerical solution.

Appendix
 Matlab Code: 

The analytical and the numerical solution along with the Error between the two are plotted as below. From the Plots we can explicitly say that the solution converges at n=9.





3) Exponential Basis Function
To follow the constraint Breaking solution we choose: Clearly, And as we need We choose, Eventually we have, From the essential boundary condition, we have, Therefore, In the identical way we have, Using the condition, Therefore Derivative of U & w will be,

Replacing in the weak form of $$\left(5.1.5\right)$$ we get,

Cancelling $$c_{i}$$ and then organizing we obtain, Performing the integration we obtain following equation,

The value of the $$a_{j}$$can be calculated from the above equation. we developed a MATLAB code to find the error between the Exact and the analytical solution. The MATLAB code is as follows.

Appendix
 Matlab Code: 

The following two figures shows the basis Function and the Error between the Analytical and the numerical solution. It can be observed that the Exponential basis Function converges somewhat Slowly than the Polynomial basis Function. For exponential the convergence was achieved with 10 terms whereas with polynomial, it was achieved with only 8 functions. Moreover the stiffness matrix obtained with exponential functions was bad-conditioned for a number of functions greater than 8.





=Problem 5.2: Weak form with appropriate functions such as Polynomial, Fourier & Exponential until convergence.=

Given
General 1-D Model 1 Data set 1. The Partial Differential equation for this data set is,

Find
Solve the Weak Form of the above general model 1 data set 1 for the Functions such as

1)Polynomial Basis Function,

2)Fourier Basis Function &

3)Exponential Basis Function

to find out the convergence point.

Solution
Solve the differential equation $$ \frac{d}{dx} \left[ \left(2+3x \right)  \frac{du}{dx} \right]+5x = 0 $$ in the interval (0,1) with the boundary conditions

$$ u\left(1 \right)= 4 ; \frac{du}{dx} \left(0 \right)= -6 $$

The weak form was obtained in the latter section:

As we do not know $$ \frac{du}{dx} \left(1\right) $$ then we chose $$ w \left(1 \right) = 0 $$ ; moreover $$ u^h \left(1 \right)= 4 $$

1) Polynomial Basis Function
Polynomial basis functions To follows the constrain braking solution, we chose:

Clearly $$ b_1=1 $$, and as we need that $$ b_i \left(1 \right) = 0, i = 1,2,... $$;we chose $$ \beta=1 $$. Then we have $$ U^h \left(x\right)= a_0+a_1 \left(x-1 \right)+a_2 \left(x-1\right)^2+... $$

From the essential boundary condition we have: $$ U^h \left(1 \right)= 4= a_0 $$

Therefore:

In a similar way we have

With the condition: $$ w \left(1 \right)= 0= c_0 $$

Therefore:

The derivative of U and w are:

Replacing in the weak form we have:

Canceling $$ c_i $$ and Organizing we obtain:

Performing the integration we obtain:

From the latter equation the coefficients $$ a_j $$ can be obtained.

The analytical general solution is:

With the boundary conditions we obtain the constants:





 Matlab Code: 

2) Fourier Basis Function
Now the differential equation is solved using Fourier basis function and the boundary conditions:

To follow the constrain braking solution we chose:

For Cosine Clearly $$ b_1=1 $$, and as we need that $$ b_j \left(1 \right) = 0,j=0,1,2,... ~$$ ; we chose $$ \beta = 1  $$.

For Sine Clearly $$ b_1=1 $$, and as we need that $$ b_k \left(1 \right) = 0,k=1,2,... ~$$ ; we chose $$ \beta = 1  $$.

Then we have $$ U^h \left(x \right)= a_0 cos \left(x-1 \right)-1 + a_1 cos2 \left( x -1 \right)-1 ... $$

$$ U^h \left(x \right)= a_1 sin \left(x-1 \right) + a_2 sin2 \left( x -1 \right) ... $$

From the essential boundary condition we have: $$ U^h \left(1 \right)= 4=a_0 $$

Therefore:

By constructing a Stiffness & Force Matrices as follows.

Here the first equation gets satisfied for even values of i & j,

Whereas 2nd equation is satisfied for odd values of i & j,

further 3rd equation is satisfied for remaining combinations of i & j.

Here the first equation gets satisfied for even values of i & j,

Whereas 2nd equation is satisfied for odd values of i & j,

We developed a MATLAB code to generate the Force and the stiffness Matrix and We get the following values.

We can now use the relationship which is, $$\displaystyle Kd=F$$ to solve for $$\displaystyle d$$.

Thus the trial solution is:



 Matlab Code: 

3) Exponential Basis Function
Now the differential equation is solved using exponential basis function and the boundary conditions:

To follow the constrain braking solution we chose:

Clearly $$ b_1=1 $$, and as we need that $$ b_i \left(1 \right) = 0,i=1,2,... ~$$ ; we chose $$ \beta =1 $$.

Then we have $$ U^h \left(x \right)= a_0+a_1 \left(e^{\left(x-1 \right)}-1\right)+a_2 \left(e^{2\left(x-1\right)}-1\right)+... $$

From the essential boundary condition we have: $$ U^h \left(1 \right)= 4=a_0 $$

Therefore:

In similar way we have:

With the condition: $$ w(1) = 0 = c_1 $$

Therefore:

The derivative of U and w are:

Replacing in the weak form we have:

Canceling $$ c_i $$ and organizing we obtain:

Performing the integration we obtain:

The following figure shows the exact and numerical solutions. Although good approximation was obtained at $$ n = 3 $$, the convergence for an error of order $$10^{-6} $$ was only achieved with $$ n=29 $$ and to $$n $$ greater than $$ 8 $$, a bad-conditioned (close to singular) stiffness matrix was obtained. Figure 2.4 shows the error.

This results shows that polynomial is better basis function because the convergence is obtained faster; moreover for other problems, the integration may not be possible to perform analytically as in this case, making necessary that the program perform the integration numerically which will increase the time of numerical calculations and decrease the accuracy.



 Matlab Code: 

=Problem 5.3: Solve the General 1-D Model 1.0/Data 1b using Weak Form with the Lagrange Interpolation Basis Functions (LIBF)=

Explain how the LIBF basis functions are used as the constraint breaking solution (CBS)
The basis functions are the following:

Solution
1.

So for the LIBF, at the essential boundary condition, there is only one term which does not equal to zero, the other terms are equal to zero. So the CBS is satisfied.

2.

Plot all the LIBF used:

For m=4:



m=6:



m=8：



3.

Solve for d, we have:

=Problem 5.4: Solve the General 1-D Model 1.0/Data 1 using Weak Form with the Lagrange Interpolation Basis Functions (LIBF)=

Explain how the LIBF basis functions are used as the constraint breaking solution (CBS)
The basis functions are the following:

Solution
1.

So for the LIBF, at the essential boundary condition, there is only one term which does not equal to zero, the other terms are equal to zero. So the CBS is satisfied.

2.

Plot all the LIBF used:

For m=4:



m=6:



m=8：



3.

Solve for d, we have:

So, we have:

4.

m=4:

m=6:

m=8：

=Problem 5.5: Calculix=

Given
Given the disk problem in homework four [| Working with Calculix]

Find

 * 5.1 Extract node information to include node numbers and coordinates, and element information to include element numbers and element nodes.
 * 5.2 Generate 3 meshes of same disk with triangular elements.
 * 5.3 Install ccx, run examples, write report for “dummies”

5.1 Node and Element information
A disk was created in Problem 4.7 in the previous homework. CGX was opened with this file, and the command “prnt se 1” was used to extract the information. This was copied to a file and is available at the following link http://en.wikiversity.org/wiki/File:Fe1.s11.team1.hw5.node.element.info.pdf

5.2 Three Meshes
Element type was changed using the command “elty all tr3”, “etly all tr3u” and “elty all tr6” and then using the “mesh all” command to mesh the disk. Refer to figures 5.1, 5.2 and 5.3

5.3 CCX
CCX was installed and opened. The manual spells out the program well; this was downloaded, but can be found at http://web.mit.edu/calculix_v2.0/CalculiX/ccx_2.0/doc/ccx/index.html. Examples were run, please see the following report for findings. http://en.wikiversity.org/wiki/File:CCX.pdf

Contributing Members
Posted by James Roark 01:59, 23 March 2011 (UTC)

=Problem 5.6: Quadratic Lagrangian Elemental Basis Function.=

Given
Consider any Quadratic Function.

Find
Plot of the Individual Lagrangian Quadratic Function and the combined plot of all the three Lagrangian Quadratic Functions.

Solution
For any Quadratic Function the Number of nodes are always 3. Therefore we can write, We will select the global node numbers as follows, We have according to the Lagrangian Interpolation Basis Function, Which is the Polynomial of degree (n-1)

So, Substituting the values from Equation $$\displaystyle (5.6.2)$$, we get,

In the Similar way,

&

We Developed a MATLAB code to plot all these Quadratic Lagrangian Basis Function. This very MATLAB code is as follows.

Appendix
 Matlab Code: 

The plots of all the Individual Lagrangian Quadratic Function and the Combined plot of all the Three Functions are as follows.









=Problem 5.7: Linear Lagrangian Element Basis Function= Solving a ordinary differential equation using Linear Lagrangian Element Basis Function. From Mtg.30

Given
The partial differential equation is:

The boundary conditions are:

$$ u\left(0\right) = 4$$

$$\frac{du}{dx}\left(1\right) = -6 $$

Find

 * 1. Explain how Linear Lagrangian Element Basis Function (LLEBF) are used as Costrain Breaking Solution
 * 2. Plot all LLEBF
 * 3. Solve the system dividing the interval [0,1] in 3 divisions (m = 4) and continue dividing until convergence of 10^-6 of the error at x = 0.5
 * 4. Plot the exact and the aproximated solution
 * 5. Plot the error as a function of the number of points

Explain how LLEBF are used as CBS
A function $$b_{i}$$ is used as CBS if: $$b_{0}\left(\beta\right) = 1$$, and $$b_{i}\left(\beta\right) = 0; i = 1, 2, ...$$.

A funtion u(x) can be interpolated by using linear lagrange polynomial considering an element of two points (nodes): 1 and 2; the function will be given by:


 * {| 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) $$
 * }
 * }

Where $$u_{1}, u_{2}$$ are the nodal values of the funtion at the nodes 1 and 2 respectively.

The Lagrange Polynomial has the property of have the value of 1 at the respective node an 0 at the others; for example in equation (5.7.2), if an esential boundary condition is applied on node 1, whe have that $$N_{1} = b_{1} = 1$$ and $$N_{2} = b_{2} = 0$$ at the node 1 where the essential boundary condition is applied. In this case we only have two fuctions (b1 and b2); therefore this polynomial follows CBS. Symilar aproach can be used if the essential boundary condition is at node 2.

Plot all LLEBF used
Fig. 5.7.1 shows the two basis functions described in Eq. (5.7.2); in the figure, clearly can be seen that N1 is equall to 1 at node 1 and 0 at node 2, and N2 is equal to 0 at node 1 and equal to 1 at node 2. The plot are obviously straight lines because the functions are linear.



In the present problem we have to start working with 4 poins (3 elements) and increase the number of points until that the error at x = 0.5 is smaller than 10^-6. We are asked to plot all the basis functions used, but the basis funtion used are similar for all elements does not matter the size of the element. When the number of elements is increased, the lenght of the elements is decreased, but the basis fuction are the same as that shown in Fig. 5.7.1 for all the elements. Obviously a function is valid only inside its respective element

Solve the system dividing the interval [0,1] in 3 divisions (m = 4) and continue dividing until convergence of 10^-6 of the error at x = 0.5
Firs we obtain the weak form of the differential equation; to do that we multiply it for a weight funtion and integrate from 0 to 1:


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

$$\int_{0}^{1}w\left ( \frac{d}{dx} \left [ \left ( 2+3x \right ) \frac{du}{dx}\right ]+5x\right )dx=0 $$
 * $$\displaystyle (5.7.3) $$
 * }
 * }

Performing integration by parts in the first term we obtain:

As we do not know the derivative of u at x = 0, we chose w(0) = 0. The global stiffness matrix will be assembled from the stiffness matrix of each element. To get the element stiffness matrix we change the cordinate x to the cordinate $$\epsilon$$ which is equal to -1 at the node 1 and is equal to 1 at the node 2, as indicated in the Fig. 5.7.2

The new cordinate $$\epsilon$$ is given by:


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

$$\varepsilon = \frac{2\left ( x-x_{1} \right )}{x_{2}-x_{1}}-1=\frac{2\left ( x-x_{1} \right )}{l}-1 $$ In the last step, it was used the fact that $$x_{2}-x_{1} = l$$; the lengt of the element. Taking the derivative on both sides of equation (5.7.5) we have:
 * $$\displaystyle (5.7.5) $$
 * }
 * }


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

$$ d\varepsilon =\frac{2dx}{l}$$
 * $$\displaystyle (5.7.6) $$
 * }
 * }

Replacing the cordinate x in the lagrange polynomial we obtain:


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

$$u\left ( \varepsilon \right )=\begin{bmatrix} \frac{1-\varepsilon }{2} & \frac{1+\varepsilon }{2} \end{bmatrix}\begin{bmatrix} u_{1}\\

u_{2}\end{bmatrix}=\left [ N \right ]\begin{bmatrix} u_{1}\\

u_{2}\end{bmatrix} $$
 * $$\displaystyle (5.7.7) $$
 * }
 * }

In symilar way w can be written as:


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

$$\left [w\left(x\right)= N \right ]\begin{bmatrix} c_{1}\\

c_{2}\end{bmatrix}=\begin{bmatrix} c_{1} & c_{2} \end{bmatrix}\left [N \right ]^{T} $$ As in the weak form, we need the derivatives of u and w with respect to x, we calculathe those derivatives by using the chain rule:
 * $$\displaystyle (5.7.8) $$
 * }
 * }


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

$$\frac{du}{dx}=\frac{du}{d\varepsilon }\frac{d\varepsilon }{dx}=\frac{du}{d\varepsilon }\frac{2}{l} $$
 * $$\displaystyle (5.7.9) $$
 * }
 * }


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

$$\frac{dw}{dx}=\frac{dw}{d\varepsilon }\frac{d\varepsilon }{dx}=\frac{dw}{d\varepsilon }\frac{2}{l} $$ From Eq. (5.7.7) we can take the derivative with respect to \epsilon to obtain:
 * $$\displaystyle (5.7.10) $$
 * }
 * }
 * {| style="width:100%" border="0"

$$\frac{du}{d\varepsilon }=\begin{bmatrix} -\frac{1}{2} & \frac{1}{2} \end{bmatrix}\begin{bmatrix} u_{1}\\

u_{2}\end{bmatrix}=\left [ B \right ]\begin{bmatrix} u_{1}\\

u_{2}\end{bmatrix} $$ In similar way the derivative of w is:
 * $$\displaystyle (5.7.11) $$
 * }
 * }


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

$$\frac{dw}{d\varepsilon }=\begin{bmatrix} c_{1} & c_{2} \end{bmatrix}\left [ B \right ]^{T} $$
 * $$\displaystyle (5.7.12) $$
 * }
 * }

Making the change of variables on the stiffness matrix and on the force vector we get:

To assemble the system we build a Matlab code (shown below) wich basically follow the algoritm described on Fig. 5.7.3. Inside the code, the integration was done by using the Gauss quadrature method. In that method, the integral of a funtion is given by:


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

$$ \int_{-1}^{1}f\left ( \varepsilon \right )d\varepsilon =\sum p_{i}f\left ( \varepsilon _{i} \right ) $$
 * $$\displaystyle (5.7.15) $$
 * }
 * }

Where $$p_{i}$$ are weighting coefficients. As in Equation (5.7.14) we have a polynomial of order 2, we need at least 2 coefficients to exactly calculate the integral. In this case the integral is given by:


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

$$ \int_{-1}^{1}f\left ( \varepsilon \right )d\varepsilon =f\left ( \frac{1}{\sqrt{3}} \right )+f\left ( \frac{-1}{\sqrt{3}} \right ) $$
 * $$\displaystyle (5.7.16) $$
 * }
 * }



Plot the exact and the aproximated solution
Fig. 5.7.4 and Fig. 5.7.5 shows the exact and aproximated solution obtained with 4 and 6 nodes respectively. It can be observed that even with 4 nodes the numerical solution aproximates well to the analytical.



To calculate the error, at x = 0.5 is the mid point of the domain, we obtain the value at x = 0.5 interpolating the two neighboring values. The interpolation was done by using the same Lagrange polynomial described before. We found that even though a good aproximation was obtained even with low number of elements (n = 4, Fig. 7.5.4) the convergence to a tolerance of 10^-6 was not obtained with a number of elements lower than 100. Therefore we decrease the tolerance to 10^-3. With this this tolerance the convergence still was achieve with a large number of elements (19) Fig.7.5.6 shows the error as a funtion of the number of elements.

Appendix
 Matlab Code: 

=Problem 5.8: Linear Lagrangian Element Basis Function= Solving a ordinary differential equation using Linear Lagrangian Element Basis Function. From Mtg.30

Given
The partial differential equation is:

The boundary conditions are:

$$ u\left(1\right) = 4$$

$$\frac{du}{dx}\left(0\right) = -6 $$

Find

 * 1. Explain how Linear Lagrangian Element Basis Function (LLEBF) are used as Costrain Breaking Solution
 * 2. Plot all LLEBF
 * 3. Solve the system dividing the interval [0,1] in 3 divisions (m = 4) and continue dividing until convergence of 10^-6 of the error at x = 0.5
 * 4. Plot the exact and the aproximated solution
 * 5. Plot the error as a function of the number of points

Solution
This problem is similar to 5.7 but only change the boundary conditions, therefore, the weak form is the same exept that now we chose w(1)=0 and impose the boundary condition u(1) = 4. As the problem is similar, we use the same Matlab routine developed for porblem 5.7. As all the develpment is the same that in problem 5.7 here we only present the results. Fig. 5.8.1 and Fig. 5.8.2 show the analytical and numerical solution obtained with 4 and 6 nodes respectively; again, although good aproximation was obtained even with 4 nodes, no convergence was achieved even with 100 nodes; therefore, the tolerance was decreased to 10^-3. With this new tolerance, convergence was achieved with 8 nodes. Fig. 5.8.3 shows the error as funtion of the number of nodes.





=Contributing Members=