# Making 2-D Plots

## Random Dots on an X-Y Plot

The following commands generate random numbers for x and y between 0 and 1 and then plots a point at the location of the coordinate pair.

### Using R

```	n <- 1000 # Set the number of coordinate pairs to 1000
x <- runif( n, min=0, max=1) # Generate a random number with uniform distribution between 0 and 1
y <- runif( n, min=0, max=1)
plot( x, y, type="p", pch=".", col="red", main="Random Points") # Plot type is "p" point using the character "."
``` ### Using Octave

```	n = 1000 # Set the number of coordinate pairs to 1000
x = unifrnd(0, 1, n, 1) # Generate a random number with uniform distribution between 0 and 1.  The last two parameters give the number of rows and columns.
y = unifrnd(0, 1, n, 1)
plot( x, y, "r." )
title( "Random Points" )
xlabel( "x" )
ylabel( "y" )
``` # Making 3-D Plots

## Surface Plot of Regularly Spaced Points

A product of cosine and sine functions is plotted bewteen 0 and 5 on both the x and y axis.

### Using R

```	library( lattice )
x <- seq(0, 5, by=0.25)
y <- seq(0, 5, by=0.25)
grid <- expand.grid(x=x, y=x)
for( i in 1:(length(x)*length(y))) {
grid\$z[i] <- sin(grid\$x[i])*cos(2*grid\$y[i])
}
levelplot(z~x*y, grid, contour=TRUE, regions=TRUE, col.regions=rainbow(100), cuts=10)
wireframe( z~x*y, grid)
wireframe(z~x*y, grid, drape=TRUE, col.regions=rainbow(100))
#	z <- matrix( 1:(length(x)*length(y)), length(x), length(y))
#	for( i in 1:length(x) ) {
#		for( j in 1:length(y) ) {
#			z[i,j] <- sin(x[i])*cos(2*y[j])
#		}
#	}
#	filled.contour( x, y, z, col=rainbow(20) )
#	contour( x, y, z)
```

The three plots formed are a level plot, wireframe and surface plot. The level plot would look much better if a larger gird of points were used. There are additional packages for R that allows for smoothing, which would improve the quality of the level plot.   If the data is placed in a different structure, then the filled.contour function can be used instead of levelplot. This code generates the following plot. If the contour function is called it generates the second plot with no filling.

```	x <- seq(0, 5, by=0.25)
y <- seq(0, 5, by=0.25)
z <- matrix( 1:(length(x)*length(y)), length(x), length(y))
for( i in 1:length(x) ) {
for( j in 1:length(y) ) {
z[i,j] <- sin(x[i])*cos(2*y[j])
}
}
filled.contour( x, y, z, col=rainbow(20) )
contour( x, y, z)
```  ### Using Octave

```	x = linspace(0, 5, 21);		# Generate range of x values
y = linspace(0, 5, 21);		# Generate range of y values
[xx,yy] = meshgrid(x,y);	# Generate a matrix of both the x and y axis for plotting the surface
for i = 1:21
for j = 1:21
z(i,j) = sin(xx(i,j))*cos(2*yy(i,j))	# Calculate the value for each point corresponding to the coordinates stored in the matrix
endfor
endfor
contourf(x,y,z)			# Make a filled contour plot
mesh( xx, yy, z)		# Plot a wireframe mesh
surf( xx, yy, z)		# Make a surface plot
```

The following are a contour, wireframe and surface plot.   # Find Solutions for Multiple Linear Equations

## Resistor Network This network of resistors requires the use of Kirchoff's Law. The equations for this network are given above. To solve these equations the six variables must be expressed in order as a series of rows in a matrix. The constant voltages are placed on the right side of the equation and must be expressed as a vector. Finally the matrix and vector can be used to determine the solution.

```		I1 - I2 - I5 = 0
I3 - I4 + I5 = 0
I2 + I4 - Ia = 0
15 I1 + 10 I2 = 25
15 I1 - 12 I3 + 5 I5 = 0
10 I2 - 12 I4 - 5 I5 = 0
```

### Using R

Notice that when entering coefficients for the different currents, that missing currents in the equation imply their coefficients are zero.

```	a <- rbind( c(1, -1, 0, 0, -1, 0),
c(0, 0, 1, -1, 1, 0),
c(0, 1, 0, 1, 0, -1),
c(15, 10, 0, 0, 0, 0),
c(15, 0, -12, 0, 5, 0),
c(0, 10, 0, -12, -5, 0))
b <- c(0, 0, 0, 25, 0, 0)
solve(a, b)
```

Running this script gives the following results

```	I1 = 0.9411765
I2 = 1.0882353
I3 = 1.1151961
I4 = 0.9681373
I5 = -0.1470588
Ia = 2.0563725
```

### Using Octave

```	A = [1, -1, 0, 0, -1, 0;
0, 0, 1, -1, 1, 0;
0, 1, 0, 1, 0, -1;
15, 10, 0, 0, 0, 0;
15, 0, -12, 0, 5, 0;
0, 10, 0, -12, -5, 0]
B = [0;
0;
0;
25;
0;
0]
A\B		# This is the same as A divided into B or A-1*B
```

You get the same result with fewer significant digits.

# Simple Linear Regression

## Linear Fit of Noisy Data

To demonstrate how linear regression is done a series points are generated. The x values range from 0 to 100 and the y values fit an equation of y = 5 - 0.14x with superimposed noise.

### Using R

```	x <- c(0:100)		# Generate series of integers from 0 to 100
noise <- runif( 101, -5.0, 5.0 )	# Generate uniformly distributed noise with a range from -5 to 5
y <- 5 - 0.14*x + noise		# Generate the data from the line equation and the noise
fit <- lm( y ~ x )	# Perform a linear regression on the data
summary( fit )	# Report the fit statistics
plot( x, y )	# Plot the noisy data
abline( 5, -0.14, col='red' ) 	# Draw a line corresponding to the equation y = 5 - 0.14x (Make it red)
abline( fit\$coefficients, fit\$coefficients, col='blue' )	# Draw a line corresponding to the regression line (Make it blue)
```

Running the script results in the following statistics being displayed. In addition the following plot was generated.

```	lm(formula = y ~ x)

Residuals:
Min      1Q  Median      3Q     Max
-5.3680 -1.6426  0.1286  2.0272  4.8218

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  5.140862   0.560969   9.164 7.27e-15 ***
x           -0.135499   0.009692 -13.980  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 2.84 on 99 degrees of freedom
Multiple R-squared:  0.6638,    Adjusted R-squared:  0.6604
F-statistic: 195.5 on 1 and 99 DF,  p-value: < 2.2e-16
``` The red line is the original line before the noise is added. The blue line is the best fit of the noisy data. Notice the Estimate of the y-intercept is 5.14 with a standard deviation of 0.56. The estimate of the slope is -0.135 with a standard deviation of 0.0097. The R-squared value is 0.66 because the data is so noisy.

### Using Octave

```	pkg load optim
x = 0:100;	# Generate series of integers from 0 to 100
noise = unifrnd( -5.0, 5.0, 1, 101 );	# Generate uniformly distributed noise with a range from -5 to 5.  This matrix has 1 row and 101 columns
y = 5 - 0.14*x + noise;		# Generate the data from the line equation and the noise
F = [ones(101,1), x(:)];	# Set up the factors multiplied by the y-intecept (ones) and the slope (x)
[p, e_var, r, p_var, y_var] = LinearRegression(F, y');	# Perform the regression
yFit = F*p;		# Generate a line based on the regression line
yLine = F*[5, -0.14]';	# Generate a line corresponding to the original line
plot(x, y, 'ok', x, yLine, '-r', x, yFit, '-b')	# Plot the data points, actual line and regression line
```

Running the script results in the following statistics. To display them type p and p_var to show the values of the fitting parameters. In addition the following plot was generated.

```	p =
4.75512
-0.13901
p_var =
3.4966e-001
1.0438e-004
``` The red line is the original line before the noise is added. The blue line is the best fit of the noisy data. Notice the Estimate of the y-intercept is 4.76 with a standard deviation of 0.35. The estimate of the slope is -0.13901 with a standard deviation of 0.0001. This deviation is not large enough to include the value of the actual slope.

# Solve ODEs

## Projectile in Free Fall

Solving ordinary differential equations with these programs is done numerically. Using initial conditions for spatial coordinate and its higher time derivatives, progression of the coordinate is calculated with small increments in time. For more accuracy the time increments are made smaller, thus increasing the computational time for the solution. The solution is a series of numbers which can be plotted against time. A functional representation of the solution is not possible with numerical solvers, but some classes of problems have analytic solutions which can be found using algebraic solvers such as Mathematica, Maple and Maxima.

### Using R

```	library(deSolve)				# Load in the differential equation solver package (This package uses LSODA solver written in fortran.)

my.atol <- c(1e-6, 1e-6)		# Set up the tolerances for accuracy of the solution
times <- seq(from=0, to=10, by=.01)	# Generate a sequence of numbers from 0 to 10 seconds at which the DE will be solved
initval <- c(x=50, v=25)		# Set the initial position to 50 and the initial velocity to 25
parms <- c(g=9.8, m=1.0, c=1)	# Set gravity to 9.8, mass to 1.0 and drag coefficient to 1

diffeq <- function(t, q0, p) {	# t passes in times, q0 passes in initial values, p passes in parameters
g <- p["g"]					# Assign the parameter for gravity to g
m <- p["m"]					# Assign the parameter for mass to m
c <- p["c"]					# Assign the parameter for drag coefficient to c
x <- q0					# Assign the initial position to x (q0 is the initial velocity)

xd <- q0					# xd is the first time derivative and is set to the initial velocity
xdd <- -g - c*xd/m			# xdd is defined as the graviational and drag forces divided by the mass

list( c(xd, xdd) )			# Group the solution vector for velocity and acceleration into a collection
}								# This function definition makes it useable by the LSODA package

out <- as.data.frame( lsoda( initval, times, diffeq, parms, rtol=1e-4, atol=my.atol))	# Solve the differential equation

plot( out\$time, out\$x, xlab="Time (s)", ylab="Position (m)" )			# Use the solution to plot the position versus time
plot( out\$time, out\$v, xlab="Time (s)", ylab="Velocity (m/s)" )			# Use the solution to plot the velocity versus time
```  The first plot shows how the projectile reaches some maximum height and then begins to descend. The second demonstrates a decreasing velocity which eventually becomes negative. The velocity flattens out at the end because the projectile has reached terminal velocity.

### Using Octave

```	pkg load odepkg				# Load the differential equation solver package

vopt = odeset( "RelTol", 1e-4, "AbsTol", [1e-6 1e-6] )	# Set up the tolerances for accuracy of the solution
times = 0:.01:10				# Generate a sequence of numbers from 0 to 10 seconds at which the DE will be solved
initval = [50, 25]				# Initial position is 50 and the initial velocity is 25
parms = [9.8, 1.0, 1]			# Set gravity to 9.8, mass to 1.0 and drag coefficient to 1

function dy = freefall(t, y, p)
g = p(1);
m = p(2);
c = p(3);
dy = zeros(2,1);
dy(1) = y(2)
dy(2) = -g - c*dy(1)/m
end
[t,y] = ode45(@freefall, times, initval, vopt, parms)

plot( t, y(:,1), 'o' )
xlabel( "Time (s)" );
ylabel( "Position (m)" );
plot( t, y(:,2), 'o' )
xlabel( "Time (s)" );
ylabel( "Velocity (m/s)" );
```  The first plot shows how the projectile reaches some maximum height and then begins to descend. The second demonstrates a decreasing velocity which eventually becomes negative. The velocity flattens out at the end because the projectile has reached terminal velocity.