Repository navigation
Expand file tree
/
Copy pathsimple-pendulum-array.lua
More file actions
129 lines (112 loc) · 3.83 KB
/
Copy pathsimple-pendulum-array.lua
File metadata and controls
129 lines (112 loc) · 3.83 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
--
--------------------------------------------------------------------------------
-- File: simple-pendulum.lua
--
-- Usage: ./simple-pendulum.lua
--
-- Description: Solves the problem of the double pendulum using Runge-Kutta 4th order method
-- minth is 7.344e-12 if N = 1e8 and tN = 1.5 (took 114.26s).
-- Options: ---
-- Requirements: ---
-- Bugs: ---
-- Notes: ---
-- Author: Brenton Horne (), <brentonhorne77@gmail.com>
-- Organization:
-- Version: 1.0
-- Created: 16/10/17
-- Revision: ---
--------------------------------------------------------------------------------
-- Define parameters
local t0 = 0 -- Initial time
local tN = 4 -- Finishing time
local theta0 = 0 -- Angle from positive x axis (theta) at t0
local dtheta0 = 0 -- Change rate in theta at t0
local N = 1e7 -- Number of steps
local g = 9.8 -- Acceleration due to gravity
local l = 1 -- Length of pendulum
-- Initiate problem variables
local t = {} -- Initialize t array
local theta = {} -- Initialize theta array
local dtheta = {} -- Initialize dtheta array
t[1] = t0 -- Initiate t[1] variable
theta[1] = theta0 -- Initiate theta[1] variable
dtheta[1] = dtheta0 -- Initiate dtheta[1] variable
h = (tN - t0) / N -- define step size
minth = theta0 + math.pi -- Define minth at t0
-- Define the d2theta/dt2 = f(g, l, t, theta, dtheta) function
-- Essentially the RHS of the problem
function f(g, l, t, theta, dtheta)
K = - (g/l) * math.cos(theta)
return K
end
local k1, l1, k2, l2, k3, l3, k4, l4
-- Loop over time
for i = 1,N do
-- First approximation
k1 = h * f(g, l, t[i], theta[i], dtheta[i])
l1 = h * dtheta[i];
-- Second approximation
k2 = h * f(g, l, t[i] + h/2, theta[i] + 1/2 * l1, dtheta[i] + 1/2 * k1)
l2 = h * (dtheta[i] + 1/2 * k1)
-- Third approximation
k3 = h * f(g, l, t[i] + h/2, theta[i] + 1/2 * l2, dtheta[i] + 1/2 * k2)
l3 = h * (dtheta[i] + 1/2 * k2)
-- Fourth approximation
k4 = h * f(g, l, t[i] + h, theta[i] + l3, dtheta[i] + k3)
l4 = h * (dtheta[i] + 1/2 * k3)
-- Updating variables
t[i+1] = t[i] + h
dtheta[i+1] = dtheta[i] + (1 / 6) * (k1 + 2 * k2 + 2 * k3 + k4)
theta[i+1] = theta[i] + (1 / 6) * (l1 + 2 * l2 + 2 * l3 + l4)
-- Determining pi + theta; at minima it is 0
diff = math.abs(theta[i] + math.pi)
-- Checking if diff is smaller than the smallest minth so far and updating it if it is
if (diff < minth) then
minth = diff
end
end
print(minth)
--Plot theta against t
gp = require("gnuplot") -- Uses this library https://bitbucket.org/lucashnegri/lua-gnuplot/raw/34136a285a8820f31e1c63b7c81aa1b70a4b60ec/gnuplot.lua
local g = gp{
width = 1600,
height = 900,
xlabel = "X axis",
ylabel = "Y axis",
key = "top left",
consts = {
gamma = 2.5
},
data = {
gp.array {
{
t,
theta
},
},
}
}:plot('01-simple-pendulum-theta-by-t.png')
local g = gp{
width = 1600,
height = 900,
xlabel = "X axis",
ylabel = "Y axis",
key = "top left",
consts = {
gamma = 2.5
},
data = {
gp.array {
{
t,
dtheta
},
},
}
}:plot('02-simple-pendulum-dtheta-by-t.png')
--[[Plot theta against t
lp = require("lua-plot")
p = lp.plot{}
p:AddSeries(t, theta)
p:Show()
io.read()]]--