© 2008--2013 Krishna Myneni

- Basic Info
- Functions
- Complex Numbers
- Vectors, Matrices, and Higher Dimensional Arrays
- Plotting Data
- Matrix Computations
- R Packages

R is a free and powerful computing environment, designed for numerical computations and visualizing data. It is available for a wide variety of systems, and may be obtained from http://www.r-project.org. R is also a language for expressing, and programming numerical calculations. The language of R is object-oriented and elegant, with a simple and intuitive syntax. While many of the built-in functions of R are useful for performing calculations in statistics, the large variety of functions, as well as package extensions, allow R to be applied to a wide range of problems in the physical sciences and applied mathematics. These notes illustrate some common computational tasks encountered by students, teachers, and researchers in the physical sciences.

R version 2.6.2 (2008-02-08) Copyright (C) 2008 The R Foundation for Statistical Computing ISBN 3-900051-07-0 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. Natural language support but running in an English locale R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. >

The

> source("filename")where

> help("seq")If you are unsure of the exact function name, you may also search the documentation for a particular topic, using

> help.search("sequence")will show relevant functions for sequence generation.

> myfunc <- function(x) { y <- cos(x) + sin(x) return (y) } >

> myfunc(0.5) [1] 1.357008 >

> myfunc(0.5*1i) [1] 1.127626+0.521095i >

> myfunc(-2 + 5*1i) [1] -98.36115+36.59336i >

> x <- seq(-4, 4, by=0.01) > y <- myfunc(x)

> x = seq(-4, 4, by=0.01) > y = myfunc(x)

> plot(myfunc, 0, 2*pi)

> plot(myfunc2, 0, 2*pi, add=TRUE)

> plot(myfunc2, 0, 2*pi, add=TRUE, col="red")

We will take as an example, the problem of finding the roots
of the function,

f(x) = f_g(x) - f_l(x) where f_g(x) = exp(-a*x^2) and f_l(x) = b/(x^2 + b)For the example, we will take a=1, b=0.5.

The problem of finding the roots of *f(x)* is the same as the
problem of finding the values of *x* at which the two functions,
*f_g(x)* and *f_l(x)*, intersect. The problem can be solved
approximately by plotting the two functions, and visually determining
the values of *x* at which they intersect --- the graphical method.
A more precise method is to use an iterative search of $x$, to find
the values at which *f(x) = 0*. R provides such a search method with
the function **uniroot**. However, we must first obtain the approximate
roots with the graphical method, so that the search interval can be bounded
for **uniroot**.

> a = 1 > b = 0.5Define the functions,

> fg <- function(x) { return( exp(-a*x^2) ) } > fl <- function(x) { return( b/(x^2 + b) ) } > plot(fg, -5, 5, col='blue') > plot(fl, -5, 5, col='red', add=TRUE)The following graph will appear.

From the graph we see that there appear to be three points of intersection, one near

> f <- function(x) { return( fg(x) - fl(x) ) }Now, find the precise value of a given root (value of

> uniroot(f, c(-1.2, -0.8)) $root [1] -1.120907 $f.root [1] -4.398283e-08 $iter [1] 3 $estim.prec [1] 0.0001094043

Additional arguments can be given to the

> root1 = uniroot(f, c(-1.2, -0.8))$root > root1 [1] -1.120907

> z <- 3.2+4*1ior

> z = 3.2+4*1i

> Re(z) [1] 3.2 > Im(z) [1] 4 > Mod(z) [1] 5.122499 > abs(z) [1] 5.122499 > Arg(z) [1] 0.8960554The

> x = 3 + 1i*2 > Conj(x) [1] 3-2i >

> x <- c(1, 2, 3, 11, 12, 13) > x [1] 1 2 3 11 12 13 >

> mat1 <- matrix(c(1,2,3,11,12,13), nrow=2, byrow=TRUE) > mat1 [,1] [,2] [,3] [1,] 1 2 3 [2,] 11 12 13 >

> mat2 <- matrix(c(1+1i, 0, 0, 1-1i), nrow=2, byrow=TRUE) > mat2 [,1] [,2] [1,] 1+1i 0+0i [2,] 0+0i 1-1i >

> x <- seq(-4, 4, by=0.01) > length(x) [1] 801shows that we have created an 801 element vector using the

> mat1 <- matrix(c(1,2,3,11,12,13), nrow=2, byrow=TRUE) > dim(mat1) [1] 2 3

> x <- 1:10 > x [1] 1 2 3 4 5 6 7 8 9 10 > y <- 11:20 > y [1] 11 12 13 14 15 16 17 18 19 20 > z <- c(x, y) > z [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20In the above example, the vectors

> x <- seq(-2, 2, by=0.1) > length(x) [1] 41 > y <- x[18:24] > y [1] -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3In this example, we create a 41 element vector,

> x <- seq(-4, 4, by=0.01) > y <- myfunc(x) > d <- array(c(x,y), c(length(x), 2))

`d` is a 2-column matrix, with column 1 containing `x`, and
column 2 containing `y`.

> x <- d[,1] > y <- d[,2]The first column of matrix

> d <- scan("filename")

> d <- matrix(scan("filename"), ncol=2, byrow=TRUE)

> write(x, ncolumns=1, file="array.dat")

> write(t(mat1), ncolumns=3, file="matrix.dat")Notice the use of the transpose function,

> x <- seq(-4, 4, by=0.01) > y <- myfunc(x) > plot(x,y)

> d <- array(c(x,y), c(length(x), 2)) > plot(d)

> plot(d, type="l", col="red", xlab="Time (s)", ylab="Vout (V)")

Other plot attributes such as `title`, `xlim`, and `ylim`
may also be specified.

> plot(d, xlim=c(-4, 4), ylim=c(-1,1)) > lines(x, y, col="blue") > lines(d2, col="green")The first

> a <- matrix(c(1, 0, -2, 0, 3, -1), nrow=2, byrow=TRUE) > b <- matrix(c(0, 3, -2, -1, 0, 4), nrow=3, byrow=TRUE) > c <- a%*%b > c [,1] [,2] [1,] 0 -5 [2,] -6 -7 > d <- b%*%a > d [,1] [,2] [,3] [1,] 0 9 -3 [2,] -2 -3 5 [3,] 0 12 -4 >

> a <- matrix(c(1, 0, -2, 4, 1, 0, 1, 1, 7), nrow=3, byrow=TRUE) > a [,1] [,2] [,3] [1,] 1 0 -2 [2,] 4 1 0 [3,] 1 1 7 > b <- solve(a) > b [,1] [,2] [,3] [1,] 7 -2 2 [2,] -28 9 -8 [3,] 3 -1 1 >

> a <- matrix(c(1, 0, -2, 4, 1, 0, 1, 1, 7), nrow=3, byrow=TRUE) > sum(diag(a)) [1] 9The

> sigma_x <- matrix(c(0,1,1,0), nrow=2, byrow=TRUE) > sigma_x [,1] [,2] [1,] 0 1 [2,] 1 0 > eigen(sigma_x) $values [1] 1 -1 $vectors [,1] [,2] [1,] 0.7071068 -0.7071068 [2,] 0.7071068 0.7071068 >

> Jy <- matrix(c(0,1,0,-1,0,1,0,-1,0), nr=3, byrow=TRUE) > Jy [,1] [,2] [,3] [1,] 0 1 0 [2,] -1 0 1 [3,] 0 -1 0 > Jy <- Jy*(0+1i)/sqrt(2) > Jy [,1] [,2] [,3] [1,] 0+0.0000000i 0+0.7071068i 0+0.0000000i [2,] 0-0.7071068i 0+0.0000000i 0+0.7071068i [3,] 0+0.0000000i 0-0.7071068i 0+0.0000000i > eigen(Jy) $values [1] 1.000000e+00 -2.376571e-16 -1.000000e+00 $vectors [,1] [,2] [,3] [1,] -0.5+0.0000000i 7.071068e-01+0.000000e+00i -0.5+0.0000000i [2,] 0.0+0.7071068i 0.000000e+00+2.498002e-16i 0.0-0.7071068i [3,] 0.5+0.0000000i 7.071068e-01+0.000000e+00i 0.5+0.0000000i >

Visit http://cran.r-project.org and
click on *Packages*. An alphabetical list of available packages for R will
be displayed. Click on a particular package name to obtain additional information
about the package, including a pdf document which describes its use.

It is not necessary to download the package source for installation if the machine on which you installed R has internet access.

> install.packages()A list of mirror sites for R packages will appear. Select an appropriate mirror site. Then a list of available packages will appear. Select the desired package. The package will be downloaded, built, and installed.

From within the R environment, or within an R script, you may load a package using,

> require("packagename")

We will use the extra R package `deSolve` to illustrate how to solve a system of ordinary
differential equations. This package provides the robust ODE solver named `lsoda` taken from
the well-known collection of solvers, ODEPACK, written by Alan Hindmarsh and others. Assuming that
you have installed the `deSolve` package using the procedure given above, you may use it
to find a solution to the following system of first order differential equations:

du/dt = -u + alpha*v dv/dt = beta*u - vwhere alpha and beta are two specified parameters. Although the above set of coupled, linear, first-order differential equations may be solved easily by linear methods, they are used to illustrate the general numerical procedure, which may be applied to a more complex system of linear or nonlinear differential equations.

First, define a function which computes the derivatives, given the current values of the variables,
the independent variable (time, for example), and the values of the parameters:

# Function to evaluate the ODEs derivs <- function(t, x, p) { udot = -x[1] + p[1]*x[2] vdot = p[2]*x[1] - x[2] return (list(c(udot, vdot))) }The function named

u = x[1] v = x[2] alpha = p[1] beta = p[2]and

params <- c( 2, 0.1 ) # alpha = 2, beta = 0.1 x0 <- c(0, 0.05) # initial conditions times <- c(0, 0.1*(1:100)) # time steps at which to compute the solution require("deSolve") # Solve the differential equations for this example out <- lsoda( x0, times, derivs, params)The solutions,

> plot(out[,1], out[,2], type='l', col="blue", xlab="Time (s)", ylab="u(t), v(t)", ylim=c(0, .05)) > lines(out[,1], out[,3], col="red")The following graph will be shown.

The full R script for the above example is given in odesolve-example.R.