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" )
		

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.

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.

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[1], fit$coefficients[2], 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.

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[1]					# Assign the initial position to x (q0[2] is the initial velocity)
		
		xd <- q0[2]					# 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.