108 lines
3.1 KiB
R
108 lines
3.1 KiB
R
# Project 3 for the University of Tulsa's CS-7863 Sci-Stat Course
|
|
# Numerical Ordinary Differential Equations
|
|
# Professor: Dr. McKinney, Spring 2023
|
|
# Noah L. Schrick - 1492657
|
|
|
|
## 1. Approximate deriv. of sin(x) at x=pi from h=10^-1 -> 10^-20
|
|
forward.approx <- function(func, x, h){
|
|
approx <- (func(x+h)-func(x))/h
|
|
}
|
|
|
|
forward.approx.table <- matrix(nrow = 0, ncol = 3)
|
|
colnames(forward.approx.table) <- c("h", "approximation", "error")
|
|
for(h in 10^(seq(-1,-20,-1))){
|
|
approx <- forward.approx(sin, pi, h)
|
|
error <- abs(cos(pi)-approx)
|
|
forward.approx.table <- rbind(forward.approx.table, c(h, approx, error))
|
|
}
|
|
|
|
plot(abs(log(forward.approx.table[,"h"],10)), forward.approx.table[,"error"],
|
|
xlab="h [1e-x]", ylab="error (logscale)", type="o", log="y",
|
|
main = "Forward Approximation of sin(x) at x=pi")
|
|
|
|
|
|
# Repeat with central difference approx
|
|
central.diff.approx <- function(func, x, h){
|
|
approx <- (func(x+h)-func(x-h))/(2*h)
|
|
}
|
|
central.diff.table <- matrix(nrow = 0, ncol = 3)
|
|
colnames(central.diff.table) <- c("h", "approximation", "error")
|
|
for(h in 10^(seq(-1,-20,-1))){
|
|
approx <- central.diff.approx(sin, pi, h)
|
|
error <- abs(cos(pi)-approx)
|
|
central.diff.table <- rbind(central.diff.table, c(h, approx, error))
|
|
}
|
|
|
|
plot(abs(log(central.diff.table[,"h"],10)), central.diff.table[,"error"],
|
|
xlab="h [1e-x]", ylab="error (logscale)", type="o", log="y",
|
|
main = "Central Difference Approximation of sin(x) at x=pi")
|
|
|
|
## 2.
|
|
# a) Runge-Kutta
|
|
my_rk4 <- function(f, y0, t0, tfinal, dt){
|
|
|
|
}
|
|
|
|
# b) Numerically solve the decay ode and compare to Euler and RK error
|
|
decay.f <- function(t,y,k){
|
|
# define the ode model
|
|
# k given value when solver called
|
|
dydt <- -k*y
|
|
as.matrix(dydt)
|
|
}
|
|
|
|
# From class docs:
|
|
my_euler <- function(f, y0, t0, tfinal, dt){
|
|
# follow ode45 syntax, dt is extra
|
|
# y0 is initial condition
|
|
tspan <- seq(t0,tfinal,dt)
|
|
npts <- length(tspan)
|
|
y<-rep(y0,npts) # initialize, this will be a matrix for systems
|
|
for (i in 2:npts){
|
|
dy <- f(t[i-1],y[i-1])*dt # t is not used in this case
|
|
y[i] <- y[i-1] + dy
|
|
}
|
|
return(list(t=tspan,y=y))
|
|
}
|
|
|
|
## 3. Use library function ode45 to solve the decay numerically
|
|
if (!require("pracma")) install.packages("pracma")
|
|
library(pracma)
|
|
# specify intial conds and params
|
|
k <- 1.21e-4
|
|
tmin <- 0
|
|
tmax <- 10000
|
|
y.init <- 10000
|
|
|
|
decay.ode45.sol <- ode45(f=function(t,y){decay.f(t,y,k=k)},
|
|
y=y.init, t0=tmin, tfinal=tmax)
|
|
plot(decay.ode45.sol$t,decay.ode45.sol$y)
|
|
par(new=T)
|
|
lines(seq(tmin,tmax,len=50),
|
|
y.init*exp(-k*seq(tmin,tmax,len=50)))
|
|
|
|
## 4. Solve the predator-prey model numerically
|
|
|
|
# a) k1=0.01, k2=0.1, k3=0.001, k4=0.05, prey(0)=50, pred(0)=15. t=0 -> 200
|
|
|
|
# plot
|
|
|
|
# b) Use Euler and RK to solve with h=10
|
|
|
|
# plot comparing Prey solutions to ode45
|
|
|
|
# c) Use k3=0.02
|
|
|
|
## 5. Solve SIR model numerically from t=0 -> 20
|
|
|
|
# a) a=0.5, b=1, S(0)=0.9, I(0)=0.1, R(0)=0
|
|
|
|
# plot
|
|
|
|
# b) a=3
|
|
|
|
## 6. Decomp of dinitrogen pentoxygen into nitrogen dioxide and molecular oxygen
|
|
|
|
# a) k1=1.0, k2=0.5, k3=0.2,k4=1.5, [N2O5]o=1, all other IC's=0, t=0 -> 10
|
|
|
|
# b) Increase k4 to make the intermediate species -> 0 |