Attached project_4 pdf
% twopBVP code function y = twopBVP(fvec,alpha,beta,L,N) h = L/(N+1); b = fvec*h^2; b(1) = b(1) - alpha; b(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 corresponding 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); beta = 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)') error(i) = norm(yvec-y); end figure(i+1) %figure(1) loglog(Ng,error,'r-o') title('Approximation Error graph') xlabel('n') ylabel('error') % % 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 corresponding % 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 boundary 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 . . . . . . 0 1 −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 solver 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 respect, 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 corresponding f . Don’t forget that boundary conditions must match. For your chosen problem, verify that your method by plotting the error in a log-log diagram in the usual way to demonstrate second order 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 solver repeatedly 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 vector 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 correct 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, but higher order methods also have sparse representations that should be 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 beam’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 reaches all the way out to the boundary. Please be careful with your SI units – steel is a strong material, so loads are on the