User:Eml5526.s11.team01/Hwk3

Problem 1: Solving using WRF

Solving a differential equation using Weight Residual Form. From lecture 13

Given

The partial differential equation is:

$\frac{\partial}{\partial x}\left(2\frac{\partial u}{\partial x}\right) + 3 = 0$

$\displaystyle (Eq.1.1)$

The boundary conditions are:

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

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

A weighting funtion should be used:

$b_{i} = b_{i}=\left ( x+1 \right )^{i}$

$\displaystyle (Eq.1.2)$

Find

1.1 Find an approximate solution of the form $U^{h}=\sum_{i=0}^{n}d_{i}b_{i}$ for n=2

1.2 Find two equations that enforce the boundary conditions

1.3 Project the weight residues

1.4 Display the equations in matrix form

1.5 Solve for $\displaystyle d$

1.6 Construct $U^{h}$ and plot $U^{h}$ and $u$

1.7 Repeat 1.1 to 1.6 for n = 4 and n = 6

Solution

1.1 Find an approximate solution with n=2

When n = 2 we have:

$d_{i}=\begin{bmatrix} d0 & d1 & d2 \end{bmatrix}$

$\displaystyle (Eq.1.3)$

and

$b_{i}=\begin{bmatrix} 1 & \left ( x+ 1\right) & \left ( x+ 1\right )^{2} \end{bmatrix}$

$\displaystyle (Eq.1.4)$

Therefore:

$U^{h}=d0+d1\left ( x+1 \right )+d2\left ( x+1\right )^{2}$

$\displaystyle (Eq. 1.5)$

1.2 Find two equations that enforce the boundary conditions

The boundary conditions are:

$u\left(1\right)=0=d0+2d1+4d2$

$\displaystyle (Eq. 1.6)$

And

$\frac{du}{dx}(0)=-4 = d1+2d2$

$\displaystyle (Eq. 1.7)$

1.3 Project the weight residues

In this case changing the function u(x) by the approximated function $U^{h}$ in the PDE (Eq.1.1), we found that P(u) is:

$P\left ( U^{h} \right )=2\frac{\partial^2 U^{h}}{\partial x^2}+3=2d2+3$

$\displaystyle (Eq.1.8)$

The (Eq.1.8) can be written as:

$P\left ( U^{h} \right )= \begin{bmatrix}0 & 0 & 2\end{bmatrix}\begin{bmatrix} d0\\ d1\\ d2\end{bmatrix}+3$

$\displaystyle (Eq.1.9)$

Projecting the residue we have:

$\int_{a}^{b}b(x)P\left ( U^{h} \right )= 0$

$\displaystyle (Eq. 1.10)$

After substitution of (Eq.1.4) and (Eq.1.9) in (Eq.1.10) the following equation is obtained

$2\int_{0}^{1}\begin{bmatrix} 1\\ \left (x+ 1\right ) \\ \left ( x+1\right )^{2} \end{bmatrix}\begin{bmatrix} 0 &0 & 2\end{bmatrix}dx\begin{bmatrix} d0\\ d1\\ d2\end{bmatrix}=-3\int_{0}^{1} \begin{bmatrix} 1\\ \left ( x+1\right ) \\ \left ( x+1 \right )^{2}\end{bmatrix}dx$

$\displaystyle (Eq. 1.11)$

1.4 Display the equations in matrix form

Performing the product and then the integration indicated in (Eq.1.11) we have:

$\begin{bmatrix} 0 &0 & 4 \\ 0 & 0& 6 \\ 0 & 0& 3\frac{28}{3}\end{bmatrix}\begin{bmatrix} d0\\ d1\\ d2\end{bmatrix}=\begin{bmatrix} -3\\ -\frac{9}{2}\\ -21\end{bmatrix}$

$\displaystyle (Eq. 1.12)$

The latter equation is clearly of the form $\displaystyle Kd=F$. We can observe that K is not symmetric.

1.5 Solve for $d$

In (Eq.1.12) we have three equations; but, as we need to enforce the boundary conditions we take the only one of them and solve it together with (Eq.1.6) and (Eq.1.7). From these equations we now have the system:

$\begin{bmatrix} 1 &2 & 4 \\ 0 &1 &2 \\ 0 & 0& 4 \end{bmatrix}\begin{bmatrix} d0\\ d1\\ d2\end{bmatrix}=\begin{bmatrix} 0\\ -4\\ -3\end{bmatrix}$

$\displaystyle (Eq.1.13)$

After solving the system in (Eq.1.13) by using a calculator, we obtain d:

$\begin{bmatrix} d0\\ d1\\ d2\end{bmatrix}=\begin{bmatrix} 8\\ -2.5\\ -0.75\end{bmatrix}$

1.6 Construct $U^{h}$ and plot $U^{h}$ and $u$

Replacing d in (Eq.1.5) we obtain the approximated solution:

$U^{h}=8-2.5\left ( x+1\right )-0.75\left ( x+1\right )^2$

$\displaystyle (Eq.1.14)$

On the other hand, the exact solution of the PDE (Eq.1.1) with the given boundary condition is:

$u\left ( x \right )=-\frac{3}{4}x^{2}-4x+\frac{19}{4}$

$\displaystyle (Eq.1.15)$

We observe that the aproximated and the exact solution are the same; therefore, increasing the number of terms in b(x) will not improve the solution

In Fig.1.1 we show the exact and the approximated solution.

Fig.1.1 Analytical and numerical solution with i=2

1.7 Repeat 1.1 to 1.6) for i = 4 and i = 6

Now we repeat the procedure using i = 4 and i=6 in Eq.1.4

As the procedure is the same shown before, we develop a routine using Matlab to perform the products, the integration and to solve the system of equations. The routine is show in the appendix. The error was also Zero for i = 4, and i = 6.

Appendix

Matlab Code:

clear all;
x = sym('x');
b = [1 x+1 (x+1)^2];
db = [0 0 2];
k = b'*db;
K = 2*int(k,0,1);
f = -3*int(b',0,1);
 
K1 = [1 2 4; 0 1 2; 0 0 4];
f1 = [0 -4 -3]';
d = inv(K1)*f1;
z = [0:0.1:1];
 
for i = 1:11
 
    f2(i) = d(1)+ d(2)*(z(i)+1) + d(3)*(z(i)+1)^2; 
    %calculates Uh
    fnal(i) = -(3/4)*z(i)*z(i)-4*z(i)+19/4; %This is the analytic solution
end
figure, plot(z,f2,'ko') %plots both solution
hold on
plot(z,fnal,'k')
xlabel('x')
ylabel('U(x)')
legend('Numerical solution','Analytical solution')
title('Numerical and analytical solution with i = 2')
err(1) = abs(fnal(6)-f2(6));
 
%Now for n = 4
x = sym('x');
b = [1 x+1 (x+1)^2 (x+1)^3 (x+1)^4];
db = [0 0 2 6*(x+1) 12*(x+1)^2];
k = b'*db;
K2 = 2*int(k,0,1);
f = -3*int(b',0,1);
 
K1 = [1 2 4 8 16; 0 1 2 3 4; 0 0 4 18 56; 0 0 6 28 90; 0 0 28/3 45 744/5];
f1 = [0 -4 -3 -9/2 -7]';
d = inv(K1)*f1;
z = [0:0.1:1];
 
for i = 1:11
 
    f2(i) = d(1)+ d(2)*(z(i)+1) + d(3)*(z(i)+1)^2+d(4)*(z(i)+1)^3+d(5)*(z(i)+1)^4; 
    %calculates Uh
    fnal(i) = -(3/4)*z(i)*z(i)-4*z(i)+19/4; %This is the analytic solution
end
figure, plot(z,f2,'ko') %plots both solution
hold on
plot(z,fnal,'k')
xlabel('x')
ylabel('U(x)')
legend('Numerical solution','Analytical solution')
title('Numerical and analytical solution with i = 2')
err(2) = abs(fnal(6)-f2(6));
 
 
%Now for n = 6
x = sym('x');
b = [1 x+1 (x+1)^2 (x+1)^3 (x+1)^4 (x+1)^5 (x+1)^6];
db = [0 0 2 6*(x+1) 12*(x+1)^2 20*(x+1)^3 30*(x+1)^4];
k = b'*db;
K2 = 2*int(k,0,1)
f = -3*int(b',0,1)
 
K1 = [1 2 4 8 16 36 64; 0 1 2 3 4 5 6; 0 0 4 18 56 150 372;
    0 0 6 28 90 248 630; 0 0 28/3 45 744/5 420 7620/7;
    0 0  15  372/5   252 5080/7   3825/2;
    0 0  124/5  126   3048/7  1275  10220/3];
f1 = [0 -4 -3 -9/2 -7 -45/4 -93/5]';
d = inv(K1)*f1
z = [0:0.1:1];
 
for i = 1:11
 
    f2(i) = d(1)+ d(2)*(z(i)+1) + d(3)*(z(i)+1)^2+d(4)*(z(i)+1)^3+d(5)*(z(i)+1)^4 +d(6)*(z(i)+1)^5+d(7)*(z(i)+1)^6; 
    %calculates Uh
    fnal(i) = -(3/4)*z(i)*z(i)-4*z(i)+19/4; %This is the analytic solution
end
figure, plot(z,f2,'ko') %plots both solution
hold on
plot(z,fnal,'k')
xlabel('x')
ylabel('U(x)')
legend('Numerical solution','Analytical solution')
title('Numerical and analytical solution with i = 2')
err(3) = abs(fnal(6)-f2(6));

Problem 2: Spring System

Given

Spring System with closed ends.

Find

2.1 Number the Element and Nodes.

2.2 Assemble the Global Stiffness and Force Matrix.

2.3 Partition the system and solve for the Nodal Displacements.

2.4 Compute the Reaction Forces.

Solution

2.1 Number the Element and Nodes.

The Spring Systems is shown in the figure.

Fig 2.1 Spring System with Element & Node numbers specified

Superscripts to each of the stiffness values denotes the element number. Furthermore, Node Numbers are also shown in the figure. As shown in the figure, there are 4 nodes.

2.2 Assemble the Global Stiffness and Force Matrix.

For Element Number 1,$I=1$ and $J=4$,

$K^{(1)}=k\begin{bmatrix} 1& -1\\ -1 & 1 \end{bmatrix}\Rightarrow \tilde{K}=\begin{bmatrix} k& 0 & 0 & -k\\ 0 & 0 & 0 & 0\\ 0& 0 & 0 & 0\\ -k& 0 & 0 & k \end{bmatrix}$

For Element Number 2,$I=4$ and $J=3$,

$K^{(2)}=2k\begin{bmatrix} 1& -1\\ -1 & 1 \end{bmatrix}\Rightarrow \tilde{K}=\begin{bmatrix} 0& 0 & 0 & 0\\ 0 & 0 & 0 & 0\\ 0& 0 & 2k& -2k\\ 0& 0 & -2k & 2k \end{bmatrix}$

For Element Number 3,$I=3$ and $J=2$,

$K^{(3)}=k\begin{bmatrix} 1& -1\\ -1 & 1 \end{bmatrix}\Rightarrow \tilde{K}=\begin{bmatrix} 0& 0 & 0 & 0\\ 0 & k & -k & 0\\ 0& -k & k& 0\\ 0& 0 & 0& 0 \end{bmatrix}$

For Element Number 4,$I=1$ and $J=3$,

$K^{(4)}=3k\begin{bmatrix} 1& -1\\ -1 & 1 \end{bmatrix}\Rightarrow \tilde{K}=\begin{bmatrix} 3k& 0 & -3k & 0\\ 0 & 0& 0& 0\\ -3k& 0 & 3k& 0\\ 0& 0 & 0& 0 \end{bmatrix}$

Assembled Global Systems Matrix will be given by,

$K=\sum_{e=1}^{4}\tilde{K}^{e}=\begin{bmatrix} k+3k & 0 & -3k & -k\\ 0 & k & -k & 0\\ -3k & -k &6k &-2k \\ -k & 0 & -2k & 3k \end{bmatrix}$

Where Superscript $'e'$ denotes the Element Number.

The Displacement and Global Force Matrices for this system are given by,

$d=\begin{bmatrix} 0\\ 0\\ u_{3}\\ u_{4} \end{bmatrix}, f= \begin{bmatrix} 0\\ 0\\ 50\\ 0 \end{bmatrix}, r= \begin{bmatrix} r_{1}\\ r_{2}\\ 0\\ 0 \end{bmatrix}$

2.3 Partition the system and solve for the Nodal Displacements.

Eventually, The Global System of Equations is Given by,

$\begin{bmatrix} k+3k & 0 & -3k & -k\\ 0 & k & -k & 0\\ -3k & -k &6k &-2k \\ -k & 0 & -2k & 3k \end{bmatrix}\begin{bmatrix} 0\\ 0\\ u_{3}\\ u_{4} \end{bmatrix}=\begin{bmatrix} r_{1}\\ r_{2}\\ 50\\ 0 \end{bmatrix}$

$\displaystyle (Eq.2.1)$

As the first two displacements are prescribed, We can Partition the system after 2 Rows and 2 Columns,

And Can form the Matrix System of Equivalent to,

$\begin{bmatrix} K_{E} & K_{EF}\\ K_{EF}^{T} & K_{F} \end{bmatrix}\begin{bmatrix} \bar{d}_{E}\\ d_{F} \end{bmatrix}=\begin{bmatrix} r_{E}\\ f_{F} \end{bmatrix}$

Where,

$K_{E}=\begin{bmatrix} 4k & 0\\ 0 & K \end{bmatrix}, K_{EF}=\begin{bmatrix} -3k & -k\\ -k& 0 \end{bmatrix},K_{EF} ^T=\begin{bmatrix} -3k & -k\\ -k & 0 \end{bmatrix}, K_{F}=\begin{bmatrix} 6k & -2k\\ -2k & 3K \end{bmatrix}, \bar{d}_{E}=\begin{bmatrix} 0\\ 0\\ \end{bmatrix}, d_{F}=\begin{bmatrix} u_{3}\\ u_{4} \end{bmatrix},r_{E}=\begin{bmatrix} r_{1}\\ r_{2} \end{bmatrix}, f_{F}=\begin{bmatrix} 50\\ 0 \end{bmatrix}$

Solving the Following Reduced System of Equations,

$K_{E}\bar{d}_{E}+K_{EF}d_{F}=r_{E}, K_{EF}^T\bar{d}_{E}+K_{F}d_{F}=f_{F}$

Using $\displaystyle (Eq.2.1)$ we can form the Following Simultaneous Equations,

$6k.u_{3}-2k.u_{4}=50$

$\displaystyle (Eq.2.2)$

$-2k.u_{3}-3k.u_{4}=0$

$\displaystyle (Eq.2.3)$

Solving These 2 Simultaneous Equations, We get,

$u_{3}=\frac{75}{7k}\, and \, u_{4}=\frac{50}{7k}$.

2.4 Compute the Reaction Forces.

Now, Utilizing the First row of the math>\displaystyle (Eq.2.1) </math> we can get

$r_{1}=\frac{-275}{7}N$

And Further, Utilizing the Second row of the math>\displaystyle (Eq.2.1) </math> we can get,

$r_{2}=\frac{-75}{7}N$

Problem 3: Truss Structure

Given

The truss Structure Shown in the Figure. Nodes A & B are fixed. A Force of 10N Acts in the Positive x-Direction at node C. Young's Modulus is $E=10^{11}Pa$ and the Cross Sectional Area for all bars are $A= 2.10^{-2}$

Find

3.1 Number the Element and Nodes

3.2 Assemble the Global Stiffness and Force Matrix

3.3 Partition the system and solve for the Nodal Displacements

3.4 Compute the Stresses & Reactions

Solution

3.1 Number the Element and Nodes

The Truss System Along with the Element Numbers and the Node Numbers Written on it, is Shown in the figure.

Fig 3.1 Truss System with Element & Node numbers specified

3.2 Assemble the Global Stiffness and Force Matrix

For Element Number 1,$I=1$ and $J=3$, Also Element 1 is inclined at 90 deg to the Positive X-direction,$\phi^{(1)} =90^{0}$, $cos(90)=0,sin(90)=1$, $l=1$,

Using the Following formula for obtaining the Elemental Stiffness Matrix for all the elements,

$K^{e}=\frac{AE}{l}\begin{bmatrix} (cos\phi ^{e})^{2} & cos\phi ^{e}.sin\phi ^{e} & -(cos\phi ^{e})^{2} & -cos\phi ^{e}.sin\phi ^{e} \\ cos\phi ^{e}.sin\phi ^{e} & (sin\phi ^{e})^{2} & -cos\phi ^{e}.sin\phi ^{e} &-(sin\phi ^{e})^{2} \\ -(cos\phi ^{e})^{2} & -cos\phi ^{e}.sin\phi ^{e} & (cos\phi ^{e})^{2} & cos\phi ^{e}.sin\phi ^{e}\\ -cos\phi ^{e}.sin\phi ^{e} & -(sin\phi ^{e})^{2} & cos\phi ^{e}.sin\phi ^{e} & (sin\phi ^{e})^{2} \end{bmatrix}$

$\displaystyle (Eq.3.1)$

Utilizing above equation for Element 1,

$K^{(1)}=\frac{AE}{l}\begin{bmatrix} 0 & 0 & 0 & 0\\ 0 &1 & 0&-1 \\ 0 & 0 & 0 & 0\\ 0 & -1 & 0 & 1 \end{bmatrix}$

For Element Number 2,$I=2$ and $J=4$, Also Element 2 is inclined at 90 deg to the Positive X-direction,$\phi^{(2)} =90^{0}$, $cos(90)=0,sin(90)=1$, $l=1$, Therefore,

$K^{(2)}=\frac{AE}{l}\begin{bmatrix} 0 & 0 & 0 & 0\\ 0 &1 & 0&-1 \\ 0 & 0 & 0 & 0\\ 0 & -1 & 0 & 1 \end{bmatrix}$

For Element Number 3,$I=2$ and $J=3$, Also Element 3 is inclined at 135 deg to the Positive X-direction,$\phi^{(3)} =135^{0}$, $cos(135)=\frac{-1}{\sqrt{2}},sin(135)=\frac{1}{\sqrt{2}}$, $l=\sqrt{2}$, Therefore,

$K^{(3)}=\frac{AE}{l}\begin{bmatrix} \frac{1}{2\sqrt{2}} & -\frac{1}{2\sqrt{2}} & -\frac{1}{2\sqrt{2}} & \frac{1}{2\sqrt{2}}\\ -\frac{1}{2\sqrt{2}} &\frac{1}{2\sqrt{2}} & \frac{1}{2\sqrt{2}}&-\frac{1}{2\sqrt{2}} \\ -\frac{1}{2\sqrt{2}} & \frac{1}{2\sqrt{2}} & \frac{1}{2\sqrt{2}} & -\frac{1}{2\sqrt{2}}\\ \frac{1}{2\sqrt{2}} & -\frac{1}{2\sqrt{2}} & -\frac{1}{2\sqrt{2}} & \frac{1}{2\sqrt{2}} \end{bmatrix}$

For Element Number 4,$I=3$ and $J=4$, Also Element 4 is inclined at 90 deg to the Positive X-direction,$\phi^{(2)} =0^{0}$, $cos(0)=1,sin(0)=0$, $l=1$, Therefore,

$K^{(4)}=\frac{AE}{l}\begin{bmatrix} 1 & 0 & -1 & 0\\ 0 &0 & 0&0\\ -1 & 0 & 1 & 0\\ 0 & 0 & 0 & 0 \end{bmatrix}$

Making the Assembly of all these 4 Elemental Matrices, we get the following Global Stiffness Mtria for the given Truss System,

$K=AE \begin{bmatrix} 0 & 0 &0 & 0 & 0 & 0 & 0& 0\\ 0 & 1 & 0 & 0 & 0 & -1 & 0 &0 \\ 0 & 0 & \frac{1}{2\sqrt{2}} & \frac{-1}{2\sqrt{2}} & -\frac{1}{2\sqrt{2}} & \frac{1}{2\sqrt{2}} &0 &0 \\ 0 & 0 & -\frac{1}{2\sqrt{2}} &1+\frac{1}{2\sqrt{2}} & \frac{1}{2\sqrt{2}} & -\frac{1}{2\sqrt{2}} & 0 & -1\\ 0 & 0 & -\frac{1}{2\sqrt{2}} & \frac{1}{2\sqrt{2}} & 1+\frac{1}{2\sqrt{2}} &-\frac{1}{2\sqrt{2}} & -1 &0 \\ 0 &-1 &\frac{1}{2\sqrt{2}} & -\frac{1}{2\sqrt{2}} & -\frac{1}{2\sqrt{2}} &1+\frac{1}{2\sqrt{2}} & 0 & 0\\ 0 & 0 & 0 & 0 & -1 & 0 & 1 & 0\\ 0 & 0 & 0 & -1 & 0 & 0 & 0 & 1 \end{bmatrix}$

$\displaystyle (Eq.3.2)$

Now, Only force is acting in positive X-Direction at C, and there will be Reactions at A & B.

Therefore, The Force Matrix Can be Assembled as,

$F=\begin{bmatrix} r_{1x}\\ r_{1y}\\ r_{2x}\\ r_{2y}\\ 0\\ 0\\ 10\\ 0 \end{bmatrix}$

$\displaystyle (Eq.3.3)$

3.3 Partition the system and solve for the Nodal Displacements

From $\displaystyle (Eq.3.2)$ & $\displaystyle (Eq.3.3)$ The Global System os Equation will be,

$AE \begin{bmatrix} 0 & 0 &0 & 0 & 0 & 0 & 0& 0\\ 0 & 1 & 0 & 0 & 0 & -1 & 0 &0 \\ 0 & 0 & \frac{1}{2\sqrt{2}} & \frac{-1}{2\sqrt{2}} & -\frac{1}{2\sqrt{2}} & \frac{1}{2\sqrt{2}} &0 &0 \\ 0 & 0 & -\frac{1}{2\sqrt{2}} &1+\frac{1}{2\sqrt{2}} & \frac{1}{2\sqrt{2}} & -\frac{1}{2\sqrt{2}} & 0 & -1\\ 0 & 0 & -\frac{1}{2\sqrt{2}} & \frac{1}{2\sqrt{2}} & 1+\frac{1}{2\sqrt{2}} &-\frac{1}{2\sqrt{2}} & -1 &0 \\ 0 &-1 &\frac{1}{2\sqrt{2}} & -\frac{1}{2\sqrt{2}} & -\frac{1}{2\sqrt{2}} &1+\frac{1}{2\sqrt{2}} & 0 & 0\\ 0 & 0 & 0 & 0 & -1 & 0 & 1 & 0\\ 0 & 0 & 0 & -1 & 0 & 0 & 0 & 1 \end{bmatrix}.\begin{bmatrix} 0\\ 0\\ 0\\ 0\\ u_{3x}\\ u_{3y}\\ u_{4x}\\ u_{4y} \end{bmatrix}=\begin{bmatrix} r_{1x}\\ r_{1y}\\ r_{2x}\\ r_{2y}\\ 0\\ 0\\ 10\\ 0 \end{bmatrix}$

As the Displacements are prescribed untill second node, Partitioning the systems after second row and column and further comparing with the following equation,

$\begin{bmatrix} K_{E} & K_{EF}\\ K_{EF}^{T} & K_{F} \end{bmatrix}\begin{bmatrix} \bar{d}_{E}\\ d_{F} \end{bmatrix}=\begin{bmatrix} r_{E}\\ f_{F} \end{bmatrix}$

we get, $\bar{d_{E}}=\begin{bmatrix} 0\\ 0\\ 0\\ 0 \end{bmatrix},\ d_{F}=\begin{bmatrix} u_{3x}\\ u_{3y}\\ u_{4x}\\ u_{4y} \end{bmatrix},\ r_{E}=\begin{bmatrix} r_{1x}\\ r_{1y}\\ r_{2x}\\ r_{2y} \end{bmatrix},\ f_{F}=\begin{bmatrix} 0\\ 0\\ 10\\ 0 \end{bmatrix},\ K_{F}=\begin{bmatrix} 1+\frac{1}{2\sqrt{2}} & -\frac{1}{2\sqrt{2}} & -1 &0 \\ -\frac{1}{2\sqrt{2}} &1+\frac{1}{2\sqrt{2}} &0 & 0\\ -1 & 0 & 1 & 0\\ 0 & 0 & 0 &1 \end{bmatrix},\ K_{EF}=\begin{bmatrix} 0 & 0 & 0 & 0\\ 0 &-1 &0 & 0\\ -\frac{1}{2\sqrt{2}} & \frac{1}{2\sqrt{2}} & 0& 0\\ \frac{1}{2\sqrt{2}} & -\frac{1}{2\sqrt{2}} &0 &-1 \end{bmatrix}$

Eventually, we can obtain,

$d_{F}=K_{F}^{-1}(f_{F}-K_{EF}^T.\bar{d}_{E})$

But $\bar{d}_{E}$ is a Zero Matrix, which will yield, $d_{F}=K_{F}^{-1}.f_{F}$ Solving this system in MATLAB, Using the following code, we get,

Matlab Code:

kf=[1+1/(2*sqrt(2)) -1/(2*sqrt(2)) -1 0;-1/(2*sqrt(2)) 1+1/(2*sqrt(2)) 0 0;-1 0 1 0;0 0 0 1]
f=[0;0;10;0]   %force matrix
A=inv(kf)*f

$d_{F}=\begin{bmatrix} u_{3x}\\ u_{3y}\\ u_{4x}\\ u_{4y} \end{bmatrix}=\frac{1}{2\times 10^{9}}\begin{bmatrix} 38.2843\\ 10\\ 48.2843\\ 0 \end{bmatrix}$

$\displaystyle (Eq.3.4)$

These are the Unknown Nodal Displacements.

3.4 Compute the Stresses & Reactions

Now the Reactions can be calculated as,

$r_{E}=K_{E}.\bar{d}_{E}+K_{EF}.d_{F}$

But again $\bar{d}_{E}$ is a Zero Matrix, Which will give, $r_{E}=K_{EF}.d_{F}$.

Using the Following MATLAB Code to obtain the Reaction Matrix,

Matlab Code:

kef=[0 0 0 0;0 -1 0 0;-1/(2*sqrt(2)) 1/(2*sqrt(2)) 0 0;1/(2*sqrt(2)) -1/(2*sqrt(2)) 0 -1]
df=[38.2843;10;48.2843;0]
re=kef*df    %Reaction Matrix

$r_{E}=\begin{bmatrix} r_{1x}\\ r_{1y}\\ r_{2x}\\ r_{2y} \end{bmatrix}=\begin{bmatrix} 0\\ -10\\ -10\\ 10 \end{bmatrix}$

$\displaystyle (Eq.3.5)$

These are the Unknown Reactions. Now We can calculate the Stress in each element using the following formula,

$\sigma ^e=\frac{E^e}{l^e}\begin{bmatrix} -cos\phi ^{e} &-sin\phi ^{e} & cos\phi ^{e} & sin\phi ^{e} \end{bmatrix}.d^e$

$\displaystyle (Eq.3.6)$

Now we can utilize the above equation to calculate the Stresses for all the elements,

For element 1, $\phi =90^0, l=1$ &

$d^1=\begin{bmatrix} u_{1x}\\ u_{1y}\\ u_{3x}\\ u_{3y} \end{bmatrix}=\frac{1}{2\times 10^{9}}\begin{bmatrix} 0\\ 0\\ 38.2843\\ 10 \end{bmatrix}$

Using the Following MATLAB code, Matlab Code:

%For element e=1 to 4 calculation of Stress
x=[0 -1e11 0 1e11]
y=5e-10*[0;0;38.2843;10]
z=x*y

Therefore,

$\sigma ^1=\frac{E}{1}\begin{bmatrix} 0 &-1 &0 & 1 \end{bmatrix}.\frac{1}{2\times 10^{9}}\begin{bmatrix} 0\\ 0\\ 38.2843\\ 10 \end{bmatrix}=500 Pa$

In the similar fashion, For Element 2, Utilizing $\displaystyle (Eq.3.6)$ we get,

$\sigma ^2=\frac{E}{1}\begin{bmatrix} 0 &-1 &0 & 1 \end{bmatrix}.\frac{1}{2\times 10^{9}}\begin{bmatrix} 0\\ 0\\ 48.2843\\ 0 \end{bmatrix}=0$

For Element 3, Utilizing $\displaystyle (Eq.3.6)$ we get,

$\sigma ^3=\frac{E}{\sqrt{2}}\begin{bmatrix} \frac{1}{\sqrt{2}} &-\frac{1}{\sqrt{2}} &-\frac{1}{\sqrt{2}} & \frac{1}{\sqrt{2}} \end{bmatrix}.\frac{1}{2\times 10^{9}}\begin{bmatrix} 0\\ 0\\ 38.2843\\ 10 \end{bmatrix}=-1000 Pa$

Finally, For Element 4, Utilizing $\displaystyle (Eq.3.6)$ we get,

$\sigma^{4}=\frac{E}{1}\begin{bmatrix} -1 & 0 & 1& 0 \end{bmatrix}.\frac{1}{2\times 10^{9}}\begin{bmatrix} 38.2843\\ 10\\ 48.2843\\ 0 \end{bmatrix}=500 Pa$

Problem 4 : Show that $\displaystyle \alpha$ does not equal $\displaystyle \beta$

Given

Given

$\alpha =b_{i}\left[ a_{2}^{\prime }b_{j}^{\prime }+a_{2}b_{j}^{\prime \prime } \right]$

$\displaystyle (Eq. 3.1)$

And

$\beta =b_{j}\left[ a_{2}^{\prime }b_{i}^{\prime }+a_{2}b_{i}^{\prime \prime } \right]$

$\displaystyle (Eq. 3.2)$

Find

Show that $\alpha \ne \beta$

Solution

Let

$b_{i}\left( x \right)=\cos \left( ix \right)$

$\displaystyle (Eq. 3.3)$

And

$b_{j}\left( x \right)=\cos \left( jx \right)$

$\displaystyle (Eq. 3.4)$

Assume that $i=1$ and $j=2$. Therefore, equations 3.3 and 3.4 become:

$b_{i}\left( x \right)=\cos \left( x \right)$

$b_{j}\left( x \right)=\cos \left( 2x \right)$

Differentiating:

$b_{i}^{\prime }\left( x \right)=-\sin \left( x \right)$

$b_{j}^{\prime }\left( x \right)=-2\sin \left( 2x \right)$

$b_{i}^{\prime \prime }\left( x \right)=-\cos \left( x \right)$

$b_{j}^{\prime \prime }\left( x \right)=-4\cos \left( 2x \right)$

Plugging in values and re-arranging

$\alpha =-2a_{2}^{\prime }\sin \left( 2x \right)\cos \left( x \right)-4a_{2}^{\prime \prime }\cos \left( 2x \right)\cos \left( x \right)$

$\displaystyle (Eq. 3.5)$

$\beta =-a_{2}^{\prime }\sin \left( x \right)\cos \left( 2x \right)-a_{2}^{\prime \prime }\cos \left( x \right)\cos \left( 2x \right)$

$\displaystyle (Eq. 3.6)$

By inspection of terms $\displaystyle a_{2}^{\prime }$ and $\displaystyle a_{2}^{\prime \prime }$ it can be shown that

$\alpha \ne \beta$

$\displaystyle ( Eq. 3.7)$

Problem 5: Equivalent Stiffness of a Spring

Given

A spring is given in Figure 5.1 as follows:

Fig.5.1 Data for Problem5

Find

Show that the equivalent stiffness of a spring aligned in the x direction for the bar of thickness t with a centered square hole shown in Figure 5.1 is can be described by

$k=\frac{5Etab}{(a+b)l}$

$\displaystyle (Eq. 5.1)$

Solution

The left end is constrained, i.e., prescribed displacement is zero at the left end.

Assume there is a force F acting on the right end, and the nodes are numbered in (Fig.5.2) .

Fig.5.2 Four Nodes and the Force

The Element stiffness matrices are the following:

$\mathbf{K}^{\left ( 1\right )} =\begin{bmatrix} k^{\left ( 1 \right )} & -k^{\left ( 1 \right )} \\ -k^{\left ( 1 \right )} & k^{\left ( 1 \right )} \end{bmatrix}, \mathbf{K}^{\left ( 2\right )} =\begin{bmatrix} k^{\left ( 2 \right )} & -k^{\left ( 2 \right )} \\ -k^{\left ( 2 \right )} & k^{\left ( 2 \right )} \end{bmatrix}, \mathbf{K}^{\left ( 3\right )} =\begin{bmatrix} k^{\left ( 3 \right )} & -k^{\left ( 3 \right )} \\ -k^{\left ( 3 \right )} & k^{\left ( 3 \right )} \end{bmatrix}$

$\displaystyle (Eq. 5.2)$

By direct assembly, the global stiffness matrix is:

$\mathbf{K}=\begin{bmatrix} k^{\left ( 1 \right )} & -k^{\left ( 1 \right )} & 0 & 0\\ -k^{\left ( 1 \right )} & k^{\left ( 1 \right )} + k^{\left ( 2 \right )} & -k^{\left ( 2 \right )} & 0\\ 0 & -k^{\left ( 2 \right )} & k^{\left ( 2 \right )} + k^{\left ( 3 \right )} & -k^{\left ( 3 \right )} \\ 0 & 0 & -k^{\left ( 3 \right )} & k^{\left ( 3 \right )} \end{bmatrix}$

$\displaystyle (Eq. 5.3)$

The displacements and force matrices for the system are:

$\mathbf{d}=\begin{bmatrix} 0\\ u_{2}\\ u_{3}\\ u_{4} \end{bmatrix}, \mathbf{f}=\begin{bmatrix} 0\\ 0\\ 0\\ F \end{bmatrix}, \mathbf{r}=\begin{bmatrix} r_{1}\\ 0\\ 0\\ 0\end{bmatrix}$

$\displaystyle (Eq. 5.4)$

The global system of equations is given by:

$\begin{bmatrix} k^{\left ( 1 \right )} & -k^{\left ( 1 \right )} & 0 & 0\\ -k^{\left ( 1 \right )} & k^{\left ( 1 \right )} + k^{\left ( 2 \right )} & -k^{\left ( 2 \right )} & 0\\ 0 & -k^{\left ( 2 \right )} & k^{\left ( 2 \right )} + k^{\left ( 3 \right )} & -k^{\left ( 3 \right )} \\ 0 & 0 & -k^{\left ( 3 \right )} & k^{\left ( 3 \right )} \end{bmatrix} \begin{bmatrix} 0\\ u_{2}\\ u_{3}\\ u_{4} \end{bmatrix}= \begin{bmatrix} r_{1}\\ 0\\ 0\\ F \end{bmatrix}$

$\displaystyle (Eq. 5.5)$

Solve the matrix equation above; we can have the following four equations

$\left\{\begin{matrix} -k^{\left ( 1 \right )}u_{2}=r_{1} \\ u_{2}\left ( k^{\left (1 \right )} + k^{\left (2 \right)} \right )- u_{3}k^{\left (2 \right)}=0 \\ -u_{2}k^{\left (2 \right )} + u_{3}\left ( k^{\left (2 \right )} + k^{\left (3 \right)} \right ) - u_{4}k^{\left (3 \right)}=0 \\ -u_{3}k^{\left (3 \right)}+u_{4}k^{\left (3 \right)}=F \end{matrix}\right.$

$\displaystyle (Eq. 5.6)$

Because we have the following relationship:

$k=\frac{AE}{l}$

$\displaystyle (Eq. 5.7)$

Then we can get the following results:

$k^{\left ( 1 \right )}=\frac{10atE}{l}, k^{\left ( 2 \right )}=\frac{5btE}{l}, k^{\left ( 3 \right )}=\frac{10atE}{l}$

$\displaystyle (Eq. 5.8)$

Substitute (Eq.5.8) into (Eq.5.6), we can have the following result:

$u_{4}=\frac{Fl\left ( a+b \right )}{5abtE}$

$\displaystyle (Eq. 5.9)$

where $u_{4}$ is the displacement of node 4.

And then,

$\frac{F}{u_{4}}=\frac{5Etab}{\left ( a+b \right )l}$

$\displaystyle (Eq. 5.10)$

which yields:

$k=\frac{5Etab}{\left ( a+b \right )l}$

$\displaystyle ( Eq. 5.11)$

Problem 6: Three Bar Structure

Given

Given a three bar structure subjected to the prescribed load at point C equal to $10^3 ~N$ as shown in Figure 2.19 of the book. The Young's modulus is $E = 10^{11} ~Pa$, the cross-sectional area of the bar BC is $2 \times 10^{-2} ~m^2$ and that of BD and BF is $10^{-2} ~m^2$. Note that point D is free to move in the x-direction. Coordinates of joints are given in meters.

Fig.6.1 Truss System with Nodes and Elements Noted

Find

6.1 Construct the global stiffness matrix and load matrix.

6.2 Partition the matrices and solve for the unknown displacements at Point B and displacement in the x-direction at point D.

6.3 Find stresses in the three bars.

6.4 Find the reactions at nodes C,D,F.

Solution

6.1 Construct the global stiffness matrix and load matrix

For Element 3, $\phi = 45^o,~~I = 3,~~J =4$

$K^{\left(3\right)} = \frac{10^{-2} \times 10^{11}}{\sqrt{2}} \begin{bmatrix} \frac{1}{2} &\frac{1}{2} &-\frac{1}{2} &-\frac{1}{2}\\ \frac{1}{2} &\frac{1}{2} &-\frac{1}{2} &-\frac{1}{2}\\ -\frac{1}{2} &-\frac{1}{2} &\frac{1}{2} &\frac{1}{2}\\ -\frac{1}{2} &-\frac{1}{2} &\frac{1}{2} &\frac{1}{2} \end{bmatrix}$

For Element 2, $\phi = 90^o,~~I = 2,~~J = 4$

$K^{\left(2\right)} = \frac{2 \times 10^{-2} \times 10^{11}}{1} \begin{bmatrix} 0 &0 &0 &0\\ 0 &1 &0 &-1\\ 0 &0 &0 &0 \\ 0 &-1 &0 &1 \end{bmatrix}$

For Element 1, $\phi = 135^o,~~I = 1,~~J =4$

$K^{\left(1\right)} = \frac{10^{-2} \times 10^{11}}{\sqrt{2}} \begin{bmatrix} \frac{1}{2} &-\frac{1}{2} &-\frac{1}{2} &\frac{1}{2}\\ -\frac{1}{2} &\frac{1}{2} &\frac{1}{2} &-\frac{1}{2}\\ -\frac{1}{2} &\frac{1}{2} &\frac{1}{2} &-\frac{1}{2}\\ \frac{1}{2} &-\frac{1}{2} &-\frac{1}{2} &\frac{1}{2} \end{bmatrix}$

Global Stiffness Matrix

$K = 10^{9} \begin{bmatrix} \frac{1}{2 \sqrt{2}} &-\frac{1}{2 \sqrt{2}} &0 &0 &0 &0 &-\frac{1}{2 \sqrt{2}} &\frac{1}{2 \sqrt{2}} \\ -\frac{1}{2 \sqrt{2}} &\frac{1}{2 \sqrt{2}} &0 &0 &0 &0 &\frac{1}{2 \sqrt{2}} &-\frac{1}{2 \sqrt{2}} \\ 0 &0 &0 &0 &0 &0 &0 &0 \\ 0 &0 &0 &2 &0 &0 &0 &-2 \\ 0 &0 &0 &0 &\frac{1}{2 \sqrt{2}} &\frac{1}{2 \sqrt{2}} &-\frac{1}{2 \sqrt{2}} &-\frac{1}{2 \sqrt{2}} \\ 0 &0 &0 &0 &\frac{1}{2 \sqrt{2}} &\frac{1}{2 \sqrt{2}} &-\frac{1}{2 \sqrt{2}} &-\frac{1}{2 \sqrt{2}} \\ -\frac{1}{2 \sqrt{2}} &\frac{1}{2 \sqrt{2}} &0 &0 &-\frac{1}{2 \sqrt{2}} &-\frac{1}{2 \sqrt{2}} &\frac{1}{\sqrt{2}} &0 \\ \frac{1}{2 \sqrt{2}} &-\frac{1}{2 \sqrt{2}} &0 &-2 &-\frac{1}{2 \sqrt{2}} &-\frac{1}{2 \sqrt{2}} &0 &\frac{1}{\sqrt{2}}+2 \end{bmatrix}$

$d = \begin{bmatrix} 0\\ 0\\ 0\\ 0\\ u_{3x}\\ 0\\ u_{4x}\\ u_{4y} \end{bmatrix}~~ f = \begin{bmatrix} 0\\ 0\\ 0\\ 0\\ 0\\ 0\\ 10^3\\ 0 \end{bmatrix}~~ r = \begin{bmatrix} r_{1x}\\ r_{1y}\\ r_{2x}\\ r_{2y}\\ 0\\ r_{3y}\\ 0\\ 0 \end{bmatrix}$

Therefore

$load~matrix = \begin{bmatrix} r_{1x}\\ r_{1y}\\ r_{2x}\\ r_{2y}\\ 0\\ r_{3y}\\ 10^3\\ 0 \end{bmatrix}$

6.2 Partition the matrices and solve

Global system of equations:

$10^9 \begin{bmatrix} \frac{1}{2 \sqrt{2}} &-\frac{1}{2 \sqrt{2}} &0 &0 &\mid &0 &0 &-\frac{1}{2 \sqrt{2}} &\frac{1}{2 \sqrt{2}} \\ -\frac{1}{2 \sqrt{2}} &\frac{1}{2 \sqrt{2}} &0 &0 &\mid &0 &0 &\frac{1}{2 \sqrt{2}} &-\frac{1}{2 \sqrt{2}} \\ 0 &0 &0 &0 &\mid &0 &0 &0 &0 \\ 0 &0 &0 &2 &\mid &0 &0 &0 &-2 \\ - &- &- &- &\mid &- &- &- &- \\ 0 &0 &0 &0 &\mid &\frac{1}{2 \sqrt{2}} &\frac{1}{2 \sqrt{2}} &-\frac{1}{2 \sqrt{2}} &-\frac{1}{2 \sqrt{2}} \\ 0 &0 &0 &0 &\mid &\frac{1}{2 \sqrt{2}} &\frac{1}{2 \sqrt{2}} &-\frac{1}{2 \sqrt{2}} &-\frac{1}{2 \sqrt{2}} \\ -\frac{1}{2 \sqrt{2}} &\frac{1}{2 \sqrt{2}} &0 &0 &\mid &-\frac{1}{2 \sqrt{2}} &-\frac{1}{2 \sqrt{2}} &\frac{1}{\sqrt{2}} &0 \\ \frac{1}{2 \sqrt{2}} &-\frac{1}{2 \sqrt{2}} &0 &-2 &\mid &-\frac{1}{2 \sqrt{2}} &-\frac{1}{2 \sqrt{2}} &0 &\frac{1}{\sqrt{2}}+2 \end{bmatrix} \begin{bmatrix} 0\\ 0\\ 0\\ 0\\ --\\ u_{3x}\\ 0\\ u_{4x}\\ u_{4y} \end{bmatrix} = \begin{bmatrix} r_{1x}\\ r_{1y}\\ r_{2x}\\ r_{2y}\\ --\\ 0\\ r_{3y}\\ 10^3\\ 0 \end{bmatrix}$

Partitions

$\bar{d}_{E} = \begin{bmatrix} u_{1x}\\ u_{1y}\\ u_{2x}\\ u_{2y} \end{bmatrix} = \begin{bmatrix} 0\\ 0\\ 0\\ 0 \end{bmatrix},~~ d_{F}= \begin{bmatrix} u_{3x}\\ 0\\ u_{4x}\\ u_{4y} \end{bmatrix},~~ r_{E} = \begin{bmatrix} r_{1x}\\ r_{1y}\\ r_{2x}\\ r_{2y} \end{bmatrix} ,~~ f_{F} = \begin{bmatrix} 0\\ r_{3y}\\ 10^3\\ 0 \end{bmatrix}$

$K_{EF} = \begin{bmatrix} 0 &0 &-\frac{1}{2 \sqrt{2}} &\frac{1}{2 \sqrt{2}} \\ 0 &0 &\frac{1}{2 \sqrt{2}} &-\frac{1}{2 \sqrt{2}} \\ 0 &0 &0 &0 \\ 0 &0 &0 &-2\ \end{bmatrix},~~~ K_{EF}^{T} = \begin{bmatrix} 0 &0 &0 &0 \\ 0 &0 &0 &0 \\ -\frac{1}{2 \sqrt{2}} &\frac{1}{2 \sqrt{2}} &0 &0 \\ \frac{1}{2 \sqrt{2}} &-\frac{1}{2 \sqrt{2}} &0 &-2 \end{bmatrix},~~~ K_{F} = \begin{bmatrix} \frac{1}{2 \sqrt{2}} &\frac{1}{2 \sqrt{2}} &-\frac{1}{2 \sqrt{2}} &-\frac{1}{2 \sqrt{2}} \\ \frac{1}{2 \sqrt{2}} &\frac{1}{2 \sqrt{2}} &-\frac{1}{2 \sqrt{2}} &-\frac{1}{2 \sqrt{2}} \\ -\frac{1}{2 \sqrt{2}} &-\frac{1}{2 \sqrt{2}} &\frac{1}{\sqrt{2}} &0 \\ -\frac{1}{2 \sqrt{2}} &-\frac{1}{2 \sqrt{2}} &0 &\frac{1}{\sqrt{2}}+2 \end{bmatrix},~~~ K_{E} = \begin{bmatrix} \frac{1}{2 \sqrt{2}} &-\frac{1}{2 \sqrt{2}} &0 &0 \\ -\frac{1}{2 \sqrt{2}} &\frac{1}{2 \sqrt{2}} &0 &0 \\ 0 &0 &0 &0 \\ 0 &0 &0 &2 \end{bmatrix}$

Unknown displacement at B is solved by

$K_{EF}^{T}\bar{d}_{E} + K_{F}d_{F} = f_{F}$

Subtracting the first term from both sides of the above equation and premultiply by $K_F^{-1}$, we obtain

$d_F = K_F^{-1} \left( f_F - K_{EF}^T \bar{d}_E \right)$

Solving for the displacements we get: $d_{F}= 10^{-6} \begin{bmatrix} 1+2\sqrt{2}\\ 0\\ \frac{1}{2}+2\sqrt{2}\\ \frac{1}{2} \end{bmatrix} ~m$

6.3 Find stresses in the three bars

Element 3 Stress

$\sigma = 0 ~Pa$

Element 2 Stress

$\sigma = \frac{E}{l} \begin{bmatrix} 0 &-1 &0 &1 \end{bmatrix} 10^{-6} \begin{bmatrix} \frac{1}{2} + 2\sqrt{2}\\ \frac{1}{2}\\ 0\\ 0 \end{bmatrix} = -5 \times 10^4 ~Pa$

Element 1 Stress

$\sigma = \frac{E}{l} \begin{bmatrix} \frac{\sqrt{2}}{2} &-\frac{\sqrt{2}}{2} &\frac{\sqrt{2}}{2} &\frac{\sqrt{2}}{2} \end{bmatrix} 10^{-6} \begin{bmatrix} \frac{1}{2} + 2\sqrt{2}\\ \frac{1}{2}\\ 0\\ 0 \end{bmatrix} = \sqrt{2} \times 10^5 ~Pa$

6.4 Find the reactions at nodes C,D,F.

Reactions can be solved using this equation and the global matrix

$r_E = K_E \bar{d}_E + K_{EF} d_F$

$load~matrix = \begin{bmatrix} r_{1x}\\ r_{1y}\\ r_{2x}\\ r_{2y}\\ 0\\ r_{3y}\\ 10^3\\ 0 \end{bmatrix} = \begin{bmatrix} -10^3\\ 10^3\\ 0\\ -10^3\\ 0\\ 0\\ 10^3\\ 0 \end{bmatrix} ~ N$

Problem 7: Solving using WRF

Solving a differential equation using Weight Residual Form. From lecture 17

Given

The partial differential equation is:

$\frac{\partial}{\partial x}\left(2\frac{\partial u}{\partial x}\right) + 3 = 0$

$\displaystyle (Eq. 7.1)$

The boundary conditions are:

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

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

As weighting function should be used:

$b_{i} = b_{i}=\left ( x \right )^{i}$

$\displaystyle (Eq. 7.2)$

Find

7.1 Find an approximate solution of the form $U^{h}=\sum_{i=0}^{n}d_{i}b_{i}$ for n=2

7.2 Find two equations that enforce the boundary conditions

7.3 Project the weight residues

7.4 Display the equations in matrix form

7.5 Solve for $\displaystyle d$

7.6 Construct $U^{h}$ and plot $U^{h}$ and $u$

7.7 Repeat 7.1 to 7.6 for n = 4 and n = 6

Solution

7.1 Find an approximate solution with n=2

When n = 2 we have:

$d_{i}=\left[ \begin{matrix} d_{0} & d_{1} & d_{2} \\ \end{matrix} \right]$

$\displaystyle (Eq.7.3)$

and

$b_{i}=\left[ \begin{matrix} 1 & \left( x \right) & \left( x \right)^{2} \\ \end{matrix} \right]$

$\displaystyle (Eq.7.4)$

Therefore:

$U^{h}=d_{0}+d_{1}\left( x \right)+d_{2}\left( x \right)^{2}$

$\displaystyle (Eq. 1.5)$

7.2 Find two equations that enforce the boundary conditions

The boundary conditions are:

$u\left( 1 \right)=0=d_{0}+d_{1}+d_{2}$

$\displaystyle (Eq. 7.6)$

Differentiating

$\frac{du}{dx}=0+d_{1}+2d_{2}x$

And evaluating

$\frac{du}{dx}(0)=-4=d_{1}$

$\displaystyle (Eq. 7.7)$

7.3 Project the weight residues

In this case changing the function $u\left( x \right)$ by the approximated function $U^{h}$ in the PDE (Eq.7.1), we found that $P\left( u \right)$ is:

$P\left( U^{h} \right)=2\frac{\partial ^{2}U^{h}}{\partial x^{2}}+3=4d_{2}+3$

$\displaystyle (Eq.7.8)$

The (Eq.7.8) can be written as:

$P\left( U^{h} \right)=\left[ \begin{matrix} 0 & 0 & 4 \\ \end{matrix} \right]\left[ \begin{matrix} d_{0} \\ d_{1} \\ d_{2} \\ \end{matrix} \right]+3$

$\displaystyle (Eq.7.9)$

Projecting the residue we have:

$\int_{a}^{b}b(x)P\left ( U^{h} \right )= 0$

$\displaystyle (Eq. 7.10)$

After substitution of equations 7.4 and 7.9) in 7.10, the following equation is obtained

$\int_{0}^{1}{\left[ \begin{matrix} 1 \\ x \\ x^{2} \\ \end{matrix} \right]}\left[ \begin{matrix} 0 & 0 & 4 \\ \end{matrix} \right]dx\left[ \begin{matrix} d_{0} \\ d_{1} \\ d_{2} \\ \end{matrix} \right]=-3\int_{0}^{1}{\left[ \begin{matrix} 1 \\ x \\ x^{2} \\ \end{matrix} \right]}dx$

$\displaystyle (Eq. 7.11)$

7.4 Display the equations in matrix form

Performing the product and then the integration indicated in equation 7.11, we have

$\left[ \begin{matrix} 0 & 0 & 4 \\ 0 & 0 & 2 \\ 0 & 0 & \frac{4}{3} \\ \end{matrix} \right]\left[ \begin{matrix} d_{0} \\ d_{1} \\ d_{2} \\ \end{matrix} \right]=\left[ \begin{matrix} -3 \\ -\frac{3}{2} \\ -1 \\ \end{matrix} \right]$

$\displaystyle (Eq. 7.12)$

The latter equation is clearly of the form $\displaystyle Kd=F$. We can observe that matrix K is not symmetric.

7.5 Solve for $\displaystyle d$

In equation 7.12 we have three equations; but, as we need to enforce the boundary conditions we take the only one of them and solve it together with equations 7.6 and 7.7. From these equations we now have the system:

$\left[ \begin{matrix} 1 & 1 & 1 \\ 0 & 1 & 2 \\ 0 & 0 & 2 \\ \end{matrix} \right]\left[ \begin{matrix} d_{0} \\ d_{1} \\ d_{2} \\ \end{matrix} \right]=\left[ \begin{matrix} 0 \\ -4 \\ -3 \\ \end{matrix} \right]$

$\displaystyle (Eq.7.13)$

The latter is a linear system and can be solved by using any calculator. After solving the system in (Eq.7.13) we obtain d:

$\left[ \begin{matrix} d_{0} \\ d_{1} \\ d_{2} \\ \end{matrix} \right]=\left[ \begin{matrix} \frac{5}{2} \\ -1 \\ -\frac{4}{3} \\ \end{matrix} \right]$

7.6 Construct $U^{h}$ and plot $U^{h}$ and $u$

Replacing $d_{i}$ in equation 7.5 we obtain the approximated solution:

$U^{h}=\frac{5}{2}-x-\frac{4}{3}x^{2}$

$\displaystyle (Eq.7.14)$

On the other hand, the exact solution of the PDE (Eq.7.1) with the given boundary condition is:

$u\left ( x \right )=-\frac{3}{4}x^{2}-4x+\frac{19}{4}$

$\displaystyle (Eq.7.15)$

We observe that the approximated and the exact solution are not the same. Perhaps they would be closer if we increased the number of terms. Figure 7.1 we show the exact and the approximated solution.

Fig.1.1 Analytical and numerical solution with i=2

7.7 Repeat 7.1 to 7.6) for i = 4 and i = 6

Now we repeat the procedure using n = 4 and n=6 in Eq.7.4

As the procedure is the same shown before, we develop a routine using Matlab to perform the products, the integration and to solve the system of equations. The routine is shown below in the appendix. Figures 7.2 and 7.3 show the results of n=4 and n=6, respectively. Figure 7.4 shows the error versus the number of terms.

Fig.7.2 Analytical and numerical solution with n=4

Fig.7.3 Analytical and numerical solution with n=6

Fig.7.4 Error

Appendix

Matlab Code:

clear all;
x = sym('x');
b = [1 x (x)^2];
db = [0 0 2];
k = b'*db;
K = 2*int(k,0,1);
f = -3*int(b',0,1);
 
K1 = [1 1 1; 0 1 2; 0 0 2];
f1 = [0 -4 -3]';
d = inv(K1)*f1;
z = [0:0.1:1];
 
for i = 1:11
 
    f2(i) = d(1)+ d(2)*(z(i)) + d(3)*(z(i))^2; 
    %calculates Uh
    fnal(i) = -(3/4)*z(i)*z(i)-4*z(i)+19/4; %This is the analytic solution
end
figure, plot(z,f2,'ko') %plots both solution
hold on
plot(z,fnal,'k')
xlabel('x')
ylabel('U(x)')
legend('Numerical solution','Analytical solution')
title('Numerical and analytical solution with i = 2')
err(1) = abs(fnal(6)-f2(6));
 
 
x = sym('x');
b = [1 x (x)^2 (x)^3 (x)^4];
db = [0 0 2 6*(x) 12*(x)^2];
k = b'*db;
K2 = 2*int(k,0,1);
f = -3*int(b',0,1);
 
K1 = [1 1 1 1 1; 0 1 2 3 4; 0 0 2 6 12; 0 0 0 6 24; 0 0 0 0 24];
f1 = [0 -4 -3 -9/2 -7]';
d = inv(K1)*f1;
z = [0:0.1:1];
 
for i = 1:11
 
    f2(i) = d(1)+ d(2)*(z(i)) + d(3)*(z(i))^2+d(4)*(z(i))^3+d(5)*(z(i))^4; 
    %calculates Uh
    fnal(i) = -(3/4)*z(i)*z(i)-4*z(i)+19/4; %This is the analytic solution
end
figure, plot(z,f2,'ko') %plots both solution
hold on
plot(z,fnal,'k')
xlabel('x')
ylabel('U(x)')
legend('Numerical solution','Analytical solution')
title('Numerical and analytical solution with n = 4')
err(2) = abs(fnal(6)-f2(6));
 
 
 
%Now for n = 6
x = sym('x');
b = [1 x (x)^2 (x)^3 (x)^4 (x)^5 (x)^6];
db = [0 0 2 6*(x) 12*(x)^2 20*(x)^3 30*(x)^4];
k = b'*db;
K2 = 2*int(k,0,1);
f = -3*int(b',0,1);
 
K1 = [1 1 1 1 1 1 1; 0 1 2 3 4 5 6; 0 0 2 6 12 20 30; 0 0 0 6 24 60 120; 0 0 0 0 24 120 360; 0 0 0 0 0 120 720; 0 0 0 0 0 0 720];
f1 = [0 -4 -3 -9/2 -7 -45/4 -93/5]';
d = inv(K1)*f1
z = [0:0.1:1];
 
for i = 1:11
 
    f2(i) = d(1)+ d(2)*(z(i)) + d(3)*(z(i))^2+d(4)*(z(i))^3+d(5)*(z(i))^4 +d(6)*(z(i))^5+d(7)*(z(i))^6; 
    %calculates Uh
    fnal(i) = -(3/4)*z(i)*z(i)-4*z(i)+19/4; %This is the analytic solution
end
figure, plot(z,f2,'ko') %plots both solution
hold on
plot(z,fnal,'k')
xlabel('x')
ylabel('U(x)')
legend('Numerical solution','Analytical solution')
title('Numerical and analytical solution with n = 6')
err(3) = abs(fnal(6)-f2(6));
 n = [2 4 6]
 figure, plot(n,err)
 xlabel('n=number of functions')
 ylabel('Error')
 title('Difference between anlytical and numerical solution as function of i, at x = 0.5')

Problem 8: Weak Form Proof

Given

The Strong form is:

$\displaystyle \frac{\mathrm{d} }{\mathrm{d} x} \left ( AE\frac{\mathrm{d} u}{\mathrm{d} x} \right ) + 2x=0$

$\displaystyle (Eq. 8.1a)$

on 1<x<3,

$\displaystyle \sigma (1)=\left ( E\frac{\mathrm{d} u}{\mathrm{d} x} \right )_{x=1}=0.1,$

$\displaystyle (Eq. 8.1b)$

$\displaystyle u(3)=0.001$

$\displaystyle (Eq. 8.1c)$

Find

Show that the weak form is:

$\int_{1}^{3}\frac{\mathrm{d} w}{\mathrm{d} x}AE\frac{\mathrm{d} u}{\mathrm{d} x}dx = -0.1\left ( wA \right )_{x=1}+\int_{1}^{3}2xwdx,\forall w$ with $w(3)=0$

$\displaystyle (Eq. 8.2)$

Solution

Equation (8.1b) is a condition on the derivative of u(x), so it is a natural boundary conditon; (8.1c) is a condition on u(x), so it is an essential boundary condition. Therefore, as the weight function must vanish on the essential boundaries, we consider all smooth weight functions w(x) such that w(3)=0. The trial solutions must satisfy the essential boundary condition u(3) = 0.001.

We start by multiplying the governing equation and the natural boundary condition over the domains where they hold by an arbitrary weight function:

$\int_{1}^{3} w [\frac{\mathrm{d} }{\mathrm{d} x}(AE\frac{\mathrm{d} u}{\mathrm{d} x})+2x]dx = 0 , \forall w(x),$

$\displaystyle (Eq. 8.3)$

$[wA(E\frac{\mathrm{d} u}{\mathrm{d} x}-0.1)]_{x=1} = 0, \forall w(1).$

$\displaystyle (Eq. 8.4)$

Next we integrate (Eq.8.3) by parts, then we get the following result:

$\int_{1}^{3} [w (\frac{\mathrm{d} }{\mathrm{d} x}AE\frac{\mathrm{d} u}{\mathrm{d} x})]dx = \left (wAE\frac{\mathrm{d} u}{\mathrm{d} x} \right )\mid _{x=1}^{x=3} - \int_{1}^{3} \frac{\mathrm{d} w }{\mathrm{d} x}AE\frac{\mathrm{d} u}{\mathrm{d} x}dx$

$\displaystyle (Eq. 8.5)$

We have constructed the weight functions so that w(3)=0; therefore, the first term on the RHS of (Eq.8.5) vanishes at x=3. Substituting (Eq.8.5) into (Eq.8.3) gives:

$-\left (wAE\frac{\mathrm{d} u}{\mathrm{d} x} \right )\mid _{x=1} - \int_{1}^{3} \frac{\mathrm{d} w }{\mathrm{d} x}AE\frac{\mathrm{d} u}{\mathrm{d} x}dx + \int_{1}^{3} 2xwdx=0, \forall w(x)$

$\displaystyle (Eq. 8.6)$

with w(3)=0. Substituting (Eq.8.4) into the first term of (Eq.8.6) gives:

$\int_{1}^{3} \frac{\mathrm{d} w }{\mathrm{d} x}AE\frac{\mathrm{d} u}{\mathrm{d} x}dx = -0.1(wA)_{x=1} + \int_{1}^{3} 2xwdx, \forall w(x)$

$\displaystyle ( Eq. 8.7)$

with w(3)=0.

Problem 9: A solution to the weak form in Problem 3.1

Consider a trial (candidate) solution of the form $u\left ( x \right ) = \alpha_{0} + \alpha_{1}\left ( x-3 \right )$ and a weight function of the same form. Obtain a solution to the weak form in Problem 3.1. Check the equilibrium equation in the strong form in Problem 3.1; is it satisfied?

Check the natural boundary condition; is it satisfied?

Solution

$u\left ( x \right ) = \alpha_{0} + \alpha_{1}\left ( x-3 \right )$

$w\left ( x \right ) = \beta_{0} + \beta_{1}\left ( x-3 \right )$

Using boundary conditions from Problem 3.9

$w \left(3 \right) = 0 ~~ \therefore ~~ \beta_{0} = 0$

$u \left(3 \right) = \alpha_{0} = 0.001$

Only one unknown parameter and one arbitrary parameter remain

$u\left ( x \right ) = 0.001 + \alpha_{1}\left ( x-3 \right ) ~~~ \frac{du}{dx} = \alpha_{1}$

$w\left( x \right) = \beta_{1} \left(x - 3 \right) ~~~ \frac{dw}{dx} = \beta_{1}$

Now substitute into Problem 3.1 weak form

$\int_{1}^{3} \frac{\mathrm{d} w }{\mathrm{d} x}AE\frac{\mathrm{d} u}{\mathrm{d} x}dx = -0.1(wA)_{x=1} + \int_{1}^{3} 2xwdx$

$\int_{1}^{3} \beta_{1}\alpha_{1}EAdx + 0.1(\beta_{1} \left(x-3 \right)A)_{x=1} - \int_{1}^{3} 2\beta_{1}x \left(x-3 \right)dx = 0$

Integrating and factoring out $\beta_1$ we get:

$\beta_{1} \left( 2 \alpha_{1} EA - 0.2A + \frac{20}{3} \right) = 0$

The $\beta_1$ goes to zero and $\alpha_1$ can be solved:

$\alpha_{1} = \frac{0.1A - \frac{10}{3}}{EA}$

Therefore the solution to the weak form is

$u \left( x \right) = 0.001 + \frac{0.1A - \frac{10}{3}}{EA} \left( x-3 \right)$

Checking for equilibrium in the strong form, values are substituted into

$\frac{d}{dx} \left( AE \frac{du}{dx} \right) + 2x = 0$

This yields $2x = 0$ which does not satisfy.

Checking the natural boundary condition, this equation is used

$\sigma \left(1 \right) = \left(E \frac{du}{dx} \right)_{x=1} = 0.1$

Substituting in $\alpha_1$ we get

$0.1 - \frac{\frac{10}{3}}{A} = 0.1$

Thus, this boundary condition will satisfy only for large values of $A$

Problem 10: Fish and Belytschko Problem 3.4

Repeat Problem 3.3 with the trial solution $u\left ( x \right ) = \alpha_{0} + \alpha_{1}\left ( x-3 \right ) + \alpha_{2}\left ( x-3 \right )^{2}$.

Find

Consider a trial (candidate) solution of the form $u\left ( x \right ) = \alpha_{0} + \alpha_{1}\left ( x-3 \right ) + \alpha_{2}\left ( x-3 \right )^{2}$ and a weight function of the same form. Obtain a solution to the weak form in Problem 3.1. Check the equilibrium equation in the strong form in Problem 3.1; is it satisfied?

Check the natural boundary condition; is it satisfied?

Given

The weak from is given by $\int_{1}^{3} \frac{\mathrm{d} w }{\mathrm{d} x}AE\frac{\mathrm{d} u}{\mathrm{d} x}dx = -0.1(wA)_{x=1} + \int_{1}^{3} 2xwdx, \forall w$ with $w \left(3 \right) = 0$

Solution

The trial solution $u\left ( x \right )$ and the weight function $w\left ( x \right )$ can be written as follows:

$u\left ( x \right ) = \alpha_{0} + \alpha_{1}\left ( x-3 \right ) + \alpha_{2}\left ( x-3 \right )^{2}$

$w\left ( x \right ) = \beta_{0} + \beta_{1}\left ( x-3 \right ) + \beta_{2}\left ( x-3 \right )^{2}$

When:

$w \left(3 \right) = 0; \Rightarrow \beta_{0} = 0$

and for:

$u \left(3 \right) \Rightarrow \alpha_{0} = 0.001$

$u\left ( x \right )$ and $w\left ( x \right )$ can therefore be re-written as follows:

$u\left ( x \right ) = 0.001 + \alpha_{1}\left ( x-3 \right )+\alpha_{2}\left ( x-3 \right )^{2}$

$w\left( x \right) = \beta_{1} \left(x - 3 \right)+\beta_{2}\left ( x-3 \right )^{2}$

Whith their respective derivatives:

$\frac{du}{dx} = \alpha_{1}+2\alpha_{2}\left(x-3\right)$

$\frac{dw}{dx} = \beta_{1}+2\beta_{2}\left(x-3\right)$

Now substituting into the weak form yields:

$\int_{1}^{3} AE[\beta_{1}+2\beta_{2}\left(x-3\right)][\alpha_{1}+2\alpha_{2}\left(x-3\right)]dx = -0.1[A(\beta_{1}\left(x-3\right)+\beta_{2}\left(x-3\right)^{2})]_{x=1} + \int_{1}^{3} 2x[\beta_{1}\left(x-3\right)+\beta_{2}\left(x-3\right)^{2}]dx$

which can be re-arranged as:

$\int_{1}^{3} AE[\beta_{1}+2\beta_{2}\left(x-3\right)][\alpha_{1}+2\alpha_{2}\left(x-3\right)]dx - 0.1[A(\beta_{1}\left(x-3\right)+\beta_{2}\left(x-3\right)^{2})]_{x=1} - \int_{1}^{3} 2x[\beta_{1}\left(x-3\right)+\beta_{2}\left(x-3\right)^{2}]dx = 0$

Evaluating the integral we get:

$AE \left(2\alpha_{1}\beta_{1}-4\alpha_{1}\beta_{2}-4\alpha_{2}\beta_{1}+\frac{32\alpha_{2}\beta_{2}}{3}\right)- 0.1[A(4\beta_{2}-2\beta_{1})] - \left(8\beta_{2}-\frac {20\beta_{1}}{3}\right)=0$

Factoring out $\left( \beta_{1}\right)$ and $\left( \beta_{2}\right)$ we get:

$\beta_{1}\left[2AE\left(\alpha_{1}-2\alpha_{2}\right)-0.2A+\frac{20}{3}\right]+\beta_{2}\left[4AE\left(\frac{8}{3}\alpha_{2}-\alpha_{1}\right)+0.4A-8\right]=0$

$\displaystyle {{\beta }_{1}}$ and $\displaystyle {{\beta }_{2}}$ must be zero since they were arbitrarely chosen and because the above equation must hold. Therefore we can write both terms of the above equation equal to zeor as folows:

$\left[2AE\left(\alpha_{1}-2\alpha_{2}\right)-0.2A+\frac{20}{3}\right]=0$

and

$\left[4AE\left(\frac{8}{3}\alpha_{2}-\alpha_{1}\right)+0.4A-8\right]=0$

Solving for $\displaystyle {{\alpha }_{1}}$ and $\displaystyle {{\alpha }_{2}}$ we find that:

$\displaystyle {{\alpha }_{1}}=\frac{\frac{116}{3}+0.2A}{2AE}$

and

$\displaystyle {{\alpha }_{2}}=\frac{8}{AE}$

The solution to the weak form can be written as

$u(x)= 0.001+\frac{\frac{116}{3}+0.2A}{2AE}(x-3)-\frac{8}{AE}(x-3)^{2}$

$\displaystyle (Eq. 10.1)$

Checking if equilibrium in the strong form is satisfied, substituted into the equations:

$\frac{d}{dx} \left( AE \frac{du}{dx} \right) + 2x = 0$

Yields the following:

$\displaystyle {2x-16=0}$

which does not satisfy.

Checking if the natural boundary condition is satisfied, the following equation can be used:

$\sigma \left(1 \right) = \left(E \frac{du}{dx} \right)_{x=1} = 0.1$

Problem 11: Solving using WRF

Solving a differential equation using Weight Residual Form. From lecture 17

Given

The partial differential equation is:

$\frac{\partial}{\partial x}\left(2\frac{\partial u}{\partial x}\right) + 3 = 0$

$\displaystyle (Eq.11.1)$

The boundary conditions are:

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

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

As weighting funtion should be used:

$b_{i} =1+a_{i}\cos \left ( ix \right )+b_{i}\sin \left ( ix \right )$

$\displaystyle (Eq.11.2)$

Find

11.1 Find an approximate solution of the form $U^{h}=\sum_{i=0}^{n}d_{i}b_{i}$ for n=2

11.2 Find two equations that enforce the boundary conditions

11.3 Project the weight residues

11.4 Display the equations in matrix form

11.5 Solve for $\displaystyle d$

11.6 Construct $U^{h}$ and plot $U^{h}$ and $u$

11.7 Repeat 1.1 to 1.6 for n = 4 and n = 6

Solution

11.1 Find an approximate solution with n=2

When n = 2 we have:

$d_{i}=\begin{bmatrix} d0 & d1 & d2 & d3 & d4 \end{bmatrix}$

$\displaystyle (Eq.11.3)$

and

$b_{i}=\begin{bmatrix} 1 & \cos\left ( x\right) & \sin\left ( x\right) & \cos\left ( 2x\right) & \sin\left ( 2x\right)\end{bmatrix}$

$\displaystyle (Eq.11.4)$

Therefore:

$U^{h}=d0+d1\cos\left ( x \right )+d2\sin\left ( x\right )+ d3\cos\left ( 2x \right )+d4\sin\left ( 2x\right )$

$\displaystyle (Eq. 11.5)$

11.2 Find two equations that enforce the boundary conditions

The boundary conditions are:

$u\left(1\right)=0=d0+d1\cos\left(1\right)+d2\sin\left(1\right)+d3\cos\left(2\right)+d4\sin\left(2\right)$

$\displaystyle (Eq. 11.6)$

And

$\frac{du}{dx}(0)=-4 = d2+2d4$

$\displaystyle (Eq. 11.7)$

11.3 Project the weight residues

In this case changing the function u(x) by the approximated function $U^{h}$ in the PDE (Eq.11.1), we found that P(u) is:

$P\left ( U^{h} \right )=2\frac{\partial^2 U^{h}}{\partial x^2}+3=2d2+3$

$\displaystyle (Eq.11.8)$

The (Eq.11.8) can be written as:

$P\left ( U^{h} \right )= -2\begin{bmatrix}0 & \cos\left(x\right) & \sin\left(x\right) & 4\cos\left(2x\right) &4\sin\left(2x\right)\end{bmatrix}\begin{bmatrix} d0\\ d1\\ d2\\ d3\\ d4\end{bmatrix}+3$

$\displaystyle (Eq.11.9)$

Projecting the residue we have:

$\int_{a}^{b}b(x)P\left ( U^{h} \right )= 0$

$\displaystyle (Eq. 11.10)$

After substitution of (Eq.11.4) and (Eq.11.9) in (Eq.11.10) the following equation is obtained

$2\int_{0}^{1}\begin{bmatrix} 1\\ \cos\left (x\right ) \\ \sin\left ( x\right )\\ \cos\left (2x\right ) \\ \sin\left ( 2x\right ) \end{bmatrix}\begin{bmatrix} 0 &\cos\left(x\right) & \sin\left(x\right)&\cos\left(2x\right) & \sin\left(2x\right)\end{bmatrix}dx\begin{bmatrix} d0\\ d1\\ d2\\ d3\\ d4\end{bmatrix}=3\int_{0}^{1} \begin{bmatrix} 1\\ \cos\left ( x\right ) \\ \sin\left ( x \right ) \\ \cos\left ( 2x\right ) \\ \sin\left ( 2x \right ) \end{bmatrix}dx$

$\displaystyle (Eq. 11.11)$

11.4 Display the equations in matrix form

Performing the product and then the integration indicated in (Eq.1.11) we have:

$\begin{bmatrix} 0 &1.68 & 0.91 &3.63& 5.66\\ 0 & 1.45& 0.70 & 3.55&4.49\\ 0 &0.70 & 0.54 &0.81 & 3.17\\ 0 & 0.88& 0.20 & 3.24& 1.65\\ 0 & 1.12& 0.79& 1.65& 4.75\end{bmatrix}\begin{bmatrix} d0\\ d1\\ d2\\ d3\\ d4\end{bmatrix}=\begin{bmatrix} 3\\ 2.52\\ 1.37\\ 1.36\\ 2.12\end{bmatrix}$

$\displaystyle (Eq. 11.12)$

The latter equation is clearly of the form $\displaystyle Kd=F$. We can observe that K is not symmetric.

11.5 Solve for $d$

In (Eq.11.12) we have three equations; but, as we need to enforce the boundary conditions we take the only one of them and solve it together with (Eq.11.6) and (Eq.11.7). From these equations we now have the system:

$\begin{bmatrix} 0 &\cos\left(1\right) & \sin\left(1\right) &\cos\left(2\right) & \sin\left(2\right)\\ 0 &0 &1 &0 & 2\\ 0 & -1.68 & -0.91 & -3.63 & -5.66\\ 0 &-1.45 & -0.70 & -3.55 & -4.49\\ 0 & -0.70& -0.54 & -0.81 & -3.17 \end{bmatrix}\begin{bmatrix} d0\\ d1\\ d2\\ d3\\ d4\end{bmatrix}=\begin{bmatrix} 0\\ -4\\ -3\\ -2.52\\ -1.37\end{bmatrix}$

$\displaystyle (Eq.11.13)$

After solving the system in (Eq.11.13), by using the code shown in apedix, we obtain d:

$\begin{bmatrix} d0\\ d1\\ d2\\ d3\\ d4\end{bmatrix}=\begin{bmatrix} 0.67\\ 4.88\\ -4.72\\ -0.80\\ -0.36\end{bmatrix}$

11.6 Construct $U^{h}$ and plot $U^{h}$ and $u$

Replacing d in (Eq.11.5) we obtain the approximated solution:

$U^{h}=0.67+4\cos\left(x\right)-4.72\sin\left(x\right)-0.8\cos\left(2x\right)-0.36\sin\left(2x\right)$

$\displaystyle (Eq.11.14)$

On the other hand, the exact solution of the PDE (Eq.11.1) with the given boundary condition is:

$u\left ( x \right )=-\frac{3}{4}x^{2}-4x+\frac{19}{4}$

$\displaystyle (Eq.11.15)$

11.7 Repeat 11.1 to 11.6) for i = 4 and i = 6

Now we repeat the procedure using i = 4 and i=6 in Eq.11.4

As the procedure is the same shown before, we develop a routine using Matlab to perform the products, the integration and to solve the system of equations. The routine is show in the appendix.

In Fig.11.1 we show the exact and the approximated solution for i = 2,4, and 6. In that figure, it can be observed that all the numerical solutions fits very well the analytical solution.

We found a small error for i = 2 (0.0003); however the error lightly increases for i = 4 and decreases again for i = 6. As all the solutions are close to the analytical, we believe that this error is mainly due to round-off error which is increased with i = 4 and i = 6 due to the higher number of operations that the program has to perform. Fig. 11.2 shows the error variation as a function of i

Fig.11.1 Analytical and numerical solution with i=2, 4, and 6

Fig.11.2 Error at x = 0.5

Appendix

Matlab Code:

clear all;
x = sym('x');
 
b = [1 cos(x) sin(x) cos(2*x) sin(2*x)];
db = -[0 cos(x) sin(x) 4*cos(2*x) 4*sin(2*x)];
k = b'*db;
K = 2*int(k,0,1);
f = -3*int(b',0,1);
 
 
K1 = [1 cos(1) sin(1) cos(2) sin(2);
    0 0 1 0 2; 
    0   -1.6829   -0.9194   -3.6372   -5.6646;
         0   -1.4546   -0.7081   -3.5540   -4.4921;
         0   -0.7081   -0.5454   -0.8145   -3.1777];
f1 = [0 -4 -3 -3*sin(1) -3+3*cos(1)]';
 
 
 
d = inv(K1)*f1;
z = [0:0.1:1];
 
for i = 1:11
 
    f2(i) = d(1)+ d(2)*cos(z(i)) + d(3)*sin(z(i)) + d(4)*cos(2*z(i)) + d(5)*sin(2*z(i)); 
    %calculates Uh
    fnal(i) = -(3/4)*z(i)*z(i)-4*z(i)+19/4; %This is the analytic solution
end
plot(z,fnal,'k')
hold on
plot(z,f2,'ko') %plots both solution
 
 
% xlabel('x')
% ylabel('U(x)')
% legend('Numerical solution','Analytical solution')
% title('Numerical and analytical solution with i = 2')
err(1) = abs(fnal(6)-f2(6));
 
%Now for n = 4 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
x = sym('x');
%w = 2*pi;
w = 1;
b = [1 cos(x) sin(x) cos(2*x) sin(2*x) cos(3*x) sin(3*x) cos(4*x) sin(4*x)];
db = -[0 cos(x) sin(x) 4*cos(2*x) 4*sin(2*x) 9*cos(3*x) 9*sin(3*x) 16*cos(4*x) 16*sin(4*x)];
k = b'*db;
K = 2*int(k,0,1);
f = -3*int(b',0,1);
 
 
K1 = [1 cos(1) sin(1) cos(2) sin(2) cos(3) sin(3) cos(4) sin(4);
    0 0 1 0 2 0 3 0 4; 
    0   -1.6829   -0.9194   -3.6372   -5.6646 -0.8467  -11.9400    6.0544  -13.2291;
         0   -1.4546   -0.7081   -3.5540   -4.4921 -2.3890  -10.0934    2.3159  -12.9056;
         0   -0.7081   -0.5454   -0.8145   -3.1777  2.6520   -5.7946    8.3210   -3.8212;
         0   -0.8885   -0.2036   -3.2432   -1.6536   -5.8472   -5.4267   -6.5293  -11.4354;
         0   -1.1230   -0.7944   -1.6536   -4.7568    2.8479   -9.2993   11.2230   -8.0195;
         0   -0.2654    0.2947   -2.5987    1.2657   -8.5809   -0.0597  -14.9652   -7.9177;
         0   -1.1215   -0.6438   -2.4119   -4.1330   -0.0597   -9.4191    6.7927  -11.9619];
 
f1 = [0 -4 -3 -3*sin(1) -3+3*cos(1)  -3/2*sin(2)  -3/2+3/2*cos(2)  -sin(3)  -1+cos(3)]';
 
 
 
d = inv(K1)*f1;
z = [0:0.1:1];
 
for i = 1:11
 
    f2(i) = d(1)+ d(2)*cos(z(i)) + d(3)*sin(z(i)) + d(4)*cos(2*z(i)) + d(5)*sin(2*z(i))+ ...
        d(6)*cos(3*z(i)) + d(7)*sin(3*z(i))+d(8)*cos(4*z(i)) + d(9)*sin(4*z(i)) ; 
    %calculates Uh
    fnal(i) = -(3/4)*z(i)*z(i)-4*z(i)+19/4; %This is the analytic solution
end
plot(z,f2,'r+') %plots both solution
 
% xlabel('x')
% ylabel('U(x)')
% legend('Analytical solution','Numerical solution i = 2','Numerical solution i = 4')
%title('Numerical and analytical solution with i = 2')
err(2) = abs(fnal(6)-f2(6));
 
%Now for n = 6 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x = sym('x');
 
b = [1 cos(x) sin(x) cos(2*x) sin(2*x) cos(3*x) sin(3*x) cos(4*x) sin(4*x) cos(5*x) sin(5*x) cos(6*x) sin(6*x)];
db = -[0 cos(x) sin(x) 4*cos(2*x) 4*sin(2*x) 9*cos(3*x) 9*sin(3*x) 16*cos(4*x) 16*sin(4*x) 25*cos(5*x) 25*sin(5*x) 36*cos(6*x) 36*sin(6*x)];
k = b'*db;
K = 2*int(k,0,1)
f = -3*int(b',0,1)
 
 
K1 = [1 cos(1) sin(1) cos(2) sin(2) cos(3) sin(3) cos(4) sin(4) cos(5) sin(5) cos(6) sin(6);
    0 0 1 0 2 0 3 0 4 0 5 0 6; 
    0   -1.6829   -0.9194   -3.6372   -5.6646 -0.8467  -11.9400    6.0544  -13.2291 9.5892   -7.1634    3.3530   -0.4780;
         0   -1.4546   -0.7081   -3.5540   -4.4921 -2.3890  -10.0934    2.3159  -12.9056  5.8942  -10.5012    3.5255   -6.4233;
         0   -0.7081   -0.5454   -0.8145   -3.1777  2.6520   -5.7946    8.3210   -3.8212  10.1693    3.5658    3.8920   10.2830;
         0   -0.8885   -0.2036   -3.2432   -1.6536   -5.8472   -5.4267   -6.5293  -11.4354 -3.5224  -17.4622    2.3591  -20.0375;
         0   -1.1230   -0.7944   -1.6536   -4.7568    2.8479   -9.2993   11.2230   -8.0195  15.7044    1.1704    9.7280   11.2633;
         0   -0.2654    0.2947   -2.5987    1.2657   -8.5809   -0.0597  -14.9652   -7.9177 -14.4580  -21.2815   -3.3419  -31.5244;
         0   -1.1215   -0.6438   -2.4119   -4.1330   -0.0597   -9.4191    6.7927  -11.9619  14.1221   -8.2745   16.2354   -0.0450;
         0    0.1447    0.5201   -1.6323    2.8057   -8.4179    3.8209  -17.9787   -2.2910  -22.1815  -16.8011  -14.4089  -32.1113;
         0   -0.8066   -0.2388   -2.8588   -2.0049   -4.4537   -6.7285   -2.2910  -14.0213    6.1837  -19.8920   18.8700  -18.3258;
         0    0.2358    0.4068   -0.5636    2.5127   -5.2049    5.0840  -14.1962    3.9576  -23.6399   -4.5977  -27.0203  -19.8074;
         0   -0.4200    0.1426   -2.7940    0.1873   -7.6613   -2.9788  -10.7527  -12.7309   -4.5977  -26.3601   13.2909  -33.5657];
 
 
f1 = [0 -4 -3 -3*sin(1) -3+3*cos(1)  -3/2*sin(2)  -3/2+3/2*cos(2)  -sin(3)  -1+cos(3)   -3/4*sin(4) -3/4+3/4*cos(4) -3/5*sin(5) -3/5+3/5*cos(5)]';
 
d = inv(K1)*f1;
z = [0:0.1:1];
 
for i = 1:11
 
    f2(i) = d(1)+ d(2)*cos(z(i)) + d(3)*sin(z(i)) + d(4)*cos(2*z(i)) + d(5)*sin(2*z(i))+ ...
        d(6)*cos(3*z(i)) + d(7)*sin(3*z(i))+d(8)*cos(4*z(i)) + d(9)*sin(4*z(i)) + d(10)*cos(5*z(i))+...
        d(11)*sin(5*z(i))+ d(12)*cos(6*z(i)) + d(13)*sin(6*z(i));
 
    %calculates Uh
    fnal(i) = -(3/4)*z(i)*z(i)-4*z(i)+19/4; %This is the analytic solution
end
plot(z,f2,'bd') %plots both solution
 
xlabel('x')
ylabel('U(x)')
legend('Analytical solution','Numerical solution i = 2','Numerical solution i = 4','Numerical solution i = 6')
%title('Numerical and analytical solution with i = 2')
err(3) = abs(fnal(6)-f2(6));
 
h = [2 4 6 ];
figure, plot(h,err,'k')
xlabel('Number of funtions i')
ylabel('error at x = 0.5')

Contributing Members

James Roark

Abhijeet Bhalerao

Fernando Casanova

Yu Hao

Harrison Sheffield

Cedric Adam