% twopBVP code
function y = twopBVP(fvec,alpha,beta,L,N)
h = L/(N+1);
= fvec*h^2;
(1) = b(1) - alpha;
(N) = b(N) - beta;
e = ones(N,1);
A = spdiags([e,-2*e,e],[-1,0,1], N,N);
y = [alpha;A\b;beta];
end
%
% Project 4 - Part 1 2pBVP solvers
%
% Simple linear 2pBVP
% y'' = f(x), y(0) = alpha, y(L) = beta
%
%
% clear all variables and plot windows
clear all; close all;
% Test the solver
%
% select the function ya and differentiate twice to get the co
esponding fa.
%
%
%k=2*pi;
%% ya = inline('sin(x)/x','x');
ya = inline('sin(7*2*pi*x).*exp(x)', 'x');
%ya'=k*cos*exp + sin*exp
%ya''=(-k^2+1)*sin*exp+2*k*cos*exp
fa = inline('(-(7*2*pi)^2+1)*sin((7*2*pi)*x).*exp(x)+2*(7*2*pi)*cos(7*2*pi*x).*exp(x)','x');
L = 1;
% set boundary conditions
%alpha = 1
alpha = ya(0);
eta = ya(L);
% set up the grid size
Ng = 2.^(8:12);
for i = 1:length(Ng)
N = Ng(i)
h = L/(N+1);
% grid points
x = (0:h:L)';
yvec = ya(x); % the known analytic solution
fvec = fa(x(2:N+1)); % interior points
y = twopBVP(fvec,alpha,beta,L,N);
figure(i)
plot(x,y,'r-o',x,yvec,'b-x')
legend('approximation','actual')
xlabel('x')
ylabel('y(x)')
e
or(i) = norm(yvec-y);
end
figure(i+1)
%figure(1)
loglog(Ng,e
or,'r-o')
title('Approximation E
or graph')
xlabel('n')
ylabel('e
or')
%
% Project 4 - Part 1.2 Beam Equation
%
% Simple 2pBVP
% y'' = f(x,y), y(0) = alpha, y(L) = beta
%
%
% clear all variables and plot windows
clear all; close all;
L = 10; % unit (m)
E = 1.9E11; % unit (N/m^2)
% grid
N = 10^3-1;
h = L/(N+1);
x = (0:h:L)';
% unit (N/m)
qvec = -50*10^3*ones(N,1);
% unit (Nm)
M = twopBVP(qvec,0,0,L,N);
%
Ivec = 10^(-3)*(3-2*(cos(pi*x(2:N+1)/L)).^12);
fvec = M(2:N+1)./(E*Ivec);
u = twopBVP(fvec,0,0,L,N);
figure(1)
plot(x,M,'r-o')
xlabel('x')
ylabel('M')
figure(2)
plot(x,u,'b-o')
xlabel('x')
ylabel('u')
disp('Deflection at the midpoint in millimeter (mm)')
u((N+1)/2+1)*1E3
%
% Project 4 - Part 1 2pBVP solvers
%
L = 1;
% set up the grid size
Ng = 2.^(8:12);
for i = 1:length(Ng)
N = Ng(i)
h = L/(N+1);
h = L/(N+1);
e = ones(N,1);
A = 1/(h^2)*spdiags([e,-2*e,e],[-1,0,1], N,N);
% [V,D] = EIGS(A,3,'SM') returns a diagonal matrix D of A's 3 smallest magnitude
% eigenvalues and a matrix V whose columns are the co
esponding
% eigenvectors.
[V,D] = eigs(A,3,'SM');
Project 4
This assignment can be worked on and submitted in groups of up to three
persons.
Goals. In this assignment, the goal is to construct solvers for two-point
oundary value problems as well as Sturm–Liouville eigenvalue problems.
The objective is to learn the basic finite difference methodology by imple-
menting it in detail. As it is too complicated to make the solver adaptive (i.e.,
the mesh width varies along the solution), we only construct solvers for
equidistant grids. Likewise, because nonlinear problems present additional
difficulties, we will only consider linear problems.
Part 1. 2pBVP solvers
Theory and problem statement. Consider a simple 2pBVP
y′′ = f(x, y)
y(0) = α; y(L) = β
Introduce an equidistant grid on [0, L] with ∆x = L/(N + 1). The derivative
is discretized with a symmetric finite difference quotient, so that we obtain
yi+1 − 2yi + yi−1
∆x2
= f(xi, yi)
y0 = α; yN+1 = β.
This is a system of N equations for the N unknowns y1, y2, . . . , yN , which
approximate the exact solution values, y(x1), y(x2), . . . , y(xN). If we denote
the system by F (y) = 0, then the equations are (note how the boundary
values enter)
F1(y) =
α− 2y1 + y2
∆x2
− f(x1, y1)
Fi(y) =
yi−1 − 2yi + yi+1
∆x2
− f(xi, yi); i = 1, . . . N
FN(y) =
yN−1 − 2yN + β
∆x2
− f(xN , yN).
1
In the nonlinear case, one would have to solve F (y) = 0 iteratively, using
Newton’s method. However, trying to keep things simple, we will only study
linear problems here, so as to avoid iterative equation solving.
In fact, we are only going to consider the very simplest linear problems,
where f(x, y) ≡ f(x). This means that the right-hand side of the BVP is
independent of y. In this case, we see that our system F (y) = 0 reduces to
a linear system of equations,
1
∆x2
−2 1 0 . . .
1 −2 1
1 −2 1
. . .
XXXXXXXXXX −2
y1
y2
...
yN
=
−α/∆x2 + f(x1)
f(x2)
...
f(xn−1)
−β/∆x2 + f(xN)
This is a tridiagonal linear system. The matrix is sparse. Solving the system
has low complexity. The LU decomposition runs in O(N) time (3N opera-
tions, to be precise; compare to the N3/3 operations needed for a dense, full
matrix). The forward and back substitutions also run in O(N) time (a total
of 2N operations). Therefore, the solution effort is moderate even when N
is large, and is proportional to the number of grid points.
Task 1.1 (1p) Write a Matlab 2pBVP solve
function y = twopBVP(fvec, alpha, beta, L, N)
that solves y′′ = f(x) with boundary conditions y(0) = α and y(L) = β on an
equidistant grid having N interior points. Test your solver by implementing
any right-hand side function of your choice, so that you know what the exact
solution is. It is necessary to verify that the code functions properly in every
espect, even when changing the boundary conditions.
Note that you cannot take a “too simple” right hand side. For example,
if you take f(x) ≡ 1, the solution is a 2nd degree polynomial, and a 2nd
order method therefore solves the problem exactly. To choose a more diffi-
cult problem, select the function y itself, and differentiate twice to find the
co
esponding f . Don’t forget that boundary conditions must match. Fo
your chosen problem, verify that your method by plotting the e
o
in a log-log diagram in the usual way to demonstrate second orde
convergence. Give all details about your computational setup.
Here are some useful hints for constructing your solver:
1. To simplify your work with Task 1.2 let your solver take fvec as an
input vector. That is, before calling twopBVP, evaluate the right hand
2
side function f of the 2pBVP at theN interior grid points x1, x2, . . . , xN
and store the results in fvec.
2. Do not plot the result from within twopBVP as you will use this solve
epeatedly below. Instead, plot the result by giving a plot command in
the script (main program) you use to call twopBVP.
3. Note that when you plot the solution, you have obtained y only on
interior points. Make sure to append the boundary values at the begin-
ning and end of your solution vector, so that you can plot the solution
all the way from boundary to boundary.
4. In Matlab you will find the function toeplitz, which generates tridi-
agonal matrices of the structure you need. Optionally, you can use the
command diag, in the following manner. If you have an N − 1 vecto
sup, an N−1 vector sub and an N vector main, then you can construct
a tridiagonal N ×N matrix with subdiagonal sub, main diagonal main
and superdiagonal sup by the command
A = diag(sub,-1) + diag(main,0) + diag(sup,1);
Using this technique, all you need to do is to first generate the three
vectors sub, main, sup, and install them in their right positions in
the matrix. After constructing the matrix, using either technique, you
co
ect the elements in the first and last row of A, if necessary, so
that the boundary conditions are properly represented. The commands
toeplitz and diag are very convenient as they allow you to avoid using
for-loops and instead work with a vectorized code. However these
commands have the drawback that you work with full matrices without
exploiting the sparse tridiagonal structure. For higher performance,
you can check Matlab’s help for the commands sparse and spdiags.
Thus, using the sparsity of the tridiagonal matrices, you can choose a
matrix representation that makes your solver much faster, and also
enables you to use very large values of N . (A better alternative in high
precision computations is of course to work with higher order methods,
ut higher order methods also have sparse representations that should
e exploited.)
The Beam Equation. An elastic beam under load is deflected in accor-
dance with its material properties and the applied load. According to elas-
3
ticity theory, the deflection u is governed by the differential equations
M ′′ = q(x)
u′′ = M(x)/(EI),
where q(x) is the load density (N/m); M(x) is the bending moment (Nm);
E is Young’s modulus of elasticity (N/m2); I is the beam’s cross-section
moment of inertia (m4); and u is the beam’s centerline deflection (m).
The beam is supported at its ends, at x = 0 and x = L. This means that
there is no deflection there: u(0) = u(L) = 0. Further, assuming that the
eam’s ends do not sustain any bending moment, we also have the boundary
conditions M(0) = M(L) = 0. We have then obtained a 2pBVP, of order 4.
The task is to solve for the deflection u. Thus, given the load vector q, you
first solve a 2pBVP for the bending moment M ; then, using your calculated
M vector as data, you solve another 2pBVP for the deflection u.
Task 1.2 (2p) Use your 2pBVP solver to write a script for solving the beam
equation for a heavy-duty railway flatcar. Its support beams have length
L = 10 m and elasticity module E = 1.9 · 1011 N/m2 (construction steel).
The beam’s web is “taller” at the center section. This is modeled by a cross-
section moment of inertia that varies along the beam, according to
I(x) = 10−3 ·
(
3− 2 cos12 πx
L
)
.
Solve the problem for this beam shape, with a load of q(x) = −50 kN/m
(i.e., a load of five metric tonnes per meter), and plot the computed de-
flection. Don’t forget to insert the boundary conditions so that your plot
eaches all the way out to the boundary. Please be careful with your SI units
– steel is a strong material, so loads are on the