doublePendulum v0.6 written by Alexander Erlich (mailto:alexander.erlich@gmail.com) Last update submitted to MATLAB File Exchange on October 5th, 2010. Lagrangian and di\u001berential equations The program double_pendulum.m animates the double pendulum's (mostly) chaotic behavior. The Lagrangian of the- double pendulum as depicted MACM 316 Assignment # 8 in \u001cgure 1 is m1 + m2 2 2 m2 2 2 l1 1 + l + m2 l1 l2 1 2 cos (1 2 ) + 2 2 2 2 + (m 1 +your m2 g report l2 cos 2(to Assignment 1 + m2 ) g l1 cosand You must upload both your code (to Assignment #8 scripts/codes) L Due Date: April 5th, 2017, at 11pm. = #8 computing report). The assignment due at havethe setEuler-Lagrange the due timeequations in Canvas to The equations ofismotion can11:00pm. be derivedIusing 11:05pm and if Canvas indicates that you submitted late, you will be given 0 on the assignment. L d L Your computing report must be exactly 1 page. There will be a penalty , igiven = 1, 2if your report is dt i i longer than one page. One obtains two ordinary di\u001berential equations of second order: Please read the Guidelines for Assignments first. \u0003 g m l \u0002 2 2 sin 1 + cos (2 1 ) 2 sin (2 1 ) 22 = l1 open forums. m1 + m2 l1 Keep in mind that Canvas discussions are \u0003 g l1 \u0002 2 + cos (2 1 ) 1 + sin (2 1 ) 21 = sin 2 + l Acknowledge any collaborations and assistance lfrom colleagues/TAs/instructor. 2 2 1 + 0 0 It is now possible to rewrite this system of two second order ODEs into a system of four \u001crst ord ODEs. De\u001cning e.g. the x1 =double 1 and xpendulum 1 = x 1 as well as x3 = 2 and x4 = 2 = x 3 and 2 = Computing Assignment - Simulating troducing a vector x = (x1 , x2 , x3 , x4 )T one obtains a system x = f (x) of \u001crst order ODEs. T results of such a manipulation are presented in double_pendulum_ODE.m. A Mathematica notebo This assignment requires you to do quite a bit of coding onallyour own soismake sure double_pendulum_ODE_deduction.nb containing manipulations also provided. that you start early. Also, you need to present a lot of results in your report so think about how you are going to organize the report. Do not use large figures but make Running double_pendulum sure that the figure fonts are large enough for the marker to be able to read them. The most simple way to run the program is >\u0017> double_pendulum_init. The parameters 1 , 1 , 2 , In this assignment you will simulate the l2motion of a double pendulum i.e. a pendulum g, m can be adapted in the double_pendulum_init.m \u001cle. Thewith parameters duration, 1 , m2 , l1 and and movie de\u001cne the duration and framerate of the animation and whether the animation is suppos another pendulum attached to its end. The figure below shows a schematic depiction of the setup be shown in realtime or rendered into assignment. a movie (.avi \u001cle). and introduces our notationtofor the different variables in this See Further and more technical information is given in the comments preceding the source code of the m-\u001cl https://www.myphysicslab.com/pendulum/double-pendulum/double-pendulum-en.html for an animation that you can play around with. y 3 x 1 g= l1 \u0012 0 g \u0013 2 1 0 -1 m1 2 l2 -2 m2 -3 -3 -2 -1 0 1 Figureof1:the Scheme drawing (left) and Matlab \u001cgure (right) Figure 1: Schematic drawing double pendulum. Assuming that the lengths l1 and l2 are constant and the connecting strings are massless one can show that the motion of the pendulums is governed by the following system of second order ODEs: \u0003 m2 l2 \u0002 g 001 (t) + sin(1 (t)) + cos(2 (t) 1 (t))002 (t) sin(2 (t) 1 (t))0 2 (t)2 = 0, l1 m1 + m2 l1 (?) \u0003 g l1 \u0002 00 00 0 2 2 (t) + sin(2 (t)) + cos(2 (t) 1 (t)) 1 (t) + sin(2 (t) 1 (t)) 1 (t) = 0. l2 l2 1 2 3 In order to solve this system we introduce a new variable x(t) = (x1 (t), x2 (t), , x9 (t))T (column vector) where x1 (t) = 1 (t), x5 (t) = g, x2 (t) = 0 1 (t), x6 (t) = m1 , x3 (t) = 2 (t), x7 (t) = m2 , x4 (t) = 0 2 (t), x8 (t) = l1 , x9 (t) = l2 . After some algebraic manipulation (which is not directly related to this assignment) we can rewrite (?) as a system of first order ODEs of the form x(t) = f (t, x(t)). (??) The script double pendulum ODE.m contains a function that evaluates the right hand side f at a given value of t and x. (a) Your first task in the assignment is to approximate the solution of this ODE system using the Euler's method for 0 t tend = 10. Use the initial condition x(0) = (/6, 0, /6, 0, 9.8, 2, 1, 2, 1)T and plot 1 (t) and 2 (t) (that is x1 (t) and x3 (t)) as a function of time using time steps of size h = 101 , 103 . Make sure that you get the code runing with the smaller step size and a short time interval. Describe your obersvations and comment on robustness of your code with respect to the step size h. I suggest you write your code by completing the script double pendulum skeleton.m that contains the skeleton of the code that you need to write. Note that there are instructions in the skeleton code as well. (b) Your next task is to repeat the above calculation using MATLAB's ode45. I suggest that you begin by reading the manual page for the ode45 function (just google \"MATLAB ode45\" or search in MATLAB's help documentation). You should use the same initial conditions and time interval as in part (a). However, ode45 does not require you to choose a time step. Plot 1 (t) and 2 (t) as before and compare your solution to those obtained using Euler's method. Which method is more accurate? which one is more efficient? (c) The kinetic and potential energies of the double pendulum are given by the following expressions: 1 1 K(t) = m1 01 (t)2 l12 + m2 [01 (t)2 l12 + 02 (t)2 l22 + 201 (t)l1 02 (t)l2 cos(1 (t) 2 (t))] 2 2 P (t) = m1 l1 g cos(1 (t)) m2 g[l1 cos(1 (t)) + l2 cos(2 (t))] Repeat the calculations in parts (a) and (b) for tend = 200 and plot the total enery of this system E(t) = K(t) + P (t) for the solution of the ode45 solver and your Euler's solver with h = 103 . Describe your observations. Which method is more trustworthy? 2 %% this is a fully worked out script for this assignment %% MACM 316 - Spring 2017, Brenda Davison %% Bamdad Hosseini clear all close all clc % initiate the problem % ivp=[phi1; phi1_dot; phi2; phi2_dot; g; m1; m2; l1; l2] - vector of % parameters and solution of ODE % tend - termination time of simulations x0 = [pi/6; 0; pi/6; 0; 9.8; 2; 1; 2; 1]; tend = 10; %% ODE45 solution - determine for yourself how to call ODE45 %[t45, sol_ode45] = ode45(...); % plot solution using ODE45 figure(1001); %improve the title so your figures are clearly labelled. plot(t45, sol_ode45(:,1), 'b', t45, sol_ode45(:,3), 'r' ); xlabel('Time (t)') ylabel('angle') legend('\\phi_1', '\\phi_2') %% Forward Euler solution h = 0.1; sol_Euler(:,1) = x0; time = 0:h:tend; % Put code here to solve using forward Euler method. figure(1002); %improve the title so your figures are clearly labelled. plot(time, sol_Euler(1,1:end-1), 'b', time, sol_Euler(3,1:end-1), 'r' ); xlabel('Time (t)') ylabel('angle') legend('\\phi_1', '\\phi_2') % compute total energy of solutions m1 = x0(6); m2 = x0(7); l1= x0(8); l2=x0(9); g = x0(5); % Here you need to correctly compute the potential and kinetic energy with % your solution from ode45 %P_ode45 = ... %K_ode45 = ... figure(2001) %improve the title so your figures are clearly labelled. plot(t45, P_ode45 + K_ode45) xlabel('t') ylabel('Total energy') title('ode45') % Here you need to correctly compute the potential and kinetic energy with % your solution from forward Euler. %P_Euler = ... %K_Euler = ... figure(2002) %improve the title so your figures are clearly labelled. plot(time, P_Euler(1:end-1) + K_Euler(1:end-1)) xlabel('t') ylabel('Total energy') title('Eulers method') %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% MACM 316 - Spring 2017 %% Brenda Davison %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % this function defines the double pendulum ODE function xdot = double_pendulum_ODE(t,x) % % % % % % % % % author: Alexander Erlich (alexander.erlich@gmail.com) parameters: t xdot Column vector of time points Solution array. Each row in xdot corresponds to the solution at a time returned in the corresponding row of t. g=x(5); m1=x(6); m2=x(7); l1=x(8); l2=x(9); xdot=zeros(9,1); xdot(1)=x(2); xdot(2)=-((g*(2*m1+m2)*sin(x(1))+m2*(g*sin(x(1)-2*x(3))+2*(l2*x(4)^2+... l1*x(2)^2*cos(x(1)-x(3)))*sin(x(1)-x(3))))/... (2*l1*(m1+m2-m2*cos(x(1)-x(3))^2))); xdot(3)=x(4); xdot(4)=(((m1+m2)*(l1*x(2)^2+g*cos(x(1)))+l2*m2*x(4)^2*cos(x(1)-x(3)))*... sin(x(1)-x(3)))/(l2*(m1+m2-m2*cos(x(1)-x(3))^2))