please use MATLAB to solve this problem, thank you very much! Part 1. 2pBVP solv
ID: 3701370 • Letter: P
Question
please use MATLAB to solve this problem, thank you very much!
Part 1. 2pBVP solvers Theory and problem statement. Consider a simple 2pBVP Introduce an equidistant grid on 10,L) with ?-L/(N + 1). The derivative is discretized with a symmetric finite difference quotient, so that we obtain This is a system of N equations for the N unknowns y, yN, which approximate the exact solution values, y(i), y2(rN). If we denote the systein by F(y) = 0, then the equations are (note how the boundary values enter -22 f) F,(y) = 12 PrtScn Home F3 F5 F6 4Explanation / Answer
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 side function f of the 2pBVP at the N 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 beginning 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 tridiagonal matrices of the structure you need. In Python, you have scipy.linalg.toeplitz for this. Optionally, you can use the command diag (numpy.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. For Python, check the scipy.sparse module with scipy.sparse.diags. 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 accordance with its material properties and the applied load. According to elasticity theory, the deflection u is governed by the differential equations
M00 = q(x) u 00 = 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.
2.The Sturm–Liouville problem does just that, and and has applications in e.g. materials science (oscillations in beams; buckling modes in structural analysis) as well as in microphysics (quantum mechanics). Sturm–Liouville problems are eigenvalue problems for differential operators. The standard formulation is
d /dx (p(x) dy /dx) ? q(x)y = ?y
y(a) = 0; y(b) = 0
where p(x) > 0 and q(x) ? 0. The problem is to find eigenvalues ? and eigenfunctions y(x) satisfying this 2pBVP. Eigenfunctions are often called “eigenmodes” or just “modes,” because they can be interpreted as oscillation modes of various frequencies. In quantum mechanics the eigenmodes are called “wave functions,” and the eigenvalues represent “energy levels” of the corresponding “states.” Note that the boundary conditions in a Sturm–Liouville problem are always homogeneous; apart from the Dirichlet condition above, one also encounters homogeneous Neumann conditions, say y 0 (a) = 0. Discretization. The Sturm–Liouville problem is discretized by standard symmetric 2nd order finite differences to obtain a matrix eigenvalue problem T?xy = ??xy where T?x is a symmetric, tridiagonal N × N matrix, which depends on the functions p, q and ?x. The discretization has exactly N real eigenvalues ??x = ? + O(?x 2 ), due to symmetry. This reflects the fact that the differential operator is self–adjoint and therefore possesses an infinite sequence of real eigenvalues and orthogonal eigenfunctions. You generate the matrix using the diag command as described before. The eigenvalues of a matrix A are computed in Matlab using the command eig(A). Depending on how you call this function, you can compute either eigenvalues only, or eigenvalues and the corresponding eigenvectors. The latter are discrete versions of the eigenfunctions, and can be plotted to visualize the eigenfunctions. Use help in Matlab to find out more
Related Questions
Navigate
Integrity-first tutoring: explanations and feedback only — we do not complete graded work. Learn more.