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.
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 "."
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" )
A product of cosine and sine functions is plotted bewteen 0 and 5 on both the x and y axis.
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)
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.
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
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
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.
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.
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.
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.
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.
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.
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.