text:
## First steps with monads in Scala

### Introduction

In the previous post I gave a quick introduction to some important concepts in functional programming, such as HOFs, closures, currying and partial application, and hopefully gave some insight into why these concepts might be useful in the context of scientific computing. Another concept that is very important in modern functional programming is that of the monad. Monads are one of those concepts that turns out to be very simple and intuitive once you “get it”, but completely impenetrable until you do! Now, there zillions of monad tutorials out there, and I don’t think that I have anything particularly insightful to add to the discussion. That said, most of the tutorials focus on problems and examples that are some way removed from the interests of statisticians and scientific programmers. So in this post I want to try and give a very informal and intuitive introduction to the monad concept in a way that I hope will resonate with people from a more scientific computing background.

The term “monad” is borrowed from that of the corresponding concept in category theory. The connection between functional programming and category theory is strong and deep. I intend to expore this more in future posts, but for this post the connection is not important and no knowledge of category theory is assumed (or imparted!).

#### Maps and Functors

All of the code used in this post in contained in the first-monads directory of my blog repo. The best way to follow this post is to copy-and-paste commands one-at-a-time from this post to a Scala REPL or sbt console. Note that only the numerical linear algebra examples later in this post require any non-standard dependencies.

The map method is one of the first concepts one meets when beginning functional programming. It is a higher order method on many (immutable) collection and other container types. Let’s start by looking at how map operates on Lists.

val x = (0 to 4).toList
// x: List[Int] = List(0, 1, 2, 3, 4)
val x2 = x map { x => x * 3 }
// x2: List[Int] = List(0, 3, 6, 9, 12)
val x3 = x map { _ * 3 }
// x3: List[Int] = List(0, 3, 6, 9, 12)
val x4 = x map { _ * 0.1 }
// x4: List[Double] = List(0.0, 0.1, 0.2, 0.30000000000000004, 0.4)

The last example shows that a List[T] can be converted to a List[S] if map is passed a function of type T => S. Of course there’s nothing particularly special about List here. It works with other collection types in the same way, as the following example with (immutable) Vector illustrates:

val xv = x.toVector
// xv: Vector[Int] = Vector(0, 1, 2, 3, 4)
val xv2 = xv map { _ * 0.2 }
// xv2: scala.collection.immutable.Vector[Double] = Vector(0.0, 0.2, 0.4, 0.6000000000000001, 0.8)
val xv3 = for (xi <- xv) yield (xi * 0.2)
// xv3: scala.collection.immutable.Vector[Double] = Vector(0.0, 0.2, 0.4, 0.6000000000000001, 0.8)

Note here that the for comprehension generating xv3 is exactly equivalent to the map call generating xv2 – the for-comprehension is just syntactic sugar for the map call. The benefit of this syntax will become apparent in the more complex examples we consider later.

Many collection and other container types have a map method that behaves this way. Any parametrised type that does have a map method like this is known as a Functor. Again, the name is due to category theory, but that doesn’t matter for this post. From a Scala-programmer perspective, a functor can be thought of as a trait, in pseudo-code as

trait F[T] {
def map(f: T => S): F[S]
}

with F representing the functor. In fact it turns out to be better to think of a functor as a type class, but that is yet another topic for a future post… Also note that to be a functor in the strict sense (from a category theory perspective), the map method must behave sensibly – that is, it must satisfy the functor laws. But again, I’m keeping things informal and intuitive for this post – there are plenty of other monad tutorials which emphasise the category theory connections.

Once we can map functions over elements of containers, we soon start mapping functions which themselves return values of the container type. eg. we can map a function returning a List over the elements of a List, as illustrated below.

val x5 = x map { x => List(x - 0.1, x + 0.1) }
// x5: List[List[Double]] = List(List(-0.1, 0.1), List(0.9, 1.1), List(1.9, 2.1), List(2.9, 3.1), List(3.9, 4.1))

Clearly this returns a list-of-lists. Sometimes this is what we want, but very often we actually want to flatten down to a single list so that, for example, we can subsequently map over all of the elements of the base type with a single map. We could take the list-of-lists and then flatten it, but this pattern is so common that the act of mapping and then flattening is often considered to be a basic operation, often known in Scala as flatMap. So for our toy example, we could carry out the flatMap as follows.

val x6 = x flatMap { x => List(x - 0.1, x + 0.1) }
// x6: List[Double] = List(-0.1, 0.1, 0.9, 1.1, 1.9, 2.1, 2.9, 3.1, 3.9, 4.1)

The ubiquity of this pattern becomes more apparent when we start thinking about iterating over multiple collections. For example, suppose now that we have two lists, x and y, and that we want to iterate over all pairs of elements consisting of one element from each list.

val y = (0 to 12 by 2).toList
// y: List[Int] = List(0, 2, 4, 6, 8, 10, 12)
val xy = x flatMap { xi => y map { yi => xi * yi } }
// xy: List[Int] = List(0, 0, 0, 0, 0, 0, 0, 0, 2, 4, 6, 8, 10, 12, 0, 4, 8, 12, 16, 20, 24, 0, 6, 12, 18, 24, 30, 36, 0, 8, 16, 24, 32, 40, 48)

This pattern of having one or more nested flatMaps followed by a final map in order to iterate over multiple collections is very common. It is exactly this pattern that the for-comprehension is syntactic sugar for. So we can re-write the above using a for-comprehension

val xy2 = for {
xi <- x
yi <- y
} yield (xi * yi)
// xy2: List[Int] = List(0, 0, 0, 0, 0, 0, 0, 0, 2, 4, 6, 8, 10, 12, 0, 4, 8, 12, 16, 20, 24, 0, 6, 12, 18, 24, 30, 36, 0, 8, 16, 24, 32, 40, 48)

This for-comprehension (usually called a for-expression in Scala) has an intuitive syntax reminiscent of the kind of thing one might write in an imperative language. But it is important to remember that <- is not actually an imperative assignment. The for-comprehension really does expand to the pure-functional nested flatMap and map call given above.

Recalling that a functor is a parameterised type with a map method, we can now say that a monad is just a functor which also has a flatMap method. We can write this in pseudo-code as

trait M[T] {
def map(f: T => S): M[S]
def flatMap(f: T => M[S]): M[S]
}

So far the functors and monads we have been working with have been collections, but not all monads are collections, and in fact collections are in some ways atypical examples of monads. Many monads are containers or wrappers, so it will be useful to see examples of monads which are not collections.

One of the first monads that many people encounter is the Option monad (referred to as the Maybe monad in Haskell, and Optional in Java 8). You can think of it as being a strange kind of “collection” that can contain at most one element. So it will either contain an element or it won’t, and so can be used to wrap the result of a computation which might fail. If the computation succeeds, the value computed can be wrapped in the Option (using the type Some), and if it fails, it will not contain a value of the required type, but simply be the value None. It provides a referentially transparent and type-safe alternative to raising exceptions or returning NULL references. We can transform Options using map.

val three = Option(3)
// three: Option[Int] = Some(3)
val twelve = three map (_ * 4)
// twelve: Option[Int] = Some(12)

But when we start combining the results of multiple computations that could fail, we run into exactly the same issues as before.

val four = Option(4)
// four: Option[Int] = Some(4)
val twelveB = three map (i => four map (i * _))
// twelveB: Option[Option[Int]] = Some(Some(12))

Here we have ended up with an Option wrapped in another Option, which is not what we want. But we now know the solution, which is to replace the first map with flatMap, or better still, use a for-comprehension.

val twelveC = three flatMap (i => four map (i * _))
// twelveC: Option[Int] = Some(12)
val twelveD = for {
i <- three
j <- four
} yield (i * j)
// twelveD: Option[Int] = Some(12)

Again, the for-comprehension is a little bit easier to understand than the chaining of calls to flatMap and map. Note that in the for-comprehension we don’t worry about whether or not the Options actually contain values – we just concentrate on the “happy path”, where they both do, safe in the knowledge that the Option monad will take care of the failure cases for us. Two of the possible failure cases are illustrated below.

val oops: Option[Int] = None
// oops: Option[Int] = None
val oopsB = for {
i <- three
j <- oops
} yield (i * j)
// oopsB: Option[Int] = None
val oopsC = for {
i <- oops
j <- four
} yield (i * j)
// oopsC: Option[Int] = None

This is a typical benefit of code written in a monadic style. We chain together multiple computations thinking only about the canonical case and trusting the monad to take care of any additional computational context for us.

#### IEEE floating point and NaN

Those with a background in scientific computing are probably already familiar with the NaN value in IEEE floating point. In many regards, this value and the rules around its behaviour mean that Float and Double types in IEEE compliant languages behave as an Option monad with NaN as the None value. This is simply illustrated below.

val nan = Double.NaN
3.0 * 4.0
// res0: Double = 12.0
3.0 * nan
// res1: Double = NaN
nan * 4.0
// res2: Double = NaN

The NaN value arises naturally when computations fail. eg.

val nanB = 0.0 / 0.0
// nanB: Double = NaN

This Option-like behaviour of Float and Double means that it is quite rare to see examples of Option[Float] or Option[Double] in Scala code. But there are some disadvantages of the IEEE approach, as discussed elsewhere. Also note that this only works for Floats and Doubles, and not for other types, including, say, Int.

val nanC=0/0
// This raises a runtime exception!

#### Option for matrix computations

Good practical examples of scientific computations which can fail crop up frequently in numerical linear algebra, so it’s useful to see how Option can simplify code in that context. Note that the code in this section requires the Breeze library, so should be run from an sbt console using the sbt build file associated with this post.

In statistical applications, one often needs to compute the Cholesky factorisation of a square symmetric matrix. This operation is built into Breeze as the function cholesky. However the factorisation will fail if the matrix provided is not positive semi-definite, and in this case the cholesky function will throw a runtime exception. We will use the built in cholesky function to create our own function, safeChol (using a monad called Try which is closely related to the Option monad) returning an Option of a matrix rather than a matrix. This function will not throw an exception, but instead return None in the case of failure, as illustrated below.

import breeze.linalg._
def safeChol(m: DenseMatrix[Double]): Option[DenseMatrix[Double]] = scala.util.Try(cholesky(m)).toOption
val m = DenseMatrix((2.0, 1.0), (1.0, 3.0))
val c = safeChol(m)
// c: Option[breeze.linalg.DenseMatrix[Double]] =
// Some(1.4142135623730951  0.0
// 0.7071067811865475  1.5811388300841898  )

val m2 = DenseMatrix((1.0, 2.0), (2.0, 3.0))
val c2 = safeChol(m2)
// c2: Option[breeze.linalg.DenseMatrix[Double]] = None

A Cholesky factorisation is often followed by a forward or backward solve. This operation may also fail, independently of whether the Cholesky factorisation fails. There doesn’t seem to be a forward solve function directly exposed in the Breeze API, but we can easily define one, which I call dangerousForwardSolve, as it will throw an exception if it fails, just like the cholesky function. But just as for the cholesky function, we can wrap up the dangerous function into a safe one that returns an Option.

import com.github.fommil.netlib.BLAS.{getInstance => blas}
def dangerousForwardSolve(A: DenseMatrix[Double], y: DenseVector[Double]): DenseVector[Double] = {
val yc = y.copy
blas.dtrsv("L", "N", "N", A.cols, A.toArray, A.rows, yc.data, 1)
yc
}
def safeForwardSolve(A: DenseMatrix[Double], y: DenseVector[Double]): Option[DenseVector[Double]] = scala.util.Try(dangerousForwardSolve(A, y)).toOption

Now we can write a very simple function which chains these two operations together, as follows.

def safeStd(A: DenseMatrix[Double], y: DenseVector[Double]): Option[DenseVector[Double]] = for {
L <- safeChol(A)
z <- safeForwardSolve(L, y)
} yield z

safeStd(m,DenseVector(1.0,2.0))
// res15: Option[breeze.linalg.DenseVector[Double]] = Some(DenseVector(0.7071067811865475, 0.9486832980505138))

Note how clean and simple this function is, concentrating purely on the “happy path” where both operations succeed and letting the Option monad worry about the three different cases where at least one of the operations fails.

Let’s finish with a monad for parallel and asynchronous computation, the Future monad. The Future monad is used for wrapping up slow computations and dispatching them to another thread for completion. The call to Future returns immediately, allowing the calling thread to continue while the additional thread processes the slow work. So at that stage, the Future will not have completed, and will not contain a value, but at some (unpredictable) time in the future it (hopefully) will (hence the name). In the following code snippet I construct two Futures that will each take at least 10 seconds to complete. On the main thread I then use a for-comprehension to chain the two computations together. Again, this will return immediately returning another Future that at some point in the future will contain the result of the derived computation. Then, purely for illustration, I force the main thread to stop and wait for that third future (f3) to complete, printing the result to the console.

import scala.concurrent.duration._
import scala.concurrent.{Future,ExecutionContext,Await}
import ExecutionContext.Implicits.global
val f1=Future{
1 }
val f2=Future{
2 }
val f3=for {
v1 <- f1
v2 <- f2
} yield (v1+v2)
println(Await.result(f3,30.second))

When you paste this into your console you should observe that you get the result in 10 seconds, as f1 and f2 execute in parallel on separate threads. So the Future monad is one (of many) ways to get started with parallel and async programming in Scala.

## Parallel Monte Carlo using Scala

### Introduction

In previous posts I have discussed general issues regarding parallel MCMC and examined in detail parallel Monte Carlo on a multicore laptop. In those posts I used the C programming language in conjunction with the MPI parallel library in order to illustrate the concepts. In this post I want to take the example from the second post and re-examine it using the Scala programming language.

The toy problem considered in the parallel Monte Carlo post used $10^9$ $U(0,1)$ random quantities to construct a Monte Carlo estimate of the integral

$\displaystyle I=\int_0^1\exp\{-u^2\}du$.

A very simple serial program to implement this algorithm is given below:

import java.util.concurrent.ThreadLocalRandom
import scala.math.exp
import scala.annotation.tailrec

object MonteCarlo {

@tailrec
def sum(its: Long,acc: Double): Double = {
if (its==0)
(acc)
else {
sum(its-1,acc+exp(-u*u))
}
}

def main(args: Array[String]) = {
println("Hello")
val iters=1000000000
val result=sum(iters,0.0)
println(result/iters)
println("Goodbye")
}

}

Note that ThreadLocalRandom is a parallel random number generator introduced into recent versions of the Java programming language, which can be easily utilised from Scala code. Assuming that Scala is installed, this can be compiled and run with commands like

scalac monte-carlo.scala
time scala MonteCarlo

This program works, and the timings (in seconds) for three runs are 57.79, 57.77 and 57.55 on the same laptop considered in the previous post. The first thing to note is that this Scala code is actually slightly faster than the corresponding C+MPI code in the single processor special case! Now that we have a good working implementation we should think how to parallelise it…

### Parallel implementation

Before constructing a parallel implementation, we will first construct a slightly re-factored serial version that will be easier to parallelise. The simplest way to introduce parallelisation into Scala code is to parallelise a map over a collection. We therefore need a collection and a map to apply to it. Here we will just divide our $10^9$ iterations into $N=4$ separate computations, and use a map to compute the required Monte Carlo sums.

import java.util.concurrent.ThreadLocalRandom
import scala.math.exp
import scala.annotation.tailrec

object MonteCarlo {

@tailrec
def sum(its: Long,acc: Double): Double = {
if (its==0)
(acc)
else {
sum(its-1,acc+exp(-u*u))
}
}

def main(args: Array[String]) = {
println("Hello")
val N=4
val iters=1000000000
val its=iters/N
val sums=(1 to N).toList map {x => sum(its,0.0)}
val result=sums.reduce(_+_)
println(result/iters)
println("Goodbye")
}

}

Running this new code confirms that it works and gives similar estimates for the Monte Carlo integral as the previous version. The timings for 3 runs on my laptop were 57.57, 57.67 and 57.80, similar to the previous version of the code. So far so good. But how do we make it parallel? Like this:

import java.util.concurrent.ThreadLocalRandom
import scala.math.exp
import scala.annotation.tailrec

object MonteCarlo {

@tailrec
def sum(its: Long,acc: Double): Double = {
if (its==0)
(acc)
else {
sum(its-1,acc+exp(-u*u))
}
}

def main(args: Array[String]) = {
println("Hello")
val N=4
val iters=1000000000
val its=iters/N
val sums=(1 to N).toList.par map {x => sum(its,0.0)}
val result=sums.reduce(_+_)
println(result/iters)
println("Goodbye")
}

}

That’s it! It’s now parallel. Studying the above code reveals that the only difference from the previous version is the introduction of the 4 characters .par in line 22 of the code. R programmers will find this very much analagous to using lapply() versus mclapply() in R code. The function par converts the collection (here an immutable List) to a parallel collection (here an immutable parallel List), and then subsequent maps, filters, etc., can be computed in parallel on appropriate multicore architectures. Timings for 3 runs on my laptop were 20.74, 20.82 and 20.88. Note that these timings are faster than the timings for N=4 processors for the corresponding C+MPI code…

### Varying the size of the parallel collection

We can trivially modify the previous code to make the size of the parallel collection, N, a command line argument:

import java.util.concurrent.ThreadLocalRandom
import scala.math.exp
import scala.annotation.tailrec

object MonteCarlo {

@tailrec
def sum(its: Long,acc: Double): Double = {
if (its==0)
(acc)
else {
sum(its-1,acc+exp(-u*u))
}
}

def main(args: Array[String]) = {
println("Hello")
val N=args(0).toInt
val iters=1000000000
val its=iters/N
val sums=(1 to N).toList.par map {x => sum(its,0.0)}
val result=sums.reduce(_+_)
println(result/iters)
println("Goodbye")
}

}

We can now run this code with varying sizes of N in order to see how the runtime of the code changes as the size of the parallel collection increases. Timings on my laptop are summarised in the table below.

 N     T1     T2     T3
1   57.67  57.62  57.83
2   32.20  33.24  32.76
3   26.63  26.60  26.63
4   20.99  20.92  20.75
5   20.13  18.70  18.76
6   16.57  16.52  16.59
7   15.72  14.92  15.27
8   13.56  13.51  13.32
9   18.30  18.13  18.12
10   17.25  17.33  17.22
11   17.04  16.99  17.09
12   15.95  15.85  15.91

16   16.62  16.68  16.74
32   15.41  15.54  15.42
64   15.03  15.03  15.28

So we see that the timings decrease steadily until the size of the parallel collection hits 8 (the number of processors my hyper-threaded quad-core presents via Linux), and then increases very slightly, but not much as the size of the collection increases. This is better than the case of C+MPI where performance degrades noticeably if too many processes are requested. Here, the Scala compiler and JVM runtime manage an appropriate number of threads for the collection irrespective of the actual size of the collection. Also note that all of the timings are faster than the corresponding C+MPI code discussed in the previous post.

However, the notion that the size of the collection is irrelevant is only true up to a point. Probably the most natural way to code this algorithm would be as:

import java.util.concurrent.ThreadLocalRandom
import scala.math.exp

object MonteCarlo {

def main(args: Array[String]) = {
println("Hello")
val iters=1000000000
val sums=(1 to iters).toList map {x => ThreadLocalRandom.current().nextDouble()} map {x => exp(-x*x)}
val result=sums.reduce(_+_)
println(result/iters)
println("Goodbye")
}

}

or as the parallel equivalent

import java.util.concurrent.ThreadLocalRandom
import scala.math.exp

object MonteCarlo {

def main(args: Array[String]) = {
println("Hello")
val iters=1000000000
val sums=(1 to iters).toList.par map {x => ThreadLocalRandom.current().nextDouble()} map {x => exp(-x*x)}
val result=sums.reduce(_+_)
println(result/iters)
println("Goodbye")
}

}

Although these algorithms are in many ways cleaner and more natural, they will bomb out with a lack of heap space unless you have a huge amount of RAM, as they rely on having all $10^9$ realisations in RAM simultaneously. The lesson here is that even though functional languages make it very easy to write clean, efficient parallel code, we must still be careful not to fill up the heap with gigantic (immutable) data structures…

## Introduction

R is different to many “easy to use” statistical software packages – it expects to be given commands at the R command prompt. This can be intimidating for new users, but is at the heart of its power. Most powerful software tools have an underlying scripting language. This is because scriptable tools are typically more flexible, and easier to automate, script, program, etc. In fact, even software packages like Excel or Minitab have a macro programming language behind the scenes available for “power users” to exploit.

### Programming from the ground up

It is natural to want to automate (repetitive) tasks on a computer, to automate a “work flow”. This is especially natural for computational tasks, as all software tools are built from programming language components, anyway. In R, you do stuff by executing a sequence of commands. By putting a bunch of commands one after another into a text file, we can source the file, and script R. Scripting is the simplest form of programming – automating a sequence of tasks. Indeed, in Unix (including Linux and MacOS), we can put a bunch of Unix shell commands together in a shell script. In Windows, you can put a bunch of terminal commands together in a batch file.

Next, one can add in simple control structures, to support looping, branching and conditional execution. Looping allows repetition of very similar tasks. Branching and conditional execution allow decisions to be made depending on what has already happened. Most scripting languages support simple control structures – this allows carrying out of tasks which we could do in principle, but perhaps not in practice, due to the laborious and repetitive nature of some work-flows. We can go a long way with this, but…

Although scripting is a simple form of programming, it isn’t “real” programming, or software engineering. Software engineering is about developing flexible, modular, robust, re-usable, generic program components, and using them to build large, complex software systems – modularity is absolutely key here. Functions and procedures are a first step towards introducing modularity, allowing the development of “real” software. Proper support for these tends to distinguish “real” programming languages from scripting languages (though many modern “scripting” languages have at least some limited support, and the distinction between scripting languages and “real” languages is now very blurred).

## Functions and procedures

Procedures (or subroutines) are re-usable pieces of code which can be called from other pieces of code when needed. They may be provided with inputs, but do not have to be. They are usually called for their “side-effects”, such as doing plots, changing global variables, or reading/writing data to/from disk.

Functions are also re-usable pieces of code, but are mainly used to obtain a return-value that is computed on the basis of the given inputs. “Pure” functions do not have any side-effects. Functions and procedures may be combined in a hierarchical way to build large, complex algorithms from much simpler modular components. Note that many languages (including R), do not make a distinction between functions and procedures in the syntax of the language, but conceptually the distinction is really quite important.

## Variable scope

Almost all programming languages allow the definition of variables which are labels or tags representing or pointing at some value that may be defined and re-defined at run-time. In most modern programming languages, functions can define local variables which can be used in addition to any inputs (formal parameters) of the function – these are very important for the development of modular, re-usable code components. In particular, they help to avoid unanticipated name clashes in the global name-space. If a function refers to a variable which is neither a formal parameter nor a local variable, then a rule is needed to find which (if any) variable with that label is in scope for the function, so that the program can know what value to use.

### Dynamic scope

Under dynamic scope, if an “unknown” variable is referred to in a function, the idea is to use the version of the variable that is in scope at the time that the function was called (and apply this rule recursively) – this is the scoping rule used by the S-PLUS implementation of the S language. Dynamic scope was common among early dynamic programming languages – including early implementations of LISP (and is still used in Emacs LISP), as it was quite intuitive and natural to implement using a stack-based approach similar to the stack-based approach to passing variables in and out of subroutines commonly used by machine code and assembly programmers.

Despite being intuitively appealing, at least initially, there are a number of problems with dynamic scope in practice. In particular, we can’t really know by code inspection whether or not a given section of code will run in all situations without actually running the code, as we can’t know whether all variable bindings will resolve correctly. This is an issue even for dynamic languages, but is particularly problematic for strongly typed compiled languages, as it becomes difficult for the compiler to figure out the types of all variables correctly and therefore generate the appropriate byte-code. It is also very difficult for a function to have associated state – to do this, you must somehow get state variables into global name-space where they then become vulnerable to masking and name clashes. See the Wikipedia page on scope for further details.

### Lexical scope

Under lexical scoping rules, if an “unknown” variable is referred to in a function, the idea is to use the version that is “in scope” in the enclosing piece of code (and apply this rule recursively) — this is the scoping rule used by R (as R is built on top of a Scheme interpreter, a LISP derivative which emphasises lexical scope). Variable bindings can be all resolved, checked and verified at compile-time – this is safer, and in many other ways better. Most modern languages adopt lexical scoping, including most functional languages, such as LISPs (including LISP-STAT) and derivatives. In fact, I first read about lexical scope, function closures and their use in statistical computing in Luke Tierney’s LISP-STAT book (Tierney, 1990) in the early 1990s. That book was published over 20 years ago, so it just goes to show that there is nothing new about these functional programming approaches. In fact, although Tierney’s book describes a now obsolete system, I would nevertheless recommend reading it if you can find a copy, as I think it is still one of the best books on statistical computing ever written. It really puts the recent glut of horrible R-themed books to shame!

Given that R has been lexically scoped and has supported function closures since day one, it is reasonable to wonder why this programming style is not used more widely in R code. I think it is the difference in scoping rules between S-PLUS and R that has led to a fear of developing any R code which relies on non-local scoping rules. Certainly, in the early days of R, I would use S-PLUS at work and R at home, and I would want my code to work in exactly the same in both places! This is a shame, as lexical scoping is very powerful, and exploited widely in functional programming styles. The use of lexical scope and function closures in R is described quite nicely in Gentleman (2008), along with many other things.

To make sure that the concepts are clear, inspect the following piece of code and figure out what the result of the final function call will be. The answer is given below the code, so try not to peek before reading on…

a=1
b=2
f<-function(x)
{
a*x + b
}
g<-function(x)
{
a=2
b=1
f(x)
}
g(2)

No, really, try and figure it out before reading on for the answer! Understanding this example is key to understanding the difference between lexical and dynamic scope. Clearly the obvious answers are 4 and 5. If you didn’t get one of those, go back and try again! 😉 So, one of those is the result you get in a dynamically scoped language like S-PLUS, and the other is the result that you get in a lexically scoped language like R. But which is which? Many people when asked what this code does give the answer 5. This is the result for a dynamically scoped language. It is not the answer you get in R. In R, you get the answer 4. This is because f() was defined in the global environment, so it is the global bindings of a and b which count. Although the function g() defines its own local bindings for a and b, these have no impact on the global bindings, and are simply not relevant to the evaluation of f().

## Function closures

Some languages (including LISPs and derivatives such as Scheme, Python, and R) have functions as “first class objects”, which means that a function is able to return as its value another function. If the function (fChild) returned by the function (fParent) refers to variables not local to fChild, then scoping rules must apply to the resolution of the variable binding. If the language is lexically scoped, then the binding is determined by the variables in scope within the function fParent. The function fChild therefore has an associated environment, which provides bindings for non-local variable references – this allows maintaining of state. A function together with its environment is referred to as a function closure, and is a very powerful programming tool. Below is some more code to help illustrate what is going on. Again, try to figure out the result of the final function call before reading on for the answer and explanation…

a=1
b=2
f<-function(a,b)
{
return( function(x) {
a*x + b
})
}
g=f(2,1)
g(2)

Here, the function g(), together with its associated environment, is referred to as a function closure. See the Wikipedia page for closure for further details. So, what is the result of calling g(2) in this case? Again, some people get this wrong, and give the answer 4. This isn’t what you get in R – in R you get 5, again due to lexical scope. The point is that the function g() is created inside f(), and so it is the variable bindings in scope within f() at the time g() was created which matter. Since f() has a and b as formal arguments, these mask the global variables of the same name, so it is the 2 and 1 that are passed into f() to create g() which matter in the evaluation of g(). This is why function closures are so powerful. They are not simply functions, they are functions together with an associated environment, and the associated environment allows function closures to have associated state. Here the state corresponds to the values of a and b that were used in the creation of g(), but in principle the state can be essentially any data structure.

## Function closures for scientific computing

Function closures have numerous important applications in a variety of problems in scientific computing that involve dealing in some way with the “function environment problem”. There is quite a nice discussion of this issue in Oliveira and Stewart (2006), in the context of several strongly typed compiled languages. Consider, for example, a function that will numerically integrate a univariate function using (say) the trapezium rule. This integration function might expect that you pass in the function to be integrated, together with the limits of integration, and possibly a step size. Most likely this integration function will expect that the function passed in is univariate. However, in practice many functions have additional parameters (eg. the straight line example, above, which was a function of x, but depending on additional parameters a and b). This problem is solved by passing in a univariate function closure that contains the necessary environment to evaluate this univariate function correctly. Similar considerations apply for functions that carry out optimisation, solve ODEs by passing in the RHS, etc.

### The smfsb R package

The second edition of my textbook, Stochastic Modelling for Systems Biology, has recently been published (Wilkinson, 2011). The second edition has an associated R package, smfsb, available from CRAN – I gave a tutorial introduction in a previous post. The code makes extensive use of lexical scope and function closures, precisely to solve the function environment problem…

## References

• Oliveira, S, Stewart, D.E. (2006) Writing scientific software, CUP.
• Gentleman, R. (2008) R Programming for Bioinformatics, Chapman & Hall/CRC Press.
• Tierney, L. (1990) LISP-STAT, Wiley.
• Wilkinson (2011), Stochastic Modelling for Systems Biology, second edition, Chapman & Hall/CRC Press.
• ## Java math libraries and Monte Carlo simulation codes

### Java libraries for (non-uniform) random number simulation

Anyone writing serious Monte Carlo (and MCMC) codes relies on having a very good and fast (uniform) random number generator and associated functions for generation of non-uniform random quantities, such as Gaussian, Poisson, Gamma, etc. In a previous post I showed how to write a simple Gibbs sampler in four different languages. In C (and C++) random number generation is easy for most scientists, as the (excellent) GNU Scientific Library (GSL) provides exactly what most people need. But it wasn’t always that way… I remember the days before the GSL, when it was necessary to hunt around on the net for bits of C code to implement different algorithms. Worse, it was often necessary to hunt around for a bit of free FORTRAN code, and compile that with an F77 compiler and figure out how to call it from C. Even in the early Alpha days of the GSL, coverage was patchy, and the API changed often. Bad old days… But those days are long gone, and C programmers no longer have to worry about the problem of random variate generation – they can safely concentrate on developing their interesting new algorithm, and leave the rest to the GSL. Unfortunately for Java programmers, there isn’t yet anything quite comparable to the GSL in Java world.

I pretty much ignored Java until Java 5. Before then, the language was too limited, and the compilers and JVMs were too primitive to really take seriously for numerical work. But since the launch of Java 5 I’ve been starting to pay more interest. The language is now a perfectly reasonable O-O language, and the compilers and JVMs are pretty good. On a lot of benchmarks, Java is really quite comparable to C/C++, and Java is nicer to code, and has a lot of impressive associated technology. So if there was a math library comparable to the GSL, I’d be quite tempted to jump ship to the Java world and start writing all of my Monte Carlo codes in Java. But there isn’t. At least not yet.

When I first started to take Java seriously, the only good math library with good support for non-uniform random number generation was COLT. COLT was, and still is, pretty good. The code is generally well-written, and fast, and the documentation for it is reasonable. However, the structure of the library is very idiosyncratic, the coverage is a bit patchy, and there doesn’t ever seem to have been a proper development community behind it. It seems very much to have been a one-man project, which has long since stagnated. Unsurprisingly then, COLT has been forked. There is now a Parallel COLT project. This project is continuing the development of COLT, adding new features that were missing from COLT, and, as the name suggests, adding concurrency support. Parallel COLT is also good, and is the main library I currently use for random number generation in Java. However, it has obviously inherited all of the idiosyncrasies that COLT had, and still doesn’t seem to have a large and active development community associated with it. There is no doubt that it is an incredibly useful software library, but it still doesn’t really compare to the GSL.

I have watched the emergence of the Apache Commons Math project with great interest (not to be confused with Uncommons Math – another one-man project). I think this project probably has the greatest potential for providing the Java community with their own GSL equivalent. The Commons project has a lot of momentum, the Commons Math project seems to have an active development community, and the structure of the library is more intuitive than that of (Parallel) COLT. However, it is early days, and the library still has patchy coverage and is a bit rough around the edges. It reminds me a lot of the GSL back in its Alpha days. I’d not bothered to even download it until recently, as the random number generation component didn’t include the generation of gamma random quantities – an absolutely essential requirement for me. However, I noticed recently that the latest release (2.2) did include gamma generation, so I decided to download it and try it out. It works, but the generation of gamma random quantities is very slow (around 50 times slower than Parallel COLT). This isn’t a fundamental design flaw of the whole library – generating Gaussian random quantities is quite comparable with other libraries. It’s just that an inversion method has been used for gamma generation. All efficient gamma generators use a neat rejection scheme. In case anyone would like to investigate for themselves, here is a complete program for gamma generation designed to be linked against Parallel COLT:

import java.util.*;
import cern.jet.random.tdouble.*;
import cern.jet.random.tdouble.engine.*;

class GammaPC
{

public static void main(String[] arg)
{
DoubleRandomEngine rngEngine=new DoubleMersenneTwister();
Gamma rngG=new Gamma(1.0,1.0,rngEngine);
long N=10000;
double x=0.0;
for (int i=0;i<N;i++) {
for (int j=0;j<1000;j++) {
x=rngG.nextDouble(3.0,25.0);
}
System.out.println(x);
}
}

}

and here is a complete program designed to be linked against Commons Math:

import java.util.*;
import org.apache.commons.math.*;
import org.apache.commons.math.random.*;

class GammaACM
{

public static void main(String[] arg) throws MathException
{
RandomDataImpl rng=new RandomDataImpl();
long N=10000;
double x=0.0;
for (int i=0;i<N;i++) {
for (int j=0;j<1000;j++) {
x=rng.nextGamma(3.0,1.0/25.0);
}
System.out.println(x);
}
}

}

The two codes do the same thing (note that they parameterise the gamma distribution differently). Both programs work (they generate variates from the same, correct, distribution), and the Commons Math interface is slightly nicer, but the code is much slower to execute. I’m still optimistic that Commons Math will one day be Java’s GSL, but I’m not giving up on Parallel COLT (or C, for that matter!) just yet…


text:
# Math Help - Simplifying This Expression

1. ## Simplifying This Expression

I have the following:

${(4 - \frac{1}{(x+1)^{2}})^{1/2}}$

How do you go about simplifying this? DO you simplify each part individually by the 1/2 exponent? Thank you for any tips.

2. ## Re: Simplifying This Expression

you can't do it individually because it is the same as:
$\sqrt{4-\frac{1}{(x+1)^2}}$
$\sqrt{a-b}\neq\sqrt{a}-\sqrt{b}$
to do this you have to subtract $4-\frac{1}{(x-1)^2}$.can you do that?

3. ## Re: Simplifying This Expression

Originally Posted by anonimnystefy
you can't do it individually because it is the same as:
$\sqrt{4-\frac{1}{(x+1)^2}}$
$\sqrt{a-b}\neq\sqrt{a}-\sqrt{b}$
to do this you have to subtract $4-\frac{1}{(x-1)^2}$.can you do that?
To get $\frac{4(x-1)^{2} - 1}{(x-1)^{2}}$ ... ?

If that's correct, am I then able to take each term to the 1/2 exponent because it's all expressed with the same denominator?

4. ## Re: Simplifying This Expression

Originally Posted by YoungMarbleGiant
To get $\frac{4(x-1)^{2} - 1}{(x-1)^{2}}$ ... ?

If that's correct, am I then able to take each term to the 1/2 exponent because it's all expressed with the same denominator?
As you were told above, the square root of a sum/difference is NOT equal to the sum/difference of the square roots. So you can't take each term in the numerator to the 1/2 power.

However, the square root of a product/quotient DOES equal the product/quotient of the square roots. So you CAN take the square root of the numerator and the square root of the denominator.

5. ## Re: Simplifying This Expression

Originally Posted by Prove It
As you were told above, the square root of a sum/difference is NOT equal to the sum/difference of the square roots. So you can't take each term in the numerator to the 1/2 power.

However, the square root of a product/quotient DOES equal the product/quotient of the square roots. So you CAN take the square root of the numerator and the square root of the denominator.
So it simplifies to the square root of the whole numerator, over (x-1)... is that correct? Sorry if I've not followed you here.

6. ## Re: Simplifying This Expression

Originally Posted by YoungMarbleGiant
So it simplifies to the square root of the whole numerator, over (x-1)... is that correct? Sorry if I've not followed you here.
The denominator should be |x - 1|, otherwise you're correct.

format: plain


text:
CAS：78-93-3

CAS：108-94-1

CAS：67-64-1

CAS：64-19-7

CAS：141-78-6

CAS：108-88-3

CAS：71-43-2

CAS：64-17-5

CAS：67-56-1

### Why is the enol form of acetoacetic ester less in percentage than the enol form of …

So, enolic form of aceto acetic ester is less stable but in case of acetyl acetone its more stable 9 7 Anik Chowdhury Read chemistry at college · Upvoted by Anubhab Halder, MSc Chemistry, Indian Institute of Technology, Kanpur (2023) and Łukasz Golon, 83

### Butanone - an overview ScienceDirect Topics

Different strains have been engineered for the production of butanone and 2-butanol (Table 12.2).Because meso-2,3-butanediol is the direct precursor of butanone and 2-butanol, a strain that can natively produce meso-2,3-butanediol is preferred. There are several industrial strains that can produce meso-2,3-butanediol with high tiers, such as K. pneumoniae, K. oxytoca, …

### The enolic [C5H10O]+· and [C3H6O]+· ions formed from aliphatic ketones - Beynon - 1974 - Organic Mass Spectrometry - Wiley Online Library

The structures of reactive and nonreactive ions formed by McLafferty rearrangement in a variety of ketones are studied using metastable ion abundances and the kinetic energy released in metastable ion fragmentations. The loss of methyl radicals is studied in [C 5 H 10 O] +. ions and the mechanism delineated for ions formed in an initial configuration corresponding both to …

### Bonds and lone pairs present in enolic form of butanone are explain

23/10/2018· Answer: There are $$12\sigma,9\pi \text{ and }2\text{ lone pairs}$$ present in enolic form of butanone.Explanation: Butanone is a chemical compound havin … ANKITRAI4561 ANKITRAI4561 23.10.2018 Chemistry Secondary School answered • …

### The enolic form of ethyl acetoacetate as below has

The correct option is A. 18 sigma bonds and 2 pi-bonds. The enolic form of ethyl acetoacetate has 16 single bonds i.e. 16 σ - bonds and 2 double bonds i.e. 2 σ - bonds and 2π - bonds. Hence, the given structure has 18 σ - bonds and 2π - bonds. Suggest Corrections.

### The enolic form of acetone contains

The enolic form of acetone contains a) 9 sigma bonds, 1 pi-bond and 2 lone pairs b) 8 sigma bonds, 2 pi-bonds and 2 lone pairs c) 10 sigma bonds, 1 pi-bond and 1 lone pair d) 9 sigma bonds, 2 pi-bonds and 1 lone pair Correct answer is option ''A''. …

### The enolic form of acetone contains: - Tardigrade

The enolic form of acetone contains: 1075 58 Chemical Bonding and Molecular Structure Report Error A 9 sigma, 1 pi bond and 2 lone pairs B 8 sigma, 2 pi bond and 2 lone pairs C 10 sigma, 1 pi bond and 1 lone pairs D 9 sigma, 2 pi bond and 1 lone pairs Solution: σ - bond = 9 π - bond = 1 L.P. = 2 KCET 2022 Chemistry Questions

### Ethyl acetoacetate ligand distribution in the course of titanium n-butoxide …

15/10/2014· The signals due to an enolic form of Eaa bonded to Ti were as follows: methyl carbons at 13.9 ppm (E1) and 24.2 ppm (E2), methylene carbon at 59.5 ppm (E4), methyne carbon at 87.9 ppm (E3) and quaternary carbons at 171.4 ppm (E5) and 184.1 ppm (E6).

### The enolic form of acetone contains(A) 9 sigma bonds, 1 pi bond, …

So the correct answer is option (A). We can cross check the nuer of lone pairs of electrons as well. There are 2 lone pairs present on the oxygen atom. Thus the correct answer is A). Note: We know that a pi bond can be formed only after the formation of a sigma bond. So in a double bond, we should count one bond as sigma and the other as pi.

### The enolic form of butanone contains: (a) \$$12 \\sigma \$$ bonds,

The enolic form of butanone contains:(a) $$12 \sigma$$ bonds, $$1 \pi$$ bond and 2 lone pairs of electrons(b) $$11 \sigma$$ bonds, $$1 \pi$$ bond and

### The enolic form of acetone contains:a) 9-sigma, 1-pi-bond and 2 …

Therefore, the answer is option (a) - the enolic form of acetone contains 9-sigma, 1-pi-bond and 2 lone pairs. Additional information: Acetone is also known as propanone. It is the smallest and simplest ketone. Note: Acetone is a colourless, volatile, flammable It

### The enolic form of acetone contains - Infinity Learn

Full Form Blog IIT-JEE NEET CBSE NCERT Scholarship General The enolic form of acetone contains A 9 sigma bonds, 1 pi bond and 2 lone pairs B 8 sigma bonds, 2 pi bond and 2 lone pairs C 10 sigma bonds, 1 pi bond and 1 lone pairs D 9 sigma bonds, 2 pi

### 5. The enolic form of acetone contains how many sigma, pi bonds …

5. The enolic form of acetone contains how many sigma, pi bonds and lone pairs?

### Butanone - Wikipedia

Butanone Butanone, also known as methyl ethyl ketone ( MEK ), [a] is an organic compound with the formula CH 3 C (O)CH 2 CH 3. This colourless liquid ketone has a sharp, sweet odor reminiscent of acetone. It is produced industrially on a large scale, but occurs in nature only in trace amounts. [7]

• The enolic form of butanone contains 12 σ bonds 1 π bon/cite>

The enolic form of butanone contains 12 σ bonds 1 π bond and 2 lone pairs of electrons 11 σ bonds 1 π bond and 2 lone pairs of electrons 12 σ bonds 2 π bond Please scroll down to see the correct answer and solution guide.

### The enolic form of ethyl acetoacetate as below has

The enolic form of ethyl acetoacetate as below has

### Butanone - Wikipedia

Butanone, also known as methyl ethyl ketone ( MEK ), [a] is an organic compound with the formula CH 3 C (O)CH 2 CH 3. This colourless liquid ketone has a sharp, sweet odor reminiscent of acetone. It is produced industrially on a large scale, but occurs in nature only in trace amounts. [7] It is partially soluble in water, and is commonly used

• Mexico''s Butanone Market Report 2022 - Prices, Size, Forecast, …/cite>

14/10/2021· Mexico - Butanone (Methyl Ethyl Ketone) - Market Analysis, Forecast, Size, Trends and Insights Update: COVID-19 Impact $499$2 390 Basic Edition (Single Report)

### The enolic form of ethyl acetoacetate as below has

Solution The correct option is A 18 sigma bonds and 2 pi-bonds The enolic form of ethyl acetoacetate has 16 single bonds i.e. 16 σ - bonds and 2 double bonds i.e. 2 σ - bonds and 2π - bonds. Hence, the given structure has 18 σ - bonds and 2π - bonds. Suggest Corrections 6 Similar questions Q. In C4− 3, the nuer of σ and π bonds are: Q.

### The enolic form of acetone contains

No. of σ bonds in enolic form : 3 + 1 + 1 + 1 + 1 + 2 = 9No. of π bonds in enolic form : 1No. of lone pairs of electrons in enolic form = 2 The enolic form of acetone containsa)9 sigma bonds, 1 pi-bond and 2 lone pairsb)8 sigma bonds, 2 pi-bonds and 2 lone pairsc)10 sigma bonds, 1 pi-bond and 1 lone paird)9 sigma bonds, 2 pi-bonds and 1 lone pairCorrect answer is option ''A''.

### Which is the more stable enol form? - Chemistry Stack Exchange

28/10/2014· The 2,4-pentanedione enol form shown is the more stable enol, since the lone pair on the enol oxygen can delocalize across the five atoms to the electronegative carbonyl oxygen. Hydrogen bonding is also a minor factor in its favor. The enol form of cyclopent-2,3-en-1-one is indeed (also) a conjugated system, but the lone pair on the enol oxygen

### The enolic form of acetone contains - Sarthaks eConnect

The enolic form of acetone contains A. 9 sigma bonds, 1 pi bond and 2 lone pairs B. 8 sigma bonds, 2 pi bonds and 2 lone pairs C. 10 sigma bonds, 1 pi bond and 1 lone pair D. 9 sigma bonds, 2 pi bonds and 1 lone pair. class-12 chemical-bonding Share It On 1

### The enolic form of butanone contains 11 CHEMICAL …

11/11/2021· The enolic form of butanone contains Class:11Subject: CHEMISTRYChapter: CHEMICAL BONDINGBook:R SHARMABoard:IIT JEEYou can ask any doubt from class 6-12, JEE, The enolic form of …

: Doubtnut: 7
• The enolic form of butanone contains 12 σ bonds 1 π bon/cite>

The enolic form of butanone contains 12 σ bonds 1 π bond and 2 lone pairs of electrons 11 σ bonds 1 π bond and 2 lone pairs of electrons 12 σ bonds 2 π bond Please scroll down to see the correct answer and solution guide.

### The enolic form of butanone contains - Sarthaks eConnect

The enolic form of butanone contains A. 12 sigma bonds, 1pi bond, and 2 lone pairs of electrons B. 11sigma bonds, 1pi bond, and 2 lone pairs of electrons C. 12sigma bonds, 1pi, and 1 lone pair of electrons D. 10sigma bonds, 2pi bond, and 2 lone pairs

### increasing order of enol content, CH3COCH2COCH3, CH3COCH3, CH3CHO…

Usually, the keto form is more stable than enol form. Enol form is stable when there is hydrogen bonding. It is stable when there is a carbon carbon double bond. The more substituted alkene makes the enol bond stable. Hence, CH3CHO

### The enolic form of acetone contains - Tardigrade

IIT JEE 1990: The enolic form of acetone contains (A) 9 sigma bonds, 1 pi bond and 2 lone pairs (B) 8 sigma bonds, 2 pi bonds and 2 lone pairs (C) 10

### SOLVED:The enolic form of acetone contains (a) 9 sigma bonds, 1 pi bond and 2 lone pair of electrons …

VIDEO ANSWER: when is the knowledge form of a certain contains what? How many Sigmund by? So first we have to draw what is the structure of the knowledge form of the asset? So in order to see Tony structured. Yeah.

### The enolic form of acetone contains - Tardigrade

Q. The enolic form of acetone contains 1291 64 AMU AMU 2012 Report Error A 9 sigma bonds, 2 pi bonds and 1 lone pair B 9 sigma bonds, 1 pi bond and 2 lone pairs C 8 sigma bonds, 2 pi bonds and 2 lone pairs D 10 sigma bonds, 1 pi bond and 1 lone pair Solution: (keto form)C H 3 −C ∣∣O − C H 3 ⇌ (enol form)C H 3 −C ∣∣O = C H 2

### The enolic [C5H10O]+· and [C3H6O]+· ions formed from aliphatic ketones - Beynon - 1974 - Organic Mass Spectrometry - Wiley Online Library

The structures of reactive and nonreactive ions formed by McLafferty rearrangement in a variety of ketones are studied using metastable ion abundances and the kinetic energy released in metastable ion fragmentations. The loss of methyl radicals is studied in [C 5 H 10 O] +. ions and the mechanism delineated for ions formed in an initial configuration corresponding both to …

Guan di miao, Yulong town

XingYang, Henan, China

7*24 Hours 365 Days


text:
# The Area Problem

The use of rectangles and triangles to calculate the areas of curved planar shapes goes back thousands of years.

For example, the area of a triangle is pretty easy to calculate: it's just one half the base times the height. But imagine taking a right triangle (red) and deforming the hypotenuse into a parabola (blue).

The blue region fits nicely inside of the red one, which allows us to come up with an upper bound on the area of the blue region. What upper bound on the area does this observation give us?

## The Area Problem

### Introduction to Calculus

# The Area Problem

In the last problem, we used a triangle to estimate a parabolic area, but it turns out that limits help us compute the area exactly!

In this unit, we'll explore in modern terms a technique similar to what the ancient mathematician Archimedes used to compute the parabolic area. We'll actually do this calculation much later in the course, but we won't leave you in suspense: the answer will turn out to be $$\frac13.$$

Instead of tackling this difficult problem right now, we'll learn the basic ideas underlying integral calculus, the second major pillar of our course, by solving an easier area problem. Specifically, we're going to compute the area of a circle from scratch.

This unit will rely quite a bit on basic geometry and trigonometry.

## The Area Problem

### Introduction to Calculus

# The Area Problem

Here's the bird's-eye view of the unit. We start with a regular 3-gon (or triangle) inscribed inside a large blue circle of radius $$r_{\text{outer}}.$$ This triangle has a smaller circle (red) inscribed within. The radius of this circle is $$r_{\text{inner}}.$$ (See the gif below if you have a hard time visualizing this!)

As we add more vertices to make a regular 4-gon, 5-gon, etc., the circles get closer and closer together and the $$n$$-gon looks more and more like a circle itself.

Over the next few problems, we'll find the area formula of the $$n$$-gon and use it to find the area formula for a circle in the limit as $$n$$ gets very, very large.

## The Area Problem

### Introduction to Calculus

# The Area Problem

Let's start with the triangle pictured below. From the options, select the correct value for the angle $$\theta.$$

## The Area Problem

### Introduction to Calculus

# The Area Problem

For the same figure, if the area of the solid black triangle is $$A_{\triangle},$$ find $\frac{A_{\triangle}}{r_{\text{outer}}^2}.$

Hint: Trigonometry tells us that the side opposite $$\theta$$ has length $$2 r_{\text{outer}} \sin \left( \frac{\theta}{2} \right),$$ and the height of the triangle formed by the solid black line opposite $$\theta$$ and the two blue dashed lines is $$r_{\text{outer}} \cos \left( \frac{\theta}{2} \right).$$

## The Area Problem

### Introduction to Calculus

# The Area Problem

Since the solid black triangle is “between” the red and blue circles, we know that $A_{\text{red circle}} \leq A_{\triangle} = \frac{3 \sqrt{3}}{4} r_{\text{outer}}^2 \leq A_{\text{blue circle}}.$

## The Area Problem

### Introduction to Calculus

# The Area Problem

For any $$n$$-gon, we'll have $A_{\text{red circle}} \leq A_{n} \leq A_{\text{blue circle}}.$

Furthermore, $$A_{\text{red circle}} \to A_{\text{blue circle}}$$ as $$n \to \infty.$$ Here, “$$\to$$” means “approaches.” So if we take $$n$$ very large, $$A_{n}$$ will give us the area formula for a circle!

## The Area Problem

### Introduction to Calculus

# The Area Problem

Now that we understand the triangle case, let's look at the figure for $$n = 5.$$

If the area of the solid black pentagon is $$A_{5},$$ find $\frac{A_{5}}{r_{\text{outer}}^2}.$

Hint: Use $$\sin(36^{\circ}) \approx 0.6$$ and $$\cos(36^{\circ}) \approx 0.8.$$

## The Area Problem

### Introduction to Calculus

# The Area Problem

Here's what we know so far: \begin{align} \frac{A_{\textcolor{red}{3}}}{r_{\text{outer}}^2} & = \textcolor{red}{3} \sin \left( \frac{180^{\circ}}{\textcolor{red}{3}} \right) \cos \left( \frac{180^{\circ}}{\textcolor{red}{3}} \right) \\ \frac{A_{\textcolor{red}{5}}}{r_{\text{outer}}^2} & = \textcolor{red}{5} \sin \left( \frac{180^{\circ}}{\textcolor{red}{5}} \right) \cos \left( \frac{180^{\circ}}{\textcolor{red}{5}} \right). \end{align} It looks like the formula for the area of the $$n$$-gon is shaping up to be $\frac{A_{\textcolor{red}{n}}}{r_{\text{outer}}^2} = \textcolor{red}{n} \sin \left( \frac{180^{\circ}}{\textcolor{red}{n}} \right) \cos \left( \frac{180^{\circ}}{\textcolor{red}{n}} \right).$ Is this correct?

## The Area Problem

### Introduction to Calculus

# The Area Problem

From the gif, we know that $$A_{\text{red circle}} \leq A_{n} \leq A_{\text{blue circle}}$$ and we expect that $$A_{\text{red circle}} \to A_{\text{blue circle}}$$ as $$r_{\text{outer}} \to r_{\text{inner}}$$ for $$n$$ very, very large.

The Python codex below calculates $$\frac{A_{3}}{r_{\text{outer}}^2} , \frac{A_{3}}{r_{\text{outer}}^2}, \dots , \frac{A_{N_{max}}}{r_{\text{outer}}^2}$$ given $$N_{max}$$. Remember, $\frac{A_{n}}{r_{\text{outer}}^2} = n \sin \left( \frac{180^{\circ}}{n} \right) \cos \left( \frac{180^{\circ}}{n} \right).$ Choose various (large!) values for $$N_{max}$$ and run the code.

What value do these numbers appear to approach?

import sys
import numpy

NMAX = 10;
n = 3;

while n <= NMAX:
n = n+1;
print(c)

## The Area Problem

### Introduction to Calculus

# The Area Problem

It's amazing to think that Eudoxus and Archimedes used similar methods long before Newton and Leibniz introduced the world to calculus.

It's tempting to put $$n = \infty$$ into the formula $$\frac{A_{n}}{r_{\text{outer}}^2} = n \sin \left( \frac{180^{\circ}}{n} \right) \cos \left( \frac{180^{\circ}}{n} \right)$$, but as we'll see, $$\infty$$ is not a number and we can't just plug it into a function as if it were. In other words, writing $\infty \sin \left( \frac{180^{\circ}}{\infty} \right) \cos \left( \frac{180^{\circ}}{\infty} \right) = \pi$ makes no sense at all.

However, we'll learn that a proper way of summarizing what we've uncovered here is to write $\lim\limits_{n \to \infty} n \sin\left(\frac{180^{\circ}}{n}\right) \cos \left( \frac{180^{\circ}}{n} \right) = \pi,$ which just means that the pattern in the numbers $$n \sin \left(\frac{180^{\circ}}{n}\right) \cos \left( \frac{180^{\circ}}{n} \right)$$ as $$n$$ tends to large values is $$\pi.$$

## The Area Problem

### Introduction to Calculus

# The Area Problem

We've uncovered the area formula for a circle $$A = \pi r^2.$$

Let's see if we can find the circumference of the circle in terms of the radius.

Using the same sort of reasoning and assuming that $$\lim\limits_{n \to \infty} n \sin\left(\frac{180^o}{n}\right) = \pi,$$ find the formula for the circumference of a circle in terms of its radius.

Note: You may wonder about what happened to the $$\cos\left(\frac{180^o}{n}\right)$$ term in the limit $$\lim\limits_{n \to \infty} n \sin\left(\frac{180^o}{n}\right) = \pi.$$ It turns out that this term approaches 1 as $$n \to \infty.$$ We'll prove it later in the course, but experimenting with some large values of $$n$$ should be pretty convincing at this stage.

## The Area Problem

### Introduction to Calculus

# The Area Problem

Each of the $$n$$-gons can be thought of as a collection of $$n$$ congruent triangles. We found the area of the circle by taking the sum of these triangular areas and then letting $$n$$ tend to $$\infty.$$

It's natural to wonder if a similar idea can find other areas, like the area bounded by $$y=x^2, y = 0, x = 0, x = 1$$ that we considered at the beginning of this unit.

With just a few rectangles we can find a pretty decent approximation to the area. With a few more, the approximation gets even better, just as adding more and more vertices to our $$n$$-gon gave us a shape closer and closer to a circle!

## The Area Problem

### Introduction to Calculus

# The Area Problem

This unit laid out the approach we'll take in later chapters to solve the area problem: use $$n$$ straight-edge figures to approximate the total area, then take the limit of this sum as $$n \rightarrow \infty.$$

The area problem opens up the vast vistas of integral calculus, which we'll explore in Chapters 6-9. It's surprising to realize how this simple notion allows us to compute a wide range of interesting quantities from physics, finance, chemistry, and engineering.

We've reached the end of the beginning. This short chapter gave us a glimpse of the fundamental ideas of calculus. Through calculus, we can compute rates of change, which give us the chance of finding extreme function values. Calculus also opens up the way to computing areas and understanding what an infinite sum truly means.

Behind all of these marvelous notions is the limit, to which we turn our attention in the next chapter.

## The Area Problem

### Introduction to Calculus

×


text:
MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

# Thread Subject: problem with cumtrapz

 Subject: problem with cumtrapz From: Sanaa Date: 7 Dec, 2012 17:51:08 Message: 1 of 4 Hi, I want to integrate the discrete function a*x_n*(1-x_n) from nr to t (take r=0.005) > > where a:0:0.001:4 and n=0,1,2,3,...., and t\in[nr,(n+1) r] I have tried the follwoing % define the vector of colors, to plot the dat aof reach value of r in % different color color_vec = ['b']; r=0.05; for rho = 0: 0.001:5 r_vec=r; % do a loop of values of r for j = 1 : size(r_vec,2) % define delay r = r_vec(j); % define coefficient rho %rho = 0.5; % define the number of deiscrete times in interval [k*r, (k+1)*r] n = 50; % define the nimber of iterations k = 1,.....,Nit Nit = 300; % define the initial vector which is a vector of size n x0 = 0.3*ones(1,n); % define vectors x_next and x_previous x_next = zeros(1,n); %x_(n+1) x_previous = zeros(1,n); %x_n time = zeros(1,n); % initialize x_previous x_previous = x0; time = zeros(1,n); %x_next = x_previous; time = linspace(0,r,n);          hold on for i = 1:Nit     time = linspace(i*r,(i+1)*r,Nit);     y_next=rho.*x_previous.*(1-x_previous);     x_next = x_previous+ cumtrapz(y_next);     x_previous = x_next;  end end       plot(rho,x_next(1),'.', 'Color', char(color_vec(j))) end When I execute the results, I found out that x and y go to starting from the second value for each? What's wrong in my code? Moreover, I want to integrate ''y'' from nr to time; how can I do this using cumtrapz? Thanks a lot.
 Subject: problem with cumtrapz From: someone Date: 7 Dec, 2012 18:51:08 Message: 2 of 4 "Sanaa" wrote in message ... > Hi, > I want to integrate the discrete function a*x_n*(1-x_n) from nr to t (take r=0.005) > > > where a:0:0.001:4 and n=0,1,2,3,...., and t\in[nr,(n+1) r] > I have tried the follwoing > % define the vector of colors, to plot the dat aof reach value of r in > % different color > color_vec = ['b']; > r=0.05; > for rho = 0: 0.001:5 > r_vec=r; > % do a loop of values of r > for j = 1 : size(r_vec,2) > % define delay > r = r_vec(j); > % define coefficient rho > %rho = 0.5; > % define the number of deiscrete times in interval [k*r, (k+1)*r] > n = 50; > % define the nimber of iterations k = 1,.....,Nit > Nit = 300; > % define the initial vector which is a vector of size n > x0 = 0.3*ones(1,n); > % define vectors x_next and x_previous > x_next = zeros(1,n); %x_(n+1) > x_previous = zeros(1,n); %x_n > > time = zeros(1,n); > > % initialize x_previous > x_previous = x0; > time = zeros(1,n); > %x_next = x_previous; > time = linspace(0,r,n); > > hold on > for i = 1:Nit > time = linspace(i*r,(i+1)*r,Nit); > y_next=rho.*x_previous.*(1-x_previous); > x_next = x_previous+ cumtrapz(y_next); > x_previous = x_next; > end > end > plot(rho,x_next(1),'.', 'Color', char(color_vec(j))) > end > > When I execute the results, I found out that x and y go to starting from the second value for each? What's wrong in my code? I'm not sure I understand the question and I don't have MATLAB installed on this computer, but I don't think your for loop is incrementing properly. For example, from the above, look at the section: r=0.05; for rho = 0: 0.001:5 r_vec=r; % do a loop of values of r for j = 1 : size(r_vec,2) % define delay r = r_vec(j); Since r & r_vec are scalars, the for loop will ALWAYS evaluate to: for j = 1 : 1 (I don't see where the values of r or r_vec are changed anywhere else in your code.) Perhaps you mean to use rho in place of one (or more) of the r or r_vec variables. Also, look at: > % define the initial vector which is a vector of size n > x0 = 0.3*ones(1,n); > % define vectors x_next and x_previous > x_next = zeros(1,n); %x_(n+1) > x_previous = zeros(1,n); %x_n > > time = zeros(1,n); > > % initialize x_previous > x_previous = x0; > time = zeros(1,n); > %x_next = x_previous; > time = linspace(0,r,n); Theres a lot of neeedless redundancy here. You define and overwrite time twice. Only the last statement is needed. Is that what you meant? Similar situation with x_previous. Perhaps its just extranious code leftover from debugging. You should cllean it up (or at least comment-out) the unneeded statements and resubmit your question if you still can't see an obvious problem. > Moreover, I want to integrate ''y'' from nr to time; how can I do this using cumtrapz? > > Thanks a lot.
 Subject: problem with cumtrapz From: Sanaa Date: 8 Dec, 2012 18:06:08 Message: 3 of 4 I am sorry for the miss forming of the code. Please see this % define the vector of colors, to plot the data of reach value of r in % different color color_vec = ['b']; r=0.05; for rho = 0: 0.001:10 % define the number of discrete times in interval [k*r, (k+1)*r] n = 50; % define the number of iterations k = 1,.....,Nit Nit = 150; % define the initial vector which is a vector of size n x0 = 0.3*ones(1,n); % define vectors x_next and x_previous x_next = zeros(1,n); %x_(n+1) x_previous = zeros(1,n); %x_n time = zeros(1,n); % initialize x_previous x_previous = x0; x_next = x_previous;     hold on for i = 1:Nit     time = linspace(i*r,(i+1)*r,Nit);     y_next=rho.*x_previous.*(1-x_previous);     x_next = x_previous+ cumtrapz(y_next);     x_previous = x_next;  end end       plot(rho,x_next(1),'.', 'Color', char(color_vec)) Now I have only one point in the plot!! I also want to integrate from n*r to t!? Any help please. thanks a lot in advance
 Subject: problem with cumtrapz From: Sanaa Date: 8 Dec, 2012 18:09:08 Message: 4 of 4 I am sorry for the miss forming of the code. Please see this % define the vector of colors, to plot the data of reach value of r in % different color color_vec = ['b']; r=0.05; for rho = 0: 0.001:10 % define the number of discrete times in interval [k*r, (k+1)*r] n = 50; % define the number of iterations k = 1,.....,Nit Nit = 150; % define the initial vector which is a vector of size n x0 = 0.3*ones(1,n); % define vectors x_next and x_previous x_next = zeros(1,n); %x_(n+1) x_previous = zeros(1,n); %x_n time = zeros(1,n); % initialize x_previous x_previous = x0; x_next = x_previous;     hold on for i = 1:Nit     time = linspace(i*r,(i+1)*r,Nit);     y_next=rho.*x_previous.*(1-x_previous);     x_next = x_previous+ cumtrapz(y_next);     x_previous = x_next;  end end       plot(rho,x_next(1),'.', 'Color', char(color_vec)) Now I have only one point in the plot!! I also want to integrate from n*r to t!? Any help please. thanks a lot in advance