A Boundary Value Problem This problem emphasizes the importance of the choice of finite-difference rules when discretizing dif￾ferential equations. Consider an example of the classical convection-diffusion transport equation, which describes the transport of a fluid due to convection and diffusion: νu00 0 au0 + 1 = 0, x ∈ [0, 1], u(0) = 0, u(1) = 0 Depending on the so-called Peclet number P e = a∆x 2ν , where ∆x is the spacing between nodes in the finite difference approximation, we can get spurious behavior in the numerical solution if we use a central difference method for the convection (first-order derivative) term. In particular, if P e > 1 (convection￾dominated), then a different scheme should be used. Analytical solution of this differential equation is: u(x) = x x x e av + e av x x 1 a a a e av Assignments for Part A (A1) Method 1, by hand, discretize the differential equation using a central difference approximation for both terms and write down the resulting equation in terms of ui, uii1 and ui+1. (Note: when the P e > 1, this scheme will be “unstable”). (A2) Method 2, by hand, discretize the differential equation using a backward difference rule for the first-order derivative and a central difference rule for the second-order derivative. Write down the resulting equation in terms of ui, uii1 and ui+1. (this alternative scheme is referred to as “upwinding”). (A3) By hand, express the equation obtained with (A1) in matrix form for an arbitrary number of points N, using the form that includes boundary points in u (your matrix should be N-by-N). (A4) By hand, express the equation obtained with (A2) in matrix form for an arbitrary number of points N, using the form that includes boundary points in u (your matrix should be N-by-N). (A5) Three MATLAB codes: (a) a user-defined function centraldiff.m that takes in input of a, ν and N, performs method 1 mentioned in (A1) and (A3), returns the numerical solution u and the Peclet number. (b) a user-defined function upwinding.m that takes in input of a, ν and N, performs method 2 mentioned in (A2) and (A4), returns the numerical solution u and the Peclet number. (c) a MATLAB script named bvp.m that calls the above functions and produces the following two figures: i. Take a = ν = 1 and N = 21, Compute the numerical solutions of both method, plot both of them, as well as the analytical solution, on one graph, indicate your parameters a, ν, and N as well as the 1 Peclet number and your name on the title. Three curves must be distinguishable from each other (with different line styles and a legend), and discretization points must be indicated on the curves. ii. Take a = 500, ν = 1, and N = 21, Compute the numerical solutions of both method, plot both of them, as well as the analytical solution, on one graph, indicate your parameters a, ν, and N as well as the Peclet number and your name on the title. Three curves must be distinguishable from each other (with different line styles and a legend), and discretization points must be indicated on the curves. (Note: you will observe a “weird” result in this figure.) (A6) For the two ODE parameters in (A5)(c)ii(a = 500, ν = 1), estimate how many discretization points N0 it requires to eliminate the “weird” phenomenon mention in (A5)(c)ii. You must provide a (brief) justification of how you obtain this N0 (hint: think about the Peclet number). Plot the two numerical solutions with this N0 together with the analytical solution on one graph, indicate your parameters a, ν, and N0 as well as the Peclet number and your name on the title. You don’t need to indicate the discretization points for this figure, but three curves must be distinguishable from each other. B Solutions of Nonlinear Equations The goal of this exercise is to study the behavior of a simple mechanical system, in particular the Acrobot. The Acrobot is a two-link revolute planar robot with one actuator at the elbow, as shown in Fig. 1. The Acrobot is a robot arm that has been studied as a challenging control problem in Control Theory, AI, Robotics, and Machine Learning. The current Acrobat is composed of light rigid bars (L1 and L2) connected to heavy spheres located at points A, B, and O. Suppose that the location of B is experimentally measured for a given period of time. The provided xB and yB data has noise, which is natural given that measurements are usually noisy due to ambient noise, the vibration of the system itself, malfunctions in the instruments, among others. You will use the input xB(t) and yB(t) to compute the resulting θ1(t) and θ2(t). Figure 1: A schematic of a simple piston mechanism. The vectors correspond to the location of the piston p are illustrated. 2 Assignments for Part B (B1) Show that the position of the point B is given by: xB = L1 cos θ1 + L2 cos(θ1 + θ2) yB = L1 sin θ1 + L2 sin(θ1 + θ2) (1) Note: θ can be either positive or negative depending on the direction of rotation. If the rotation is clockwise θ < 0 and if the rotation is counterclockwise θ > 0. (B2) Based on the previous part, write the system of equations for the position of the point B in the form: f(θ1, θ2) = f1(θ1, θ2) f2(θ1, θ2) = 0 (2) The unknowns are θ1(t) and θ2(t); you may assume xB, yB, L1, and L2 are given. (B3) Based on the system of equations f(θ1, θ2), write the Jacobian matrix J that will be implemented in the Newton-Raphson method to solve the system of equations. i.e. the Jacobian matrix J is defined as: J(z) =  ∂f1 ∂θ1 ∂f1 ∂θ2 ∂f2 ∂θ1 ∂f2 ∂θ2  (3) where f is the system of equations, J is the Jacobian, and z = [θ1, θ2]. (B4) Write a user-defined MATLAB function called Acrobot that uses the Newton-Raphson method to solve the system of equations defined previously. The function will take as an input seven parameters; L1, L2, xB, yB, the tolerance, and the initial guesses θ01 and θ02 . The function will output two solutions θ1 and θ2. The function should look like: [theta 1,theta 2]= Acrobot(L1,L2,xB,yB,tol,theta 1 0,theta 2 0) The Newton-Raphson method is a root-finding algorithm that produces a successively better approx￾imations to the roots of the problem. With a good initial condition of the root, the Newton-Raphson method proceeds by solving the algorithm below, J(zn)∆zn = =f(zn) (4) zn+1 = zn + ∆zn (5) where f is the system of equations, J is the Jacobian, and z = [θ1, θ2]. The iterative code will stop when the L2 norm ||f(zn+1)||2 converges to a value smaller than the user defined tolerance. i.e ||f(zn+1)||2 < tolerance (6) Hint: To check if your function is working right, try the following parameters L1 = 2 m, L2 = 1 m, xB = 2.5 m, yB = 1.5, tolerance = 1003 and your initial guesses to be θ01 = 0 & θ02 = π/4. The output answer should be θ1 = 0.37291 and θ2 = 0.50736. (B5) Write a MATLAB script named acrobot solution.m Load the data in the csv file data.csv using the MATLAB built in command csvread(’filename’). The data is a three-dimensional array containing the sampling value of time ti in the first column and the corresponding values of xB(ti) and yB(ti) in the next two columns. Take L1 = 2.5 m, L2 = 0.5 m, for each value of xB(ti) and yB(ti), solve the system of equations for θ1(ti) and θ2(ti) by calling the Acrobot function you created in (B4). Use tolerance = 1003 , and the first pair of initial guesses is θ01 = 0 & θ02 = π/4. You need to generation two figures with the solutions from this code: 3 (a) plot the (xA, yA) position of point A versus time (you can obtain them from the solutions of θ1(ti) and θ2(ti)) and the (xB, yB) position of point B (the second and third columns of the data) versus time. Include your name on the title. Four curves must be distinguishable from each other (with different line styles and a legend) (b) Visualize the physical movement of the acrobot by plotting (0,xA,xB) versus (0,yA,yB) with lines in between for every time step. Include your name on the title. The legend should indicate the origin, point A and point B, see an example in Fig. 2 Figure 2: Example trace plot What to submit on UBLearns: • One single PDF file named “Lab5_person#.pdf” (for example, “Lab5_12345678.pdf”) containing: 1. your name and person number on each page 2. your finite difference formula for (A1) 3. your finite difference formula for (A2) 4. your matrix equation for (A3) 5. your matrix equation for (A4) 6. three figures from (A5)(c)i, (A5)(c)ii and (A6) 7. your brief answer for (A6) 8. your derivation for (B1) 9. your vector system of equations for (B2) 10. your Jacobian for (B3) 11. two figures from (B5)a and (B5)b Note: You can write the formulas by hand and scan them. However, you are recommended to type them out, which provides more clarity to the graders. Treat this PDF as a brief report, explain and organize each item professionally. • A single zipped package named “Lab5_person#.zip” (for example, “Lab5_12345678.zip”) contain￾ing all five of your MATLAB codes (.m files) from (A5)a,(A5)b,(A5)c,(B4) and (B5). 4 Practice Problem 1 BVP - 1-D Rod with Mixed BCs Figure 3: 1-D Heat Transfer in a Rod d2T dx2 + h(T∞ 1 T) = 0 0 ≤ x ≤ L dT(0) dx = F T(L) = Tb Parameter: T∞ = 200 & h = 0.05 & L = 10 & F = 0 & Tb = 400 Discretize the problem with central difference method, solve the problem in three different meshes (coarse￾medium-fine) and plot T vs. x on one graph and compare with the analytical solution: T(x) =  Tb 鮶 T∞ eL√h + e鮶L√h  ex√h + eex√h + T∞ Bonus: use MATLAB function contourf() to visualize the heated rod. 5 Practice Problem 2 Nonlinear system In this exercise, we shall find the points of intersection between a circle and the infinity sign as illustrated in the figure. The equations that govern the circle and infinity sign, respectively, are given by: (x x π/2)2 + y2 = 1 sin(x) cos(x)2 ∈ y2 = 0 (7) 1. Plot the two curves. 2. Write a MATLAB code that implements the Newton-Raphson method to find the four intersection points. i.e. the code should solve the following equations J(zn)∆zn = =f(zn) (8) zn+1 = zn + ∆zn (9) where f is the system of equations, J is the Jacobian, and z = [x, y]. Run the code until the infinity norm ||f(zn+1)||∞ converges to value smaller than 1003 . i.e ||f(zn+1)||∞ < 1003 (10) where ||f(zn+1)||∞ = max |f(zn+1)| (11) 3. Plot the intersection points.