Let us now describe the motion of a double pendulum. The double pendulum is displayed in figure (figure from Wikipedia; credit to JabberWok) atitle The two masses are connected by rigid massless sticks. As a simplification, we also assume that the top stick is not connected to a ceiling, but to a nail (so it can go all the way around and we do not have to worry about collisions with the ceiling). Deriving the equation of motion is very tedious, so I save you the effort. We write the equation of motion for the angles, 1 and 2, between the two sticks and the vertical. The equations of motion for the two angles are dt2d21=L1(2m1+m2m2cos(2122))g(2m1+m2)sin1m2gsin(122)2sin(12)m2((dtd2)2L2+(dtd1)2L1cos(12))dt2d22=L2(2m1+m2m2cos(2122))2sin(12)((dtd1)2L1(m1+m2)+g(m1+m2)cos1+(dtd2)2L2m2cos(12)) The equations of motion need initial conditions. Let us assume for simplicity that the initial angular velocities are 0 , and that the initial angles are chosen by the user. Write a python function that takes as input: - mass of the mass 1m1 in kg - mass of the mass 2m2 in kg - length of the stick 1L1 in m - length of the stick 2L2 in m - initial value of angle 1(t=t0) in rad - mass of the mass 1m1 in kg - mass of the mass 2m2 in kg - length of the stick 1L1 in m - length of the stick 2L2 in m - initial value of angle 1(t=t0) in rad - initial value of angle 2(t=t0) in rad - initial time - final time - tolerance on the error on the angles The python function should output: - times at which the two angles have been calculated - angle 1 at the above times - angle 2 at the above times The python function should calculate the above quantities using adaptive time step RK4 (with step doubling). After constructing the python function above, use its output to plot the trajectories: - angle 1 in dependence of time - angle 2 in dependence of time - xy position of mass 1 (you need to figure out how it depends on 1 ) - xy position of mass 2 (you need to figure out how it depends on 1 and 2 ) You do not need to make an animation, just a plot. You can create an animation for your own extra training, but it will not be graded. Below you find the definitions of the functions you will need (notice that they have been slightly modified from class to make sure that they work with this problem). There is no need to program them again, just use them. import numpy as np import matplotlib.pyplot as plt def rungekutta4stepperError(operator, initCond, time, dt): Performs a step composed of 2 RK 4 steps an ODE of type dy(t)/dt=f(t,y(t)) dt=safetyFactordt(tolerance/error)(0.25) if dt=42102d1,I02d2,Iabs[L1sin(1,F(1,I,2,I))+L2sin(2,F(1,I,2,I))] Calculate this integral with the rectangle method. Below you find the definitions of the function you will need (notice that they have been slightly modified from class to make sure that they work with this Below you find the definitions of the function you will need (notice that they have been slightly modified from class to make sure that they work with this problem). There is no need to program them again, just use it. Notice that integrationRectangles2D expects the function to integrate to be in the form fun =lambdax,y: In [10]: def integrationRectangles2D(funcToUse, xmin, xmax, ymin, ymax, nSplitsX, nSplitsY): '' Calculates the integral of a 2D function by rectangle method ' ' dx=((xmaxxmin)SplitsX) dy=((ymaxymin)SplitsY) XRectMeth =np.linspace (xmin, xmax dx,nSplitsX) YRectMeth =np.linspace (ymin,ymaxdy,nSplitsY) XRectMeth, YRectMeth = np.meshgrid(XRectMeth, YRectMeth) XRectMeth = XRectMeth.flatten() YRectMeth = YRectMeth.flatten() totintegral =0 for x,y in zip(XRectMeth, YRectMeth): totIntegral += funcToUse (x,y) return dx * dy * totintegral Hints - The integration method needs a python function that gives the value of the mathematical function we want to integrate. This means that you need to construct it. Inside that you need to use the function you built in the previous homework, extract the result you need to produce the absolute value of the x position of the second mass and then return it. Build a function like def dist_from_axis(thetainitial1, theteInitial2): .. return calculatedDistance This will be the function that needs to be integrated. Don't get confused by the fact that it will not be an analytic expression: the only thing that matters here is that it gets two variables and returns a value
