CS-7863-Sci-Stat-Proj-2/Schrick-Noah_Project-2.R
2023-01-31 21:51:09 -06:00

128 lines
4.4 KiB
R

# Project 2 for the University of Tulsa's CS-7863 Sci-Stat Course
# Roots and Zeros of Functions
# Professor: Dr. McKinney, Spring 2023
# Noah L. Schrick - 1492657
## Part 1: Bisection
findZeroBisect <- function(func, xl, xr, tol, maxsteps=1e6){
steps <- 0
if (func(xl) * func(xr) >= 0){
stop("Bad interval: try again with different endpoints")
}
repeat{
xm <- (xl + xr)/2
if(abs(func(xm)) < tol || steps >= maxsteps){
break
}
ifelse(func(xl) * func(xm) > 0, xl<-xm, xr<-xm)
steps <- steps + 1
}
return(c(xm, func(xm), steps))
}
# Solve 1a
fun.a <- function(x) {x^2-2}
fun.a.string <- "x^2-2=0"
sqrt2.estimate <- findZeroBisect(fun.a.string,1,2,1e-6)
# Plot 1a
if (!require("ggplot2")) install.packages("ggplot2")
library(ggplot2)
ggplot(data.frame(x=seq(1,2,.1)), aes(x)) +
stat_function(fun=fun.a, aes(col=fun.string)) +
geom_hline(aes(yintercept=0, col = "y=0"), show.legend=TRUE)+
geom_vline(aes(xintercept=sqrt2.estimate[1],
col = "bisection estimate"), show.legend=TRUE)+
ggtitle("Zero Function") +
xlab("x") + ylab("Zero Function") +
theme(text = element_text(size=20), plot.title = element_text(hjust = 0.5))
# Solve 1b
fun.b <- function(x) {-3+x+2*exp(-x)}
fun.b.string <- "-3+x+2*exp(-x)=0"
fun.b.estimate.left <- findZeroBisect(fun.b,-2.5,0,1e-6)
fun.b.estimate.right <- findZeroBisect(fun.b,1,4,1e-6)
# Plot 1b
ggplot(data.frame(x=seq(-2.5,5,.1)), aes(x)) +
stat_function(fun=fun.b, aes(col=fun.b.string)) +
geom_hline(aes(yintercept=0, col = "y=0"), show.legend=TRUE)+
geom_vline(aes(xintercept=fun.b.estimate.left[1],
col = "first bisection estimate"), show.legend=TRUE)+
geom_vline(aes(xintercept=fun.b.estimate.right[1],
col = "second bisection estimate"), show.legend=TRUE)+
ggtitle("Zero Function") +
xlab("x") + ylab("Zero Function") +
theme(text = element_text(size=20), plot.title = element_text(hjust = 0.5))
# Solve 1c
fun.c <- function(x) {x^3-sin(x)^2}
fun.c.string <- "x^3-sin(x)^2=0"
fun.c.estimate <- findZeroBisect(fun.c,0.5,1,1e-6)
# Plot 1c
ggplot(data.frame(x=seq(0.5,1,.1)), aes(x)) +
stat_function(fun=fun.c, aes(col=fun.c.string)) +
geom_hline(aes(yintercept=0, col = "y=0"), show.legend=TRUE)+
geom_vline(aes(xintercept=fun.c.estimate[1],
col = "bisection estimate"), show.legend=TRUE)+
ggtitle("Zero Function") +
xlab("x") + ylab("Zero Function") +
theme(text = element_text(size=20), plot.title = element_text(hjust = 0.5))
## Part 2: Relaxation Method
findZeroRelax <- function(g, x.guess, tol=1e-6, maxsteps=1e6){
steps <- 0
x.old <- x.guess
isConverged <- FALSE
while (!isConverged){
x.new <- g(x.old)
if(steps >= maxsteps || (abs(x.new-x.old)<tol)){
isConverged <- TRUE
}
else{
x.old <- x.new
}
steps <- steps + 1
}
return(c(x.new, g(x.new), steps))
}
# Solve 2b
fun.2b <- function(x) {2-exp(-x)}
fun.2b.fx <- function(x) {x-2+exp(-x)}
fun.2b.string <- "x-2+exp(-x)"
fun.2b.estimate <- findZeroRelax(fun.2b,0)
# Find the other root with bisection
fun.2b.bisect <- findZeroBisect(fun.2b,-2,0,1e-6)
# Plot 2b
ggplot(data.frame(x=seq(-2.5,2.5,.1)), aes(x)) +
stat_function(fun=fun.2b.fx, aes(col=fun.2b.string)) +
geom_hline(aes(yintercept=0, col = "y=0"), show.legend=TRUE)+
geom_vline(aes(xintercept=fun.2b.estimate[1],
col = "relaxation estimate"), show.legend=TRUE)+
geom_vline(aes(xintercept=fun.2b.bisect[1],
col = "bisection estimate"), show.legend=TRUE)+
ggtitle("Zero Function") +
xlab("x") + ylab("Zero Function") +
theme(text = element_text(size=20), plot.title = element_text(hjust = 0.5))
# Solve 2c
fun.2c <- function(x) {3-2*exp(-x)}
fun.2c.fx <- function(x) {x+2*exp(-x)-3}
fun.2c.string <- "x+2*exp(-x)-3"
fun.2c.estimate <- findZeroRelax(fun.2c,-.5)
fun.2c.estimate[1]
# Find the other root with bisection
fun.2c.bisect <- findZeroBisect(fun.2c,-1,0,1e-6)
# Plot 2c
ggplot(data.frame(x=seq(-1,4,.1)), aes(x)) +
stat_function(fun=fun.2c.fx, aes(col=fun.2c.string)) +
geom_hline(aes(yintercept=0, col = "y=0"), show.legend=TRUE)+
geom_vline(aes(xintercept=fun.2c.estimate[1],
col = "relaxation estimate"), show.legend=TRUE)+
geom_vline(aes(xintercept=fun.2c.bisect[1],
col = "bisection estimate"), show.legend=TRUE)+
ggtitle("Zero Function") +
xlab("x") + ylab("Zero Function") +
theme(text = element_text(size=20), plot.title = element_text(hjust = 0.5))