E7 Fall 2019, UC Berkeley Homework Assignment 12 Ordinary Differential Equations The purpose of this lab is to introduce you to ODEs and ODE solvers in MATLAB. Note: Use the templates provided in the...

complete assignment and check with autogravder for 100 percent


E7 Fall 2019, UC Berkeley Homework Assignment 12 Ordinary Differential Equations The purpose of this lab is to introduce you to ODEs and ODE solvers in MATLAB. Note: Use the templates provided in the assignment download to complete this assignment. The template has a very specific format to enable auto-grading, and if the format is altered your assignment may not be properly graded (this will be apparent to you). The due time is 11:59pm (midnight) December 11, 2019. NO LATE HOME- WORK WILL BE ACCEPTED. Please upload the following files through the bCourses website: • HW12 ODE.m • HW12 ODE.pdf (created using the publish command) • All function files you are asked to create or are necessary to run your code. • All generated plot files specified in the questions. • Generated movie file specified. Directions to upload files can be found here. NOTE: This homework is designed in a sequential manner and gradually builds some con- cepts as the problems progress. Some latter problems might have dependencies on initial problems. 1. The Lorenz attractor, named after Edward N. Lorenz, was derived as a simplified model of convection currents in the atmosphere. It is known for having chaotic, unpredictable behavior. The Lorenz attractor can be described with the following three coupled ordinary differential equations: dx dt = σ(y − x) dy dt = x(ρ− z) − y dz dt = xy − βz 1 https://guides.instructure.com/m/4212/l/54353-how-do-i-upload-a-file-to-my-assignment-submission E7 Fall 2019, UC Berkeley where t is the time, x, y, and z are the state variables, and σ, β, and ρ are scalar parameters. σ is called the Prandtl number and ρ is called the Rayleigh number. Usually σ = 10, β = 8/3 and ρ is varied. The system exhibits chaotic behavior for ρ = 28, but displays knotted periodic orbits for other values of ρ. The ODEs that describe the Lorentz attractor can be written in vector form Ẏ = f(Y ) d dt  Y1Y2 Y3  =  σ (Y2 − Y1)Y1 (ρ− Y3) − Y2 Y1 Y2 − β Y3  (1) Here Y1 is x, Y2 is y and Y3 is z. Write a function named LorenzODE with the header line: function dotY = LorenzODE(Y, sigma, beta, rho) whose first input argument Y is a 3x1 column vector and sigma, beta, and rho are constants, and returns the right hand side of Eq. (1) (also a 3x1 column vector). Use ode45 to numerically integrate the ODE Eq. (1) for: σ = 10, β = 8/3, and ρ = 14, tspan = [0,20], and the following set of initial conditions: (a) Y0 = [0;1;1] (b) Y0 = [1;1;1] (c) Y0 = [1;-1;1] (d) Y0 = [1;-1;-1] (e) Y0 = [1;-1;-10] Remember that first input argument to ode45 should be an anonymous function handle, such as Lorenz = @(t,Y) LorenzODE(Y, sigma, beta, rho) You can use the commands for plotting: p lo t3 (Y( : , 1 ) ,Y( : , 2 ) ,Y( : , 3 ) ) ; view ( 0 , 0 ) ; hold on ; to do 3D plots of the solution (Y1(t), Y2(t), Y3(t)) for all five initial conditions on the same figure (using different values of color). Add the initial conditions as legend entries. Save the figure using print command as ’Lorenz.png’. 2. In this problem, you will use ode45 to simulate the motion of the double pendulum illustrated in the figure below. The equations of motion of the double pendulum are as follows: 2 E7 Fall 2019, UC Berkeley Figure 1: Double Pendulum System • θ̈ = −m2L1θ̇ 2 sin(θ−β) cos(θ−β)+gm2 sin(β) cos(θ−β)−m2L2β̇2 sin(θ−β)−(m1+m2)g sin(θ) L1(m1+m2)−m2L1 cos2(θ−β) • β̈ = m2L2β̇ 2 sin(θ−β) cos(θ−β)+g sin(θ) cos(θ−β)(m1+m2)+L1θ̇2 sin(θ−β)(m1+m2)−g sin(β)(m1+m2) L2(m1+m2)−m2L2 cos2(θ−β) Here θ and β are the angles of the two rods of the pendulum from the vertical con- figuration, as in the figure. Also, m1, m2 are the two masses (located at the ends of the two rods of length L1 and L2, respectively), and g is the gravitational acceleration. Assume that the preceding parameters take the following values: • Gravitational acceleration, g = 9.81 m·s−2 • Length of rod 1, L1 = 1 m • Length of rod 2, L2 = 1 m • Mass of joint 1, m1 = 1 kg • Mass of joint 2, m2 = 1 kg Use initial conditions defined as follows: • Angular position of rod 1, θ = 0.5 rad • Angular velocity of rod 1, θ̇ = 0 rad·s−1 • Angular position of rod 2, β = 0.5 rad • Angular velocity of rod 2, β̇ = 0 rad·s−1 Notice that there are two equations provided but have double derivatives (angular acceleration). We can, however, simplify them by identifying that θ̇ and β̇ can be set as two variables and write two equations for them, i.e. d dt θ = θ̇ 3 E7 Fall 2019, UC Berkeley d dt β = β̇ Thus the whole system can be written as Ẏ = f(Y ) as in the previous question. Start by creating an anonymous function, called dpend, for the double pendulum ODEs that you will subsequently use to call ode45. The function returns the value of f(Y) for an input of Y. Here Y = [θ; θ̇; β; β̇]. Use a time span of 0 to 20 seconds, but specify that the outputs of ode45 should be every 1 30 seconds. To do this define the time span as an evenly spaced vector from 0 to 20 with increments of 1 30 , i.e., define tspan = 0:1/30:20. Note when solving the problem, ode45 will still use an adaptive step size to control error, but will output values at the intervals specified by the time span vector. Define a function Y=DoublePendulum(I,t,C). Here I are the initial condition vector I = [θinitial; θ̇initial; βinitial; β̇initial], t is the time span vector tspan, C is the constant vector C = [L1;L2;m1;m2; g] and Y is the output vector with each column Yi defined as Yi = [θi; θ̇i; βi; β̇i] at time step ti. The total number of column vectors will be equal to the t vector, i.e. every column describes the angular position and angular velocities of the system. Remembering that θ and β are angular positions, use clf and drawnow commands to create an animation and save that animation as ’DoublePendulum.avi’ movie. Remem- ber to use suitable axis limits so the limits do not change from frame to frame. You have already done this for a previous homework. You can display the pendulums as ’o’ and rods with with dashed lines. 3. In this problem, you will come up with your own integration functions using Euler’s method and Runge-Kutta 4 (RK4) methods. Consider the following system of equa- tions: r̈ − rθ̇2 = − µ r2 rθ̈ + 2ṙθ̇ = 0 These are orbital equations. For simplicity, use µ = 1. Using same techniques as above come up with a system of equation in the form of Ẏ = f(Y ). Here order your equations so Y = [r; ṙ; θ; θ̇]. Now start by writing an anonymous function OrbitFunc(Y,mu) that takes as input the vector Y and outputs f(Y ). Once you have this function we can now solve them using two different techniques: 4 E7 Fall 2019, UC Berkeley (a) Euler’s Method: Write a function Yout = Euler(Yinitial, t,∆t, OrbitFunc) where Yinitial is your initial conditions. Use Yinitial = [2; 1.5; 1.5; 3], ∆t = 0.001 and t=5. Here Yinitial is the initial value of at t=0, t is the final time till which you need to advance your solution to, ∆t is the time step size and OrbitFunc returns f(Y). Inside your function, you will step forward in time using Euler’s method. The Euler’s equations, for value of Y at time tn+1, (where f(Y) is not a function of time), is defined as: Yn+1 = Yn + ∆tf(Yn) where Yn is the Y value at the previous time step. Define a vector tspan = 0 : ∆t : t. Your output Yout should be an array of Y where each column of Y is Y at that time step, i.e. number of columns of Y should be equal to the length of tspan. (b) RK4: Now lets use a much more stable time stepping scheme than Euler’s called Runge-Kutta 4 (RK4) method. Again a write a function Yout = RK4(Yinitial, t,∆t, OrbitFunc). Use the same inputs as before and your Yout would be similar as well as in the case of Euler’s. The method for time stepping will be changed however. RK4 method for value of Y at time tn+1, (where f(Y) is not a function of time), is defined as: Yn+1 = Yn + 1 6 (k1 + 2k2 + 2k3 + k4) where k1 = ∆tf(Yn) k2 = ∆tf(Yn + k1/2) k3 = ∆tf(Yn + k2/2) k4 = ∆tf(Yn + k3) Store your calculated variables Y from both methods as YEulerOrb and YRK4Orb. 4. In this problem, you will use ode45, Euler and RK4 to simulate the orbits of a planet around a (fixed) star and a moon around the planet, as illustrated in the figure below. The equations of motion for this problem are quite complex, so we will use instead a simplified set of equations to obtain a reasonable approximation of the system. The dependent variables of these approximate equations of motion are: the star-planet dis- tance R; the planet-moon distance r; the star-planet angle θ; and, the planet-moon angle β. The simplified equations of motion are second-order and can be expressed as: • R̈ = Rθ̇2 −GmSR−2 m·s−2 (radial acceleration of the planet) 5 E7 Fall 2019, UC Berkeley Figure 2: Moon-Planet-Star Orbit • θ̈ = −2Ṙθ̇R−1 rad·s−2 (angular acceleration of the planet) • r̈ = rβ̇2−GmP r−2−GmS(R+ r cos β)−2 m·s−2 (Radial acceleration of the moon) • β̈ = −2ṙβ̇r−1 rad·s−2 (angular acceleration of the moon) The constants that enter the preceding equations are: • Mass of the star, mS = 1 × 1030 kg • Mass of the planet, mP = 1 × 1023 kg • Mass of the moon, mM = 1 × 1022 kg • Gravitational constant, G = 6.673 × 10−11 m3·kg−1·s−2 Also, the initial conditions selected for the four variables are: • Radial position of the planet about the star, rP/S = 1 × 1011 m • Radial velocity of the planet, 0 m·s−1 • Angular position of the planet, 0 rad • Angular velocity of the planet, ωP/S = 2.25 × 10−7 rad·s−1 • Radial position of the moon about the planet, rM/P = 1 × 109 m • Radial velocity of the moon, 0 m·s−1 • Angular position of the moon, 0 rad • Angular velocity of the moon, ωM/P = 2 × 10−6