2 (part 1) CREDIT WILL BE GIVEN FOR AT MOST ONE BONUS PROBLEM (1) Using the same procedure in class (Jan. 24th lecture), construct the 9 9 system of equations for the coefficients of the cubic spline interpolation through the points (0, 7), (1, 10), (3, 3), and (6, 0). Write the system in the form Ax = b and give A, x, and b. As in class, impose the \"natural\" conditions that the double derivative of the interpolating function must vanish at the end points. Plot the resulting curve, along with the four data points. (2) MATLAB does not impose the two \"natural\" conditions of the double derivative vanishing at the end points. Instead, the two conditions it imposes are that the third derivate of the interpolating function must be continuous at the second point (x = 1) and the second to last point (x = 3). Modify your matrix from question (1) to replace the original \"natural\" conditions with these two new conditions, and give the resulting matrix (only two rows in the matrix need be modified). Plot the curve on the same figure as that from question (1), along with the four data points. (3) Define the function G(x; ) as the solution of the boundary value problem (x)2 1 d2 G e 22 , = 2 dx 2 2 0 < x < 1; G(0; ) = 0 ; G(1; ) = 0 . (1) Here, G(x, ) describes the deflection of a string held on both ends at height equal to 0 due to a downward force that is normally (i.e., according to a Gaussian distribution) distributed around the point with variance 2 . See figure below. 0 0.2 0.4 0.6 0.8 1 Figure 1: String of length one (black) held fixed at both ends under load from an external downward force (blue). The length of the arrow is proportional to the magnitude of the force. What is the resulting deflection of the string? MATLAB can compute integrals using the trapezoidal rule. For example, make a vector x = linspace(0,2,1000), R2 then define a function, e.g., y = x. 2. Then the value of 0 x2 dx is approximated by trapz(x,y). For three values of (0.1), use trapz() to verify that the integral of the function on the right-hand side of Eqn. (1) over 0 x 1 is approximately 1 when = 0.01. (a) Use finite differences to solve numerically for G(x, ) when = 0.7 and = 1, 0.1, and 0.01. Plot all solutions in one figure. As becomes small, does the solution fit your intuition? (b) Use finite differences to solve numerically for G(x, ) when = 0.01 and = 0.2, 0.5 and 0.7. Plot all solutions in one figure. 1 (c) Explain why solutions to parts (3a) and (3b) exist and are unique. You may use MATLAB to compute the determinant of a particular matrix. (d) For part (3b), find the difference in the slope of the string for x > versus x < . That is, numerically compute Gx ( + 0.1; ) Gx ( 0.1; ). This value should be independent of . (e) Compute the deflection of the same string when the loading is given by f (x) = 1 x2 . That is, replace the function on the right-hand side of Eqn. (1) by f (x). Plot the solution. (f) BONUS: Use trapz() to numerically compute the quantity Z 1 u() = G(x; )(1 x2 ) dx . 0 Plot u() over the interval (0, 1) and compare to the solution in part (3e). Submit a printout of your code. (4) Use finite differences to formulate the boundary value problem d2 u = b(x) , dx2 u0 (0) = 0 ; 0 < x < 1; u0 (1) = 0 , (2) as a linear system of equations written in the form Au = b , (3) and answer the following questions. In Eq. (2), u0 (a) denotes the derivative of u at x = a. Assume the following: the interval is divided up into N + 1 evenly distributed points (and therefore N equally sized intervals) with u0 = u(0) and uN = u(1). (a) Give the matrix A and the vectors u and b. Be sure to indicate the dimensions of each. (b) Use MATLAB to compute det A. (c) Use MATLAB to determine the dimension of the nullspace of A. Here, I would suggest using a few small values of N (e.g., 5, 10, 15). (d) Will a solution to Eqn. (3) exist for very possible vector b? Explain. (e) BONUS: What condition must b satisfy in order for there to exist a solution to Eqn. (3)? What is the corresponding condition that b(x) must satisfy for there to exist a solution to Eqn. (2)? (f) BONUS: When b satisfies the condition of part (4e), the solution to Eqn. (3) (and by extension, (2)), will not be unique. In this case, an additional constraint can be imposed to uniquely specify the solution. Often, this is an integral constraint, e.g., Z 1 u(x) dx = 0 . (4) 0 For Eqn. (2) with the constraint in Eqn. (4), formulate a linear system of equations of the form A0 u = b0 . Give the matrix A0 and its dimensions. Use this to find the unique numerical solution for u when b(x) = 1 2 2 e (x)2 2 2 ; = 0.3 , = 0.01 . (5) The command [V,D] = eigs(A,m,'lr') returns the m eigenvalues of A with largest real part along the diagonal of D, and the corresponding eigenvectors in the columns of V. For example, type the following: A = diag([1,4,6,8,-10,5,1,-15, -3, 12]); [V,D] = eigs(A,3,'lr'); V D It is clear that 12, 8, and 6 are the three largest eigenvalues of A, while the corresponding eigenvectors are e10 , e4 , and e3 . You can verify both in D and V. 2 (a) Use finite differences with 300 equally spaced points to formulate the following boundary value problem d2 u + kex u = u 0 < x < 1 ; dx2 u(0) = 0 ; u(1) = 0 , (5) as an eigenvalue problem of the form Au = u . (6) Give A and its dimensions. (b) For k = 1, find the three largest eigenvalues of Eqn. (6), and plot the corresponding eigenvectors (as a function of x) on one plot. (c) Repeat for k = 0. Can you guess what the exact answer should be? 3