More‎ > ‎

### 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:

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.

## Physical Fundamentals (optional)

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

1) 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.

2) Use of the Lagrange Multiplier Method (with a 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.

3) 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.

4) 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 soltutions 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.

## Differential-Algebraic Equations (optional)

To solve a damped pendulum we need some Differential-algebraic equations (DAEs) which describe the physical behaviour. These equations can then be solved with R or other programming languages.

One mathematical approach to get the solution is the Euler Lagrange Method with the Rayleigh Dissipation Function:

here qi 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 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. vx = dx/dt and v= 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:

 ``` 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``` ``` # 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) A 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) The 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 line 56:

`    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

 ``` 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``` ```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 = 0 ) ), 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 = 20, step = 0.1 #, 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')) ) ) ) ) ) ```

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.

### File server.R

 ``` 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 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209``` ```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") }) ```

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.