Simple Damped Pendulum - Solved and Visualised with R and Hosted as Shiny Application

BY MARKUS SPRUNCK

Although R is already more than 20 years old, in the last years this language has been becoming increasingly popular among data scientists for machine learning, statistical software, and general data analysis solutions. In Januar 2018 the language R even reached the top 10 of TIOBE index. With R you get large and active ecosystem for open source libraries and extensions.

This article describes how to solve the differential equations of a damped pendulum with the R language and visualise resulting diagrams in a web page. For this task the R library deSolve has been use to solve differential equations. A nice option to host a web page of your R application is shinyapps.io.

Expected Result

Just open the link https://markus-sprunck.shinyapps.io/pendulum/ and you should see something like the following screen shot:

https://sites.google.com/site/markussprunck/blog-1/damped-pendulum-with-language-r/Bildschirmfoto%202018-07-21%20um%2014.42.42.png

Figure 1: Screen shot of damped pendulum simulation with R

Change the parameter with the slider controls and after short time the graphs will be updated. All calculations and the rendering of the diagrams happen on server side. So, your implementation code is not exposed to the browser or user.

The eight diagrams depict the most important aspects of the physical behaviour of the damped pendulum. Notice that due to the damping the total energy is decreasing. The red dot indicates always the last state of the pendulum.

Before we get into details of the solution, some words on physical and mathematical background. In the case you are already familiar with the theory and/or just not interested, you may skip the next two chapters.

Fundamentals of Pendulum Physics

If you never have had lectures in advanced mechanics, you may like to brush up your knowledges with the following tutorials and descriptions:

Basics of Lagrangian Mechanics

Excelent video tutorials are the "Lectures in Lagrangian Mechanics" series on iLectureOnline, see here.

Start with part 1 for the basics and then you may proceed with part 6 and 7.

Use of the Lagrange Multiplier Method (friction less pendulum)

We like to calculate also the tension in the stab of the pendulum. One way to get the tension is the Lagrange Multiplier Method as described on page 7 of the lecture slides by Prof. N. Harnew; University of Oxford; HT 2017;

Classical Mechanics LECTURE 26: THE LAGRANGE EQUATION EXAMPLES; see here.

Introduce Damping with Rayleigh Dissipation Function

Now we have to describe the damping as velocity-proportional frictional forces. In combination with the Euler Lagrange Method we can use the Rayleigh Dissipation Function, see here.

Basics of Damped Mechanical Oscillations

When you try a new numerical method to solve a physical problem, it is alway a good idea to compare the solutions with an analytical model, some measurements and/or the results of a simplified model. You may compare the damping effect with Free Harmonic Oscillations, see here.

Solution

here q_i are generalised coordinates. The Lagrangian L for a system is defined as:

L = T - V

here T is the total kinetic energy of the system and V is the potential energy of the system. For the pendulum the Lagrange Function is:

here l is the length of the pendulum. The Lagrange multiplier method has been added to define the constraint (mass is moving along a circular path) of a stiff pendulum.

The Rayleigh Dissipation Function is defined in the following way:

The partial derivatives of the Lagrange Function with respect to the coordinates x, y and λ results the following three equations:

These two differential equations and the one algebraic equation describe the behavour of a simple damped pendulum.

Use of The R LIBRARY 'deSolve'

We describe the differential-algebraic equations as index 3 problem in Hessenberg form. This can be solved with the deSolve function radau, which uses Implicit Runge-Kutta.

First we introduce two new equations for speed, i.e. v_x= dx/dt and v_y= dy/dt ). Now we can do a substitution of the x deviation as a function of time.

The resulting five equations can be described in R in the following way:

# CALCULATED PARAMETER

polarCoordinateAnglele0 <- - pi / 2 + theta

x0 <- length * cos(polarCoordinateAnglele0)

y0 <- length * sin(polarCoordinateAnglele0)

vx0 <- speed * sin(polarCoordinateAnglele0)

vy0 <- speed * cos(polarCoordinateAnglele0)

lamda0 <- - mass * g * cos(theta) / length

yini <- c(x = x0, y = y0, vx = vx0, vy = vy0, lamda = lamda0)

times <- seq(from = 0, to = times_max, by = 0.001)

# SOLVE DIFFERENTIAL EQUATION

pendulum <- function (t, Y, parms) {

with (as.list(Y),

list(c( vx ,

vy ,

lamda * x - alpha * vx ,

lamda * y - g * mass - alpha * vy ,

x^2 + y^2 - length^2

))

)

}

M <- diag(nrow = 5)

M[1, 1] <- 1

M[2, 2] <- 1

M[3, 3] <- mass

M[4, 4] <- mass

M[5, 5] <- 0

index <- c(2, 2, 1)

out <- radau (y = yini, func = pendulum, parms = NULL,

time = times, mass = M, nind = index)

This description is not really intuitive, but with some experience with the R syntax it's getting more clear. You need to prepare three main things:

1) A Function,

that describes the right-hand side of the equation:

M dy/dt = f(t, y)

it is be defined as:

functionName <- function (t, y, parms) {

...

}

where t is the time point in the integration, y is the estimate of the variables, parms is a list of parameters (here not used). The return value of func should be a list, whose first element is a vector containing the derivatives of y with respect to time, and whose next elements are global values that are required at each point in times. The derivatives must be specified in the same order as the state variables y.

2) The Mass Matrix,

that describes the M on left-hand side of the equation:

M dy/dt = f(t, y)

the dimension is the number of y-values.

3) An Index,

is a three-valued vector with the number of variables of index 1, 2, 3 respectively. The equations have be defined such that the

index 1 variables precede the index 2 variables which in turn precede the index 3 variables. The sum of the variables of different index should equal N, the total number of variables. Here we have N = 5.

In the line you find the definition:

index <- c(2, 2, 1)

this means that we have the first 2 equations of index 1, the next 2 equations of index 2 and last equation of index 3.

You can find a more detailed description of the deSolve function radau, here.

Development Environment

The code has been developed in RStudo - this small IDE has good support for shinyapps.io.

Figure 2: RStudio is a free integrated development environment for R

Source Files

In the case you are not familiar with Shiny app you may start reading here.

Most Shiny application has at least two files, i.e.:

ui.R for user interface and controls and
server.R for calculation and rendering of diagrams.

File ui.R

The plots are directly inserted into the web page. If you uncomment line 66, it is possible to animate the time. This means that a small play button appears.

bootstrapPage(fluidPage(

titlePanel("DAE of Simple Damped Pendulum - Solved and Visualised with R"),

tags$br(),

sidebarLayout(

sidebarPanel(fluidRow(

column(

2,

sliderInput(

inputId = "mass",

label = "Mass [kg]:",

min = 10,

max = 100,

value = 30.0

)

),

column(

2,

sliderInput(

inputId = "length",

label = "Length [m]:",

min = 1,

max = 10,

value = 3.0

)

),

column(

2,

sliderInput(

inputId = "damping",

label = "Damping [Ns/m]:",

min = 0,

max = 16,

value = 4

)

),

column(

2,

sliderInput(

inputId = "speed",

label = "Speed [m/s]:",

min = -15.0,

max = 15.0,

value = 0.0,

step = 0.1

)

),

column(

2,

sliderInput(

inputId = "theta",

label = "Deflection [degree]:",

min = 1,

max = 90,

value = 15

)

),

column(

2,

sliderInput(

inputId = "max_time",

label = "Time [s]:",

min = 1,

max = 30,

value = 10 ,

step = 1.0

#, animate = TRUE

)

), width = 12),

fluidPage(

fluidRow(

column( 3, plotOutput('main_plot1')),

column( 3, plotOutput('main_plot2')),

column( 3, plotOutput('main_plot3')),

column( 3, plotOutput('main_plot4'))

),

fluidRow(

column( 3, plotOutput('main_plot5')),

column( 3, plotOutput('main_plot6')),

column( 3, plotOutput('main_plot7')),

column( 3, plotOutput('main_plot8'))

)

File server.R

The interesting thing here is, that the calculation happens just once and is reused for the rendering of each graph. This is possible with the reactive keyword.

library(deSolve)

library(minpack.lm)

# HELPER (converts cartesian coordinates to polar coordinate angle)

polarCoordinateAngle <- function(x, y) {

z <- x + 1i * y

res <- 90 - Arg(z) / pi * 180

res %% 360

}

# CONSTANTS

g <- 9.80665 # gravitational constant in m/s^2

function(input, output) {

result <- reactive({

# SET SIMULATION PARAMETER

length <- input$length # in m

speed <- input$speed # in m/s

mass <- input$mass # in kg

alpha <- input$damping # in Ns/m

theta <- input$theta / 180 * pi # in rad

times_max <- input$max_time # in seconds

# CALCULATED PARAMETER

polarCoordinateAnglele0 <- - pi / 2 + theta

x0 <- length * cos(polarCoordinateAnglele0)

y0 <- length * sin(polarCoordinateAnglele0)

vx0 <- speed * sin(polarCoordinateAnglele0)

vy0 <- speed * cos(polarCoordinateAnglele0)

lamda0 <- - mass * g * cos(theta) / length

yini <- c(x = x0, y = y0, vx = vx0, vy = vy0, lamda = lamda0)

times <- seq(from = 0, to = times_max, by = 0.001)

# SOLVE DIFFERENTIAL EQUATION

pendulum <- function (t, Y, parms) {

with (as.list(Y),

list(c( vx ,

vy ,

lamda * x - alpha * vx ,

lamda * y - g * mass - alpha * vy ,

x^2 + y^2 - length^2

))

)

}

M <- diag(nrow = 5)

M[1, 1] <- 1

M[2, 2] <- 1

M[3, 3] <- mass

M[4, 4] <- mass

M[5, 5] <- 0

index <- c(2, 2, 1)

out <- radau (y = yini, func = pendulum,

parms = NULL, time = times, mass = M, nind = index)

})

output$main_plot1 <- renderPlot({

# GET VALUES FOR RENDERING

out <- result()

# SET PLOT PARAMETER

par(mar=c(5, 5, 3, 2), cex.lab=1.6, cex.axis=1.6, cex.main=1.6, cex.sub=1.6)

# PLOT X-POSITION

p1 <- plot(out[, c("time", "x")], type = "l",

lwd = 1, ylab="X-Position [m]", xlab="Time [s]")

p1 <- p1 + grid(lwd = 1) + lines(out[, c("time", "x")], type = "l", lwd = 1)

p1 <- p1 + points(tail(out[, c("time", "x")], 1), type = "p", lwd = 3, col="red")

})

output$main_plot2 <- renderPlot({

# GET VALUES FOR RENDERING

out <- result()

# SET PLOT PARAMETER

par(mar=c(5, 5, 3, 2), cex.lab=1.6, cex.axis=1.6, cex.main=1.6, cex.sub=1.6)

# PLOT Y-POSITION

p2 <- plot(out[, c("time", "y")], type = "l",

lwd = 1, ylab="Y-Position [m]", xlab="Time [s]")

p2 <- p2 + grid(lwd = 1) + lines(out[, c("time", "y")], type = "l", lwd = 1)

p2 <- p2 + points(tail(out[, c("time", "y")], 1), type = "p", lwd = 3, col="red")

})

output$main_plot3 <- renderPlot({

# GET VALUES FOR RENDERING

out <- result()

xy <- out[, c("x", "y")]

# SET PLOT PARAMETER

par(mar=c(5, 5, 3, 2), cex.lab=1.6, cex.axis=1.6, cex.main=1.6, cex.sub=1.6)

# PLOT X-Y-DIAGRAM

p3 <- plot(xy, type = "l", lwd = 1, xlab="X-Position [m]", ylab="Y-Position [m]")

p3 <- p3 + grid(lwd = 1) + lines(out[, c("x", "y")], type = "l", lwd = 1)

p3 <- p3 + points(tail(xy, 1), type = "p", lwd = 3, col="red")

})

output$main_plot4 <- renderPlot({

# GET VALUES FOR RENDERING

out <- result()

time <- out[-1, c("time")]

xs <- out[-1, c("x")]

ys <- out[-1, c("y")]

theta <- 180 - polarCoordinateAngle(xs, ys)

# SET PLOT PARAMETER

par(mar=c(5, 5, 3, 2), cex.lab=1.6, cex.axis=1.6, cex.main=1.6, cex.sub=1.6)

# PLOT DEFLECTION

p4 <- plot(time, theta, type = "l", lwd = 0,

ylab="Deflection [degree]", xlab="Time [s]",

ylim=c(-max(theta)*1.1 ,max(theta)*1.1 ) )

p4 <- p4 + grid(lwd = 1) + lines(time, theta, type = "l", lwd = 1)

p4 <- p4 + points(tail(time, 1), tail(theta, 1), type = "p", lwd = 3, col="red")

})

output$main_plot5 <- renderPlot({

# GET VALUES FOR RENDERING

out <- result()

length <- input$length # in m

lamdas <- out[-1, c("lamda")]

tension <- - lamdas * length

times <- out[-1, c("time")]

# SET PLOT PARAMETER

par(mar=c(5, 5, 3, 2), cex.lab=1.6, cex.axis=1.6, cex.main=1.6, cex.sub=1.6)

# PLOT tension

p5 <- plot(times, tension, type = "l", lwd = 1,

ylab="Tension (N)", xlab="Time [s]")

p5 <- p5 + grid(lwd = 1) + lines(out[-1, c("time")], tension, type = "l")

p5 <- p5 + points(tail(times, 1),tail(tension, 1), type = "p", lwd = 3, col="red")

})

output$main_plot6 <- renderPlot({

# GET VALUES FOR RENDERING

out <- result()

mass <- input$mass # in kg

length <- input$length # in m

speed <- sqrt(out[-1, c("vx")] * out[-1, c("vx")]

+ out[-1, c("vy")] * out[-1, c("vy")])

E_k <- 0.5 * mass * speed * speed

E_p <- mass * g * (length + out[-1, c("y")])

E <- E_k + E_p

times <- out[-1, c("time")]

# SET PLOT PARAMETER

par(mar=c(5, 5, 3, 2), cex.lab=1.6, cex.axis=1.6, cex.main=1.6, cex.sub=1.6)

# PLOT ENERGY

p6 <- plot(times, E_p, type="l", lwd=0.5, col="grey",

ylab="Energy [J]", xlab="Time [s]", ylim=c(0,max(E)*1.20) )

p6 <- p6 + grid(lwd = 1)

p6 <- p6 + lines(times, E_k, type="l", lwd=0.5, col="blue" )

p6 <- p6 + lines(times, E, type="l", lwd=1.0, col="black")

p6 <- p6 + points(tail(times, 1), tail(E, 1), type = "p", lwd = 3, col="red")

})

output$main_plot7 <- renderPlot({

# GET VALUES FOR RENDERING

out <- result()

vxy <- out[-1, c("vx", "vy")]

# SET PLOT PARAMETER

par(mar=c(5, 5, 3, 2), cex.lab=1.6, cex.axis=1.6, cex.main=1.6, cex.sub=1.6)

# PLOT X-Y-DIAGRAM

p7 <- plot(vxy, type = "l", lwd = 1, xlab="X-Speed [m/s]", ylab="Y-Speed [m/s]")

p7 <- p7 + grid(lwd = 1) + lines(vxy, type = "l", lwd = 1)

p7 <- p7 + points(tail(vxy, 1), type = "p", lwd = 3, col="red")

})

output$main_plot8 <- renderPlot({

# GET VALUES FOR RENDERING

out <- result()

xs <- out[-1, c("x")]

ys <- out[-1, c("y")]

speed <- sqrt(out[-1, c("vx")] * out[-1, c("vx")]

+ out[-1, c("vy")] * out[-1, c("vy")])

theta <- 180 - polarCoordinateAngle(xs, ys)

# SET PLOT PARAMETER

par(mar=c(5, 5, 3, 2), cex.lab=1.6, cex.axis=1.6, cex.main=1.6, cex.sub=1.6)

# PLOT V-DEFLECTION-DIAGRAM

p8 <- plot(theta, speed, type = "l", lwd = 1,

xlab="Deflection [degree]", ylab="Speed [m/s]")

p8 <- p8 + grid(lwd = 1) + lines(theta, speed, type = "l", lwd = 1)

p8 <- p8 + points(tail(theta, 1), tail(speed, 1), type = "p", lwd = 3, col="red")

})

}