Problem 1: Numerical Differentiation (50 points) In this problem, we will use numerical differentiation and Euler-Bernoulli beam theory to approximate the elastic curve of a clean airplane wing under...

Problem 1: Numerical Differentiation (50 points)
In this problem, we will use numerical differentiation and Euler-Bernoulli beam theory to approximate
the elastic curve of a clean airplane wing under elliptical lift distribution. The wing has length L and,
for simplicity, is assumed to be a uniform cantilevered beam as shown in Figure 1. Recall that the
governing ODE for the deflection v of a homogeneous beam with constant cross-section is given by
EIv(4)(x) = w(x) , (1)
where
• E is the beam’s Young’s modulus,
• I is the area moment of inertia of the cross section, and
• w is the distributed transverse load with positive direction as shown in the figure.
Since the wing is fixed at the left end and free at the right, the deflection must satisfy the following
oundary conditions:
• Zero displacement at x = 0 =) v(0) = 0 , (2a)
• Zero rotation at x = 0 =) v0(0) = 0 , (2b)
• Zero bending moment at x = L =) v00(L) = 0 , (2c)
• Zero shear force at x = L =) v000(L) = 0 . (2d)
1. (40 points) Problem Set-up and Implementation
(a) (6 points) Assuming that v has as many derivatives as needed, use the data v(x � 2h) ,
v(x � h) , v(x) , v(x + h) , and v(x + 2h) to show that the 4th derivative of v at the point x
can be approximated as
v(4)(x) ⇡
v(x� 2h)� 4v(x� h) + 6v(x)� 4v(x+ h) + v(x+ 2h)
h4
. (3)
(b) (1 point) What is the order of the approximation?
We set up a uniform grid along the beam, i.e. on [0 , L], by letting h =
L
n
and
xj = jh , j = 0 , 1 , . . . , n ,
where n+ 1 is the total number of grid points (including the endpoints).
Furthermore, we let vj be the approximation of the deflection at the mesh point xj , i.e., vj ⇡ v(xj) ,
and let wj = w(xj)
h4
EI
.
2
(c) (2 points) Use the approximation of Equation (3) to evaluate Equation (1) at xj, thus relating
wj and the approximations vj�2 , vj�1 , vj , vj+1 , vj+2 of the beam deflection at the grid points
xj�2 , xj�1 , xj , xj+1 , xj+2 respectively.
Note that the above result is problematic near the endpoints as it requires knowledge of data beyond
the beam. For example, using the result of part (c) with j = 1 would require knowledge of v�1
which does not exist within the grid domain. We consider this below.
(d) (1.5 points) For j = 1 , use the result of part (c) to obtain an equation relating v1 , v2 , v3 , v�1 ,
and w1 .
Hint: Remember that v(0) = 0 from the boundary condition of Equation (2a).
(e) (1 point) Recall the central difference formula for the 1st derivative:
v0(x) ⇡
v(x+ h)� v(x� h)
2h
. (4)
State the order of the central difference scheme of Equation (4).
(f) (2 points) Use Equation (4) and the boundary condition of Equation (2b) to express v�1 in
terms of v1 .
(g) (2 points) Use the results of parts (d) and (f) to obtain an equation relating v1 , v2 , v3 , and
w1 .
(h) (1.5 points) Next, for j = n � 1 , use the result of part (c) to obtain an equation relating
vn�3 , vn�2 , vn�1 , vn , vn+1, and wn�1 .
(i) (1 point) Likewise, recall the central difference formula for the 2nd derivative:
v00(x) ⇡
v(x� h)� 2 v(x) + v(x+ h)
h2
. (5)
State the order of this central difference scheme.
(j) (2 points) Use Equation (5) and the boundary condition of Equation (2c) to express vn+1 in
terms of vn�1 and vn .
(k) (2 points) Use the results of parts (h) and (j) to obtain an equation relating vn�3 , vn�2 , vn�1 , vn ,
and wn�1 .
(l) (1.5 points) Finally, for j = n , use the result of part (c) to obtain an equation relating
vn�2 , vn�1 , vn , vn+1 , vn+2, and wn .
(m) (2 points) It can be shown that a second order approximation for the 3rd derivative is
v000(x) ⇡
�v(x� 2h) + 2v(x� h)� 2v(x+ h) + v(x+ 2h)
2h3
. (6)
Use Equation (6) and the boundary condition of Equation (2d) to express vn+2 in terms of
vn�2 , vn�1 and vn+1 .
3
(n) (2.5 points) Insert your results of parts (j) and (m) into that of part (l) to obtain an equation
elating vn�2 , vn�1 , vn , and wn .
(o) (6 points) Combine the results of parts (c), (g), (k), and (n) to obtain the system of equations
to be solved for the unknown deflections v1 , v2 , . . . , vn .
Specifically, letting v = (v1 , v2 , . . . , vn)T and w = (w1 , w2 , . . . , wn)T , determine the matrix
A 2 Rn⇥n such that Av = w .
In what follows, use values for a beam made out of steel:
E = 2.10⇥1011 N/m2, I = 9.1⇥10�6 m4, L = 5m, w(x) = w0
1�
x2
L2
with w0 = 1284 N/m .
(p) (6 points) Use Matlab® to solve for and plot the approximate deflection of the beam fo
n = 24 . Call your script Problem1.m. Hint: Remember to include the point x0 = 0 .
On the same figure, also plot the exact solution given by
v(x) =
w0
3EI
(
1
16
"
1
15
(16L4 + 83L2x2 + 6x4)
1�
x2
L2
+x(3L3 + 4Lx2) arcsin
⇣x
L
⌘i
�
⇡L
8
x3 �
L4
15
�
2. (10 points) E
or Analysis
For a given n , the infinity norm of the approximation e
or is computed as
en = max
0jn
��vj � v(xj)
�� = max
���v0 � v(x0)
��,
��v1 � v(x1)
��, . . . ,
��vn � v(xn)
�� , (7)
where, recall, v(xj) is the exact solution at xj and vj is the approximate deflection at xj .
If p is the order of the approximation, then it is expected that the infinity norm of the e
or would
obey
en = O(h
p) = O(n�p)
as n approaches infinity. In other words, as n ! 1 , we can write
en ⇡ C · n
�p for some constant C . (8)
(a) (4 points) Using Equation (7), add to your script from Part 1(p) to tabulate the infinity norm
of the approximation e
or for the following values of n : 24 , 25 , 26 , 27 , 28 , 29 , 210.
(b) (3 points) Plot the result of part (a) on a log-log plot and use Matlab® ’s polyfit function
to determine the slope of the best-fit line through the points (log n , log en) .
(c) (3 points) Noticing from Equation (8) that log en ⇡ �p log n + logC , deduce the order p of
the approximation. Is this what you expected in light of the order(s) of the approximations
given in Equations (3), (4), (5). Why or why not?
4
L
latexit sha1_base64="XPq1QFTrF6Vxuh0QgS7iDa0wsXU=">AAACbHicbVHNSsNAEN7G
XqjcRFoviqSQq6LHoxYMHBVOFNpTNdmKX
Jhd9JSQp/Aqz6cL+EzuKk9GOvAsh/f/H0zE6ZSGHTdz4qztLyyura+Ud3c2t7ZrdX32kZlmoPPlVT6JWQGpEjAR4ESXlINLA4lPIfD28L/PAJthEqecJJCELPXRESCM7TU432v1nCb7szoIvDmoEHm9tCrV/xuX/EshgS5ZMZ0PDfFIGcaBZcwrXYzAynjQ/YKHQsTFoMJ8pnSKT2xTJ9GStuXIJ2xvzNyFhsziUMbGTMcmL++gvzP18kwug5ykaQZQsJ/GkWZpKhoMTbtCw0c5cQCxrWwWikfMM042uWUuoRxaYYczMhm2orFZ4WPBQ5s2ZQyrdW4FFrIQqWkKVcIQyX75cXMBuA6yCGzn0hxag/h/V37ImifN72Lpvt42WjdzE+yTg7JMTkjHrkiLXJHHohPOAHyRt7JR+XLOXAOnaOfUKcyz9knJXNOvwEa7L76
latexit
v
latexit sha1_base64="WGIK5WkfpK+jRuSo4CjpAT0kjc8=">AAACbHicbVHLTgIxFC3jC/EF6s6YNBKNKzKjJrokunEJiTwSmJhOuUBjZzpp70DIhC9wqx/nT/gNdoCFI96k6cm5r3PvDWIpDLruV8HZ2Nza3inulvb2Dw6PypXjtlGJ5tDiSirdDZgBKSJooUAJ3VgDCwMJneDtKfN3JqCNUNELzmLwQzaKxFBwhpZqTl7LVbfmLoyuA28FqmRljddKodUfKJ6EECGXzJie58bop0yj4BLmpX5iIGb8jY2gZ2HEQjB+ulA6p5eWGdCh0vZFSBfs74yUhcbMwsBGhgzH5q8vI
z9RIcPvipiOIEIeLLRsNEUlQ0G5sOhAaOcmYB41pYrZSPmWYc7XJyXYIwN0MKZmIzbcXss8KnAse2bEyZ1mqaC81koVLS5CsEgZKD/GIWA3Dtp5DYT8Q4t4fw/q59HbRvat5tzW3eVeuPq5MUyRm5INfEI/ekTp5Jg7QIJ0DeyQf5LHw7p86Zc74MdQqrnBOSM+fqB3A8vyQ=
latexit
x
latexit sha1_base64="SfSE2EoYauMFAB2xhbKQ9duPc48=">AAACbHicbVHNSsNAEN7E
XqjcRFoviqSQq6FH04lHB2EIbZLOd2qW
NidWEvoE3jVh/MlfAY3bQ+mdWDZj2/+vpmJUikMet634y4tr6yurW9UNre2d3artb1nozLNIeBKKt2KmAEpEghQoIRWqoHFkYRmNLg
M030Eao5AlHKYQxe01ET3CGlnp8f6nWvYY3MboI/Bmok5k9vNScoNNVPIshQS6ZMW3fSzHMmUbBJYwrncxAyviAvULbwoTFYMJ8onRMTyzTpT2l7UuQTti/GTmLjRnFkY2MGfbNvK8g
O1M+xdh7lI0gwh4dNGvUxSVLQYm3aFBo5yZAHjWlitlPeZZhztckpdorg0Qw7mzWbaisVnhQ8F9m3ZlDKt1bAUWshCpaQpV4giJbvlxUwG4DrMIbOfSHFsD+HPr30RPJ83/IuG93hZv7mdnWSdHJJjckZ8ckVuyD15IAHhBMgH+SRfzo974B66R9NQ15nl7JOSuae/dEy/Jg==
latexit
latexit sha1_base64="W+wmdJ0bQhJgnDf199TUAue2MTE=">AAAB7HicbVBNSwMxEJ2tX7V+VT16CRahXsquFOqx6MVjBWsL7VKyabYNTbJLklXL0
gxYOCePUHefPfmG33oK0PBh7vzTAzL4g508Z1v53C2vrG5lZxu7Szu7d/UD48utdRoghtk4hHqhtgTTmTtG2Y4bQbK4pFwGknmFxnfueBKs0ieWemMfUFHkkWMoJNJj1Wn84H5Ypbc+dAq8TLSQVytAbl
4wIomg0hCOte55bmz8FCvDCKezUj/RNMZkgke0Z6nEgmo/nd86Q2dWGaIwUrakQXP190SKhdZTEdhOgc1YL3uZ+J/XS0x46adMxomhkiwWhQlHJkLZ42jIFCWGTy3BRDF7KyJjrDAxNp6SDcFbfnmVdC5qXr3mebf1SvMqz6MIJ3AKVfCgAU24gRa0gcAYnuEV3hzhvDjvzseiteDkM8fwB87nDx1Bjhs=
latexit
w(x)
Figure 1: Simplified (clean) airplane wing subjected to an elliptical lift distribution.
Problem 2: Partial Differential Equations (50 points)
As we saw in class, the 1-D heat equation can be used to model the diffusion of a property of a
substance (e.g. incompressible fluid) in a one-dimensional medium (rod, bar. etc.) If in addition we
consider the transport of the substance in the medium, we obtain the 1-D advection-diffusion equation
given by
ut + c ux = Duxx , (9)
where u(x, t) describes the state of the fluid (e.g.: temperature, density, energy, etc.), D > 0 is the
diffusion coefficient (assumed to be constant), and c > 0 is the advection velocity (which we will also
assume to be constant).
Furthermore, the initial state of the fluid is given by
u(x, 0) = f(x) for all x . (10)
Let wi,j be the numerical approximation to the true solution u(xi, tj) at the point (xi, tj) on a grid of
(spatial) length L , where, taking h and k as the spatial and temporal step sizes respectively, we have
• xi = ih , for i = 0, 1, 2, 3, . . . ,m+ 1 , with x0 = 0 and xm+1 = (m+ 1)h = L co
esponding to
the left and right boundaries of the grid respectively, and
• tj = jk , for j = 0, 1, 2, 3, . . .
Because Equation (9) has spatial derivatives, we must also specify one boundary conditions for the
variable x . Here, we will do so through the use of periodic boundary conditions. Specifically, we will
assume that the grid repeats indefinitely so that we may write
u(0, t) = u(L, t) for t � XXXXXXXXXX)
Physically, this means that whatever flows out of the boundary x = L must flow back into the grid at
the boundary x = 0 .
5
XXXXXXXXXXpoints) Forward in Time – Central in Space (FTCS) Scheme
In this part, we use the forward Euler in time and central difference in space for the 1st and 2nd
derivatives to discretize Equation (9) over the grid.
(a) (3 points) Using the FTCS finite difference scheme, show that wi,j satisfies
wi,j+1 = wi,j �
�
2
(wi+1,j � wi�1,j) + � (wi�1,j � 2wi,j + wi+1,j) ,
or, equivalently,
wi,j+1 =
✓
� +
�
2
◆
wi�1,j + (1� 2�)wi,j +
✓
� �
�
2
◆
wi+1,j , (12)
where � ⌘
c k
h
and � ⌘
Dk
h2
.
(b) (1.5 points) Draw the stencil for the FTCS scheme of Equation (12). Be sure to label it as
done in class. Is the scheme implicit or explicit? Why?
(c) (3 points) By considering i = 0, 1, 2, . . . ,m � 1,m,m + 1 , show that Equation (12) can be
epresented in matrix form as