Download Runge–Kutta methods HW and more Assignments Physics in PDF only on Docsity!
HOMEWORK 2, PHYS 562
Kam Modjtahedzadeh
October 15, 2020
Abstract The Runge–Kutta methods have been used to solve three different problems: 1) The double-pendulum. 2) The circular restricted 3 body problem. 3) The spring pendulum with a ring suspension. The equations of motion have been found for each individual problem and have been included in scripts used to supplement the Runge–Kutta code. All scripts have been primarily written in GFortran; a Computer Algebra System is used for the heavy mathematics and Python was used for plotting.
1 Double Pendulum^1
Figure 1: Image from Google.^2 Here the first black ball from the top has mass m 1 and the bottom black ball has mass m 2.
1.1 Equations of Motion
Let r 1 = (x 1 , y 1 ) be the position of m 1 and r 2 = (x 2 , y 2 ) be the position of m 2. To find the equations of motion I will treat the masses as point particles and start with writing Newton’s equations for net force ΣF = m¨r. Starting with m 1 ,
m 1 x¨ 1 = T 2 sin θ 2 − T 1 sin θ 1 (1a) m 1 y¨ 1 = T 1 cos θ 1 − T 2 cos θ 2 − m 1 g. (1b)
where T 1 and T 2 are the tension forces on l 1 and l 2 respectively. Also writing the same equations for m 2 :
m 2 x¨ 2 = −T 2 sin θ 2 (2a) m 2 y¨ 2 = T 2 cos θ 2 − m 2 g. (2b)
(^1) I approached this problem using Newtonian mechanics. It is much easier and perhaps better to use Lagrangian mechanics [4]. (^2) Credit to author(s)/artists(s). I forgot to check the exact website I found this image from.
Rewriting the equations in (2) for T 2 terms and substituting them into the equations for (1):
(m 1 x¨ 1 = −m 2 x¨ 2 − T 1 sin θ 1 ) × cos θ 1 (3a) (m 1 ¨y 1 = T 1 cos θ 1 − m 2 y¨ 2 − m 2 g − m 1 g) × sin θ 1. (3b)
where I have taken the liberty of also multiplying the equations in (3) by cos θ 1 and sin θ 1 respectively. This gives me the following subset of equations:
T 1 sin θ 1 cos θ 1 = −(m 1 x¨ 1 + m 2 ¨x 2 ) cos θ 1 (4a) T 1 sin θ 1 cos θ 1 = (m 1 y¨ 1 + m 2 y¨ 2 + m 2 g + m 1 g) sin θ 1. (4b)
Now equating the two equations above:
− (m 1 ¨x 1 + m 2 x¨ 2 ) cos θ 1 = (m 1 y¨ 1 + m 2 y¨ 2 + m 2 g + m 1 g) sin θ 1. (5)
Furthermore, multiplying the equations in (2) by cos θ 2 and sin θ 2 respectively I get,
T 2 sin θ 2 cos θ 2 = −m 2 x¨ 2 cos θ 2 (6a) T 2 sin θ 2 cos θ 2 = (m 2 y¨ 2 + m 2 g) sin θ 2 (6b)
− m 2 x¨ 2 cos θ 2 = (m 2 y¨ 2 + m 2 g) sin θ 2 , (7)
where in the last step I just equated the equations in (6). I also have the following equations:
x 1 = l 1 sin θ 1 y 1 = −l 1 cos θ 1 x 2 = x 1 + l 2 sin θ 2 y 2 = y 1 − l 2 cos θ 2.
So taking the second derivative (8):
¨x 1 = − θ˙^21 l 1 sin θ 1 + ¨θ 1 l 1 cos θ 1 ¨y 1 = θ˙^21 l 1 cos θ 1 + θ¨ 1 l 1 sin θ 1 ¨x 2 = − θ˙^21 l 1 sin θ 1 + ¨θ 1 l 1 cos θ 1 − θ˙^22 l 2 sin θ 2 + θ¨ 2 l 2 cos θ 2 ¨y 2 = θ˙^21 l 1 cos θ 1 + θ¨ 1 l 1 sin θ 1 + θ˙^22 l 2 cos θ 2 + θ¨ 2 l 2 sin θ 2.
Now plugging in these values into equations (5) and (7);
m 1
− θ˙^21 l 1 sin θ 1 + θ¨ 1 l 1 cos θ 1
− θ˙^21 l 1 sin θ 1 + θ¨ 1 l 1 cos θ 1 − θ˙ 22 l 2 sin θ 2 + ¨θ 2 l 2 cos θ 2
cos θ 1 = ( m 1
θ˙^21 l 1 cos θ 1 + θ¨ 1 l 1 sin θ 1
θ˙ 12 l 1 cos θ 1 + ¨θ 1 l 1 sin θ 1 + θ˙^22 l 2 cos θ 2 + θ¨ 2 l 2 sin θ 2
sin θ 1 (10)
−m 2
− θ˙^21 l 1 sin θ 1 + θ¨ 1 l 1 cos θ 1 − θ˙^22 l 2 sin θ 2 + θ¨ 2 l 2 cos θ 2
cos θ 2 = ( m 2
θ˙^21 l 1 cos θ 1 + θ¨ 1 l 1 sin θ 1 + θ˙^22 l 2 cos θ 2 + θ¨ 2 l 2 sin θ 2
sin θ 2. (11)
22 / ( l 2 ∗ ( 2 ∗m1+m2−m2∗ c o s ( 2 ∗ y ( 1 ) − 2 ∗y ( 2 ) ) ) ) 23 24 k ( 1 : n_eq ) = h∗ f ( 1 : n_eq ) 25 end f u n c t i o n kv 26 end s u b r o u t i n e r k f 4 5 s t e p
The script above is the function that goes into the RKF45 subroutine. Lines 2-4 have been excluded to save space.^5 m1, m2, l1, and l2 are the constants defined the setup module (see following script). y(1:4) are θ 1 , θ 2 , ω 1 , and ω 2 respectively. f(1:4) are the equations of motion (12)–(15) respectively.
The constants are: l 1 = 0.5 m, l 2 = 0.4 m, m 1 = 3 kg, and m 2 = 2 kg. First I will the use the initial conditions θ 1 = θ 2 = π/ 10 , π/ 6 , π/ 4 , π/ 2 as well as ω 1 = ω 2 = 0. I will start with time t = 0 and go up to 200 seconds.
0 25 50 75 100 125 150 175 200 t
2
1
0
1
2
(a) θ 1 = θ 2 = π/ 10
0 25 50 75 100 125 150 175 200 t
4
3
2
1
0
1
2
3
4
(b) θ 1 = θ 2 = π/ 6
0 25 50 75 100 125 150 175 200 t
6
4
2
0
2
4
6
(c) θ 1 = θ 2 = π/ 4
0 25 50 75 100 125 150 175 200 t
(d) θ 1 = θ 2 = π/ 2
Figure 2: ω vs t. The blue and orange colors represent the data for ω 2 and ω 1 respectively.
Next, I will use the initial conditions θ 1 = θ 2 = 0, ω 1 = π/ 10 , and ω 2 = π/ 2. The plots corresponding to these initial conditions are in Figure 3. The external code used to support the RKF45 script is included below.
1 module s e t u p 2 u s e numtype 3 i m p l i c i t none 4 5 i n t e g e r , parameter : : n_eq = 4 6 r e a l ( dp ) , parameter : : g = 9. 8 1 _dp , m1 = 3. _dp , m2 = 2. _dp , & 7 l 1 = 0. 5 _dp , l 2 = 0. 4 _dp 8 end module s e t u p 9 10 program pendulum 11 u s e s e t u p 12 i m p l i c i t none
(^5) These are the same lines included in Professor Zoltan Papp’s code from class. Nothing has been modified.
0.3 0.2 0.1 0.0 0.1 0.2 0. 2
2
(a) ω 2 vs θ 2 resembles a torus.
1.00 0.75 0.50 0.25 0.00 0.25 0.50 0.75 1. 1
2
(b) ω 2 vs ω 1 resembles a sheet.
Figure 3: Special case where I start with θ 1 = θ 2 = 0, ω 1 = π/ 10 , and ω 2 = π/ 2. Subfigure 3a represents the angular phase space of the bottom pendulum with the given initial conditions. Subfigure 3b plots the angular velocities of the two pendulums against each other.
13 14 r e a l ( dp ) : : t , dt , tmax , t h e t a ( 0 : 3 ) 15 r e a l ( dp ) , d i m e n s i o n ( n_eq ) : : y 16 i n t e g e r : : i 17 18 t h e t a = ( / p i / 1 0. _dp , p i / 6. _dp , p i / 4. _dp , p i / 2. _dp/ ) 19 20 do i = 0 , 3! l o o p i n g a l l 4 d i f f e r e n t ICs 21 t = 0. _dp! r e s e t t i n g y , dt , tmax each time 22 dt = 1. 5 _dp 23 tmax = 2 0 0. _dp 24 25 y ( 1 ) = t h e t a ( i )! t h e t a 1 26 y ( 2 ) = t h e t a ( i )! t h e t a 2 27 y ( 3 ) = 0! omega 1 28 y ( 4 ) = 0! omega 2 29 30 do w h i l e ( t < tmax ) 31 w r i t e ( i +13 ,∗) t , y ( 1 ) , y ( 2 ) , y ( 3 ) , y ( 4 ) 32 c a l l r k f 4 5 s t e p ( t , dt , y ) 33 end do 34 end do 35 36 t = 0. _dp! r e s e t t i n g them a g a i n 37 dt = 1. 5 _dp 38 tmax = 2 0 0. _dp 39 40 y ( 1 ) = 0! new ICs 41 y ( 2 ) = 0 42 y ( 3 ) = p i / 1 0. _dp 43 y ( 4 ) = p i / 2. _dp 44 45 do w h i l e ( t < tmax ) 46 w r i t e ( 2 0 , ∗ ) t , y ( 1 ) , y ( 2 ) , y ( 3 ) , y ( 4 ) 47 c a l l r k f 4 5 s t e p ( t , dt , y ) 48 end do 49 end program pendulum
Line 18 has the 4 various initial conditions for θ 1 and θ 2 , and the lines 20–34 include an algorithm to write the data separately for each initial condition. After once again resetting the time parameters and initial conditions lines 45–49 output data for a certain case of θ 1 = 0 and ω 1 > 0.
so with these new masses I now have a normalized relationship between the two bodies;
m 2 + m 1 = 1. (26) To accommodate for the stationary configurationsthe two primary large masses have fixed positions in a co-rotating frame with the origin at the center of mass. To do this first consider the derived formula for angular velocity Ω from Kepler’s third law;^7
Ω =
G(m 2 + m 1 ) r^312
where r 12 is the distance between the two primary masses [6]. Now that I am working in a non-inertial frame of reference, we know that for uniform circular motion the acceleration is,
a = v^2 r
= Ω^2 r, (28)
but since it is a rotating frame, the velocity changes by ∆v at each rotation. Therefore,
a =
(v + ∆v)^2 r
Ω^2 r^2 + 2Ωr∆v + ∆v^2 r
= Ω^2 r + 2Ω∆v + ∆v^2 r
≈ Ω^2 r + 2Ω∆v, (32)
where the last term in line (31) has been ignored because it is small relative to the other two terms.
Furthermore, in the current defined frame of reference the following are the distances squared between the smaller mass m 3 and the larger masses m 1 and m 2 respectively.
r^21 = (x + μ)^2 + y^2 (33) r^22 = (x − 1 + μ)^2 + y^2 (34)
Now I will rewrite equation (21),
x ¨ − 2Ω ˙y − Ω^2 x = −G
m 2 r^32 (x − 1 + μ) − G
m 1 r^31 (x + μ) (35)
y ¨ − 2Ω ˙x − Ω^2 y = −G
m 2 r 23
y − G
m 1 r^31
y, (36)
where I have broken up #»r into its Cartesian coordinates (x, y). To finalize the normalization of (21) I set G = Ω = 1. I also substitute m 2 and m 2 in the above equations with the expressions in lines (25) and (24) respectively.
x ¨ − 2 ˙y − x = −
μ(x − 1 + μ) r^32
(1 − μ)(x + μ) r^31
y ¨ − 2 ˙x − y = −
μy r^32
(1 − μ)y r^31
2.2 The Lagrange Points
Since I am trying to solve for the stationary configurations, the velocity and acceleration are assumed to be zero.^8 This simplifies equations (37) and (38) to the following,
x = μ(x − 1 + μ) r^32
(1 − μ)(x + μ) r^31
y =
μy r^32
(1 − μ)y r^31
(^7) Admittedly I do not fully understand Kepler’s law; I adopted this equation and its use here from Nkosi N. Trim [6]. (^8) The velocity of m 3 is not necessarily zero.
The first three Lagrange points lie on the same line connecting m 1 and m 2 , so their y components are zero. Their x components are found by solving equation (39) using a CAS. These points are:
L 1 = (+0. 9900271 , 0)
L 2 = (+1. 0100336 , 0)
L 3 = (− 1. 0000013 , 0),
note that from (23) μ is found to be
R⊕
R⊕ + R
≈ 3. 003 × 10 −^6.
The stable Lagrange points L 4 and L 5 can be obtained by remembering that they are the third vertices of equilateral triangles, the other two being the m 1 and m 2. The points are,
L 4 = (0. 5 − μ,
L 5 = (0. 5 − μ, −
2.3 RK4 Method
To use the results found above in order to find the two-dimensional motion of the comet m 3 , I will utilize the classic Runge–Kutta method (RK4). The equations of motion for the RK4 are,
vx = ˙x (41) vy = ˙y (42)
v ˙x = − μ(x − 1 + μ) r 23
(1 − μ)(x + μ) r^31
v ˙y = −
μy r 23
(1 − μ)y r^31
The Fortran 95 code for the above equations is the following,
1 s u b r o u t i n e r k 4 s t e p ( x , h , y ) 2!.. 3!.. 4!.. 5 c o n t a i n s 6 f u n c t i o n kv ( t , dt , y ) r e s u l t ( k ) 7 u s e s e t u p 8 i m p l i c i t none 9 r e a l ( dp ) , i n t e n t ( i n ) : : t , dt 10 r e a l ( dp ) , i n t e n t ( i n ) , d i m e n s i o n ( n_eq ) : : y 11 r e a l ( dp ) , d i m e n s i o n ( n_eq ) : : f , k 12 r e a l ( dp ) : : r1 , r 2 13 14 r 1 = s q r t ( ( y ( 1 ) −(1−mu) ) ∗∗ 2 + y ( 2 ) ∗ ∗ 2 ) 15 r 2 = s q r t ( ( y ( 1 )+mu) ∗∗ 2 + y ( 2 ) ∗ ∗ 2 ) 16 17 f ( 1 : 2 ) = y ( 3 : 4 ) 18 19 f ( 3 ) = −(mu∗ ( y ( 1 )−1+mu) ) / r 1 ∗∗ 3 − ((1 −mu) ∗ (mu+y ( 1 ) ) ) / r 2 ∗∗ 3 & 20 + y ( 1 ) + 2 ∗ y ( 4 ) 21 22 f ( 4 ) = − (mu∗y ( 2 ) ) / r 1 ∗∗ 3 − ((1 −mu) ∗y ( 2 ) ) / r 2 ∗∗ 3 + y ( 2 ) − 2 ∗ y ( 3 ) 23 24 k = dt ∗ f 25 end f u n c t i o n kv 26 end s u b r o u t i n e r k 4 s t e p
where similar to the RKF45 listing in section 1.3, I have excluded lines 2–4 to save space. The setup module as well as the program which included the initial conditions used to supplement the RK4 method are included below.
1.00 0.75 0.50 0.25 0.00 0.25 0.50 0.75 1. x
y
(a) L 1 , vx = 0, t from 0 to 90
1.00 0.75 0.50 0.25 0.00 0.25 0.50 0.75 1. x
y
(b) L 1 , vx = 0, t from 0 to 300
1.0 0.5 0.0 0.5 1. x
y
(c) L 1 , vx = 0, t from 0 to 200
1.0 0.5 0.0 0.5 1. x
y
(d) L 3 , vx = 0. 1 , t from 0 to 321
1.00 0.75 0.50 0.25 0.00 0.25 0.50 0.75 1. x
y
(e) L 4 , vx = 0. 1 , t from 0 to 90
1.00 0.75 0.50 0.25 0.00 0.25 0.50 0.75 1. x
y
(f) L 4 , vx = 0, t from 0 to 400
Figure 6: y vs x plot for various initial condition for different Lagrange points. This is the motion of a comet m 3 in a CR3BP with the other two masses in this instant being the Earth and the Sun. vy starts at 0 for all the above cases.
3 Driven Spring Pendulum on Ring^9
In this problem I will attempt to find the equations of motion for a spring pendulum which with angular velocity Ω along a circular path driven with radius R.^10 To solve this problem I will use Lagrangian mechanics. (^9) Unlike the previous two sections, I did not break this one up into subsections. This is partially because I ‘did it all in one step,’ partially because I adapted my energies from Starosta et al., 2008 [7], and partially because of a lack of time. (^10) I used Ω instead of ω to avoid confusion with Figure 7.
Figure 7: Image from Starosta et al., 2008 [7]. Here, xst = ∆L, X and ϕ are the generalized coordinates, where the position of the mass m is x = L 0 + ∆L + X(t) = L + X(t). Variables such as the moment M are not used in this problem. Not that ϕ is the angle made by the spring pendulum and the vertical axis; its derivative ω should not be mistaken with Ω. The kinetic energy T is:
T =
m 2
RΩ cos (tΩ) + X˙(t)(sin ϕ(t)) + cos (ϕ(t)) · (L + X(t)) · ϕ˙(t)
m 2
RΩ sin (tΩ) − X˙(t) cos ϕ(t) + sin (ϕ(t)) · (L + X(t)) · ϕ˙(t)
The potential energy V is:
V =
k 2
X(t) +
mg k
− mg (R cos (tΩ) + cos (ϕ(t)) · (L + X(t))) , (46)
The expressions for kinetic and potential in (45) and (46) have been adapted from Starosta et al., 2008 [7]. The Lagrangian is L = T − V , so the Euler–Lagrange equations are:
∂L ∂X
d dt
∂L
∂ X˙
∂L
∂ϕ
d dt
∂L
∂ ϕ ˙
Writing the first two equations of motion along with the equations found for X¨ and ϕ¨ (using a CAS) I get,
v = X˙ (49) ω = ˙ϕ (50)
v ˙ =
−mg + mRΩ^2 cos (tΩ − ϕ) + mg cos ϕ − kX + Lm ϕ˙^2 + mX ϕ˙^2 m
ω ˙ = mRΩ^2 sin (tΩ − ϕ) − mg sin ϕ − 2 m X˙ ϕ˙ m (L + X)
Using equations (49)–(52) for the RKf45 method I have the script below.
1 s u b r o u t i n e r k f 4 5 s t e p ( x , h , y ) 2!.. 3!.. 4!.. 5 c o n t a i n s 6 f u n c t i o n kv ( t , dt , y ) r e s u l t ( q )
4 2 0 2 4
x
3
2
1
0
1
2
y
Figure 9: Motion of the spring pendulum with initial velocity v = 10 m/s, t from 0 to 30 s. The origin is centered at the beginning of the expansion component X(t).
0 2 4 6 8 10 12 14 t
4
2
0
2
4
x ,^ y
(a) x, y vs t. x is orange.
1.5 1.0 0.5 0.0 0.5 1.0 1.
10
0
10
20
30
40
(b) Angular phase space.
0 2 4 6 8 10 12 14 t
(c) ϕ vs t.
Figure 10: Various plots relevant to t varying between 0 and 15 seconds whilst v 0 = 10 m/s.
6 4 2 0 2 4 6
3
2
1
0
Figure 11: Motion of the spring pendulum with v 0 = 10 m/s, t from 0 to 15 s. The origin is centered at the spring pendulum suspension, i.e. wherever on the ring which the spring pendulum is attached to.
4 Discussion
After solving the Spring Pendulum, CR3BP, and Double Pendulum problems it became quite clear that the Runge–Kutta methods are powerful computational methods to find the motion and other relevant info of difficult and/or chaotic systems, given the equations of motion.
For the double pendulum, the RKF45 method I used gives a nearly perfect prediction for the behavior of the system. After normalizing the definition of it, the circular restricted three-body problem became easy to solve generally. Using the RK4 method, the motion of the comet in the Earth–Sun was plotted in section 2.3, its motion in the Jupiter–Sun system is plotted in Appendix A. As for the Spring Pendulum, the RKF method helped the find various information regarding the phenomenally nonlinear system.
References
[1] Zoltan Papp, Mastering Computational Physics, Phys 562 Lecture Notes, CSU Long Beach, 2020. [2] Erik Neumann, Double Pendulum (2019), myPhysicsLab.com, retrieved on 09/27/2020 from https://www. myphysicslab.com/pendulum/double-pendulum-en.html. [3] Gabriela Gonz´alez, Single and Double plane pendulum, Louisiana State University. [4] John Taylor, Classical Mechanics, chapter 11.4, University Science Books, 2005. [5] Lagrange Point (2020), Wikipedia, retrieved on 10/01/2020 from https://en.wikipedia.org/wiki/Lagrange_ point. [6] Nkosi N. Trim, Visualizing Solutions of the Circular Restricted Three-Body Problem, Masters Thesis, Rutgers University, 2008. [7] Starosta, Spyniewska-Kaminska, and Awrejcewicz, Parametric and External Resonances in Kinematically and Externally Excited Nonlinear Spring Pendulum, International Journal of Bifurcation and Chaos, Vol. 21, No. 10, World Scientific Publishing Company, 2011. [8] T.S. Amer & M.A. Bek, Chaotic Responses of a Harmonically Excited Spring Pendulum Moving in Circular Path, Nonlinear Analysis: Real World Applications, Elsevier Ltd., 2008.
Appendices
A CR3BP; Comet–Jupiter–Sun
(a) L 4 with v 0 = 0 and tmax = 300.
(b) L 1 with low v 0 and tmax = 200.
(c) L 4 with low v 0 and tmax = 200.
Figure A.1: Plots for the comet in the Lagrange points of the Jupiter–Sun system. All t 0 = 0.
