222 lines
5.5 KiB
R
222 lines
5.5 KiB
R
# Project 4 for the University of Tulsa's CS-7863 Sci-Stat Course
|
|
# Higher Order Differential Equations and Shooting Method
|
|
# Professor: Dr. McKinney, Spring 2023
|
|
# Noah L. Schrick - 1492657
|
|
|
|
## 1. Transform the second order ODE into a system of two first order ODEs
|
|
|
|
# use ode45 or rk4sys to solve for the motion of the plane pendulum
|
|
if (!require("pracma")) install.packages("pracma")
|
|
library(pracma)
|
|
tmin <- 0
|
|
tmax <- 7
|
|
theta0 <- pi/4
|
|
omega0 <- 0
|
|
y.init <- c(theta0, omega0)
|
|
g <- 9.81
|
|
L <- 1.0
|
|
|
|
## Substitution Tricks
|
|
# theta_1 = theta_prime
|
|
# theta_2 = theta
|
|
# theta_1_prime = theta_pp
|
|
# theta_2_prime = theta_prime = theta_1
|
|
|
|
## Equation Replacements
|
|
ode_1.f <- function(t, y, k){
|
|
theta <- y[1]
|
|
omega <- y[2]
|
|
g <- k[1]
|
|
L <- k[2]
|
|
omega_p <- -k[1]/k[2] * sin(theta)
|
|
theta_p <- omega
|
|
as.matrix(c(theta_p, omega_p))
|
|
}
|
|
|
|
p.params <- c(g, L)
|
|
|
|
pend.sol <- ode45(f=function(t,y){ode_1.f(t,y,k=p.params)},
|
|
y=y.init, t0=tmin, tfinal=tmax)
|
|
|
|
# a. Plot
|
|
plot(pend.sol$t,pend.sol$y[,1],type="l",col="blue",
|
|
ylim=c(-2.3,2.3), xlab="time",ylab="amplitude")
|
|
par(new=T)
|
|
lines(pend.sol$t,pend.sol$y[,2],type="l",col="red")
|
|
abline(h=0)
|
|
legend("topleft", # coordinates, "topleft" etc
|
|
c("angle","omega"), # label
|
|
lty=c(1,1), # line
|
|
lwd=c(2.5,2.5), # weight
|
|
#cex=.8,
|
|
bty="n", # no box
|
|
col=c("blue","red") # color
|
|
)
|
|
|
|
# b. Overlay a plot of the small-angle approx
|
|
plot(pend.sol$t,pend.sol$y[,1],type="l",col="blue",
|
|
ylim=c(-2.3,2.3), xlab="time",ylab="amplitude")
|
|
par(new=T)
|
|
lines(seq(tmin,tmax,len=length(pend.sol$t)),
|
|
theta0*cos(sqrt(g/L)*seq(tmin,tmax,len=length(pend.sol$t))),type="l",col="red")
|
|
abline(h=0)
|
|
legend("topleft", # coordinates, "topleft" etc
|
|
c("angle","small angle approx"), # label
|
|
lty=c(1,1), # line
|
|
lwd=c(2.5,2.5), # weight
|
|
#cex=.8,
|
|
bty="n", # no box
|
|
col=c("blue","red") # color
|
|
)
|
|
|
|
## 2. Use ode45 to solve the damped harmonic oscillator
|
|
tmin <- 0
|
|
tmax <- 10
|
|
yo <- 1
|
|
ybaro <- 0
|
|
|
|
## Substitution Tricks
|
|
# y1 = y_prime
|
|
# y_2 = y
|
|
# y_1_prime = y_pp
|
|
# y_2_prime = y_prime = y_1
|
|
|
|
# Eqs:
|
|
# m*y_1_prime + c*y1 + k*y_2 = 0
|
|
# y_1_prime = (-c*y_1 - k*y_2)/m
|
|
# y_2_prime = y_1
|
|
|
|
# a
|
|
m <- 1
|
|
c <- 1
|
|
kvar <- 2
|
|
|
|
dho.f <- function(t, y, kpass){
|
|
yn <- y[1]
|
|
ydot <- y[2]
|
|
# c k m
|
|
ypp <- (-kpass[2]*yn - kpass[3]*ydot)/kpass[1]
|
|
#ypp <- -kpass[2]*y - kpass[3]*ydot
|
|
as.matrix(c(ydot,ypp))
|
|
}
|
|
|
|
y.init <- c(yo, ybaro)
|
|
dho.params <- c(m,c,kvar)
|
|
dho.sol <- ode45(f=function(t,y){dho.f(t,y,kpass=dho.params)},
|
|
y=y.init, t0=tmin, tfinal=tmax)
|
|
|
|
# plot pos vs time
|
|
pos <- (-m*dho.sol$y[,2] - c*dho.sol$y[,1])/kvar
|
|
plot(dho.sol$t, pos ,type="o",col="blue",
|
|
ylim=c(-1.0,1.0), xlab="time",ylab="amplitude")
|
|
par(new=T)
|
|
lines(dho.sol$t,dho.sol$y[,1],type="o",col="red")
|
|
lines(dho.sol$t,dho.sol$y[,2],type="o",col="green")
|
|
legend("topright", # coordinates, "topleft" etc
|
|
c("y","yprime", "ydoubleprime"), # label
|
|
lty=c(1,1), # line
|
|
lwd=c(2.5,2.5), # weight
|
|
#cex=.8,
|
|
bty="n", # no box
|
|
col=c("blue","red", "green") # color
|
|
)
|
|
|
|
# b
|
|
m <- 1
|
|
c <- 2
|
|
k <- 1
|
|
dho.params <- c(m,c,k)
|
|
|
|
y.init <- c(yo, ybaro)
|
|
dho.params <- c(m,c,kvar)
|
|
dho.sol <- ode45(f=function(t,y){dho.f(t,y,kpass=dho.params)},
|
|
y=y.init, t0=tmin, tfinal=tmax)
|
|
|
|
# plot pos vs time
|
|
pos <- (-m*dho.sol$y[,2] - c*dho.sol$y[,1])/kvar
|
|
plot(dho.sol$t, pos ,type="o",col="blue",
|
|
ylim=c(-1.0,1.0), xlab="time",ylab="amplitude")
|
|
par(new=T)
|
|
lines(dho.sol$t,dho.sol$y[,1],type="o",col="red")
|
|
lines(dho.sol$t,dho.sol$y[,2],type="o",col="green")
|
|
legend("topright", # coordinates, "topleft" etc
|
|
c("y","yprime", "ydoubleprime"), # label
|
|
lty=c(1,1), # line
|
|
lwd=c(2.5,2.5), # weight
|
|
#cex=.8,
|
|
bty="n", # no box
|
|
col=c("blue","red", "green") # color
|
|
)
|
|
|
|
## 3. Use ode45/rk4sys to solve for the traj of a projectile thrown vertically
|
|
|
|
# a. IVP
|
|
tmin <- 0
|
|
tmax <- 2.04
|
|
xo <- 0
|
|
vo <- 10
|
|
|
|
# b. BVP
|
|
yo <- 0
|
|
ymax <- 0
|
|
|
|
library(pracma)
|
|
projectile.f <- function(t,y){
|
|
g <- 9.81
|
|
v <- y[2] # y1dot
|
|
a <- -g # could be any f, y'' = f(t,y)
|
|
matrix(c(v,a))
|
|
}
|
|
|
|
proj.obj <- function(v0, y0=0, tfinal){
|
|
# minimize w.r.t. v0
|
|
proj.sol <- ode45(projectile.f,
|
|
y=c(y0, v0), t0=0, tfinal=tfinal)
|
|
final_index <- length(proj.sol$t)
|
|
yf <- proj.sol$y[final_index,1] # want equal to right boundary
|
|
log(abs(yf)) # minimize this
|
|
}
|
|
|
|
# user specifies tfinal and yfinal for BVP
|
|
v_best <- optimize(proj.obj,
|
|
interval=c(1,100), #bisect-esque interval
|
|
tol=1e-10,
|
|
y0=0, tfinal=10) # un-optimized obj params
|
|
v_best$minimum # best v0
|
|
|
|
best.sol <- rk4sys(projectile.f, a=0, b=10, y0=c(0, v_best$minimum),
|
|
n=20) # 20 integration stepstmax
|
|
|
|
# c. Shooting method for damped oscillator with perturbation parameter
|
|
yo <- 0
|
|
y1 <- 1
|
|
ymin <- 0
|
|
ymax <- 2
|
|
|
|
pfirst <- 0.5
|
|
psec <- 0.05
|
|
|
|
# analytical sol
|
|
pthird <- 0
|
|
|
|
## 4. Position of the earth and moon
|
|
# a. Plotly
|
|
G <- 6.673e-11 # m^3 kg-1 s^-2
|
|
M_S <- 1.9891e30 # sun kg
|
|
M_E <- 5.98e24 # earth kg
|
|
M_m <- 7.32e22 # moon kg
|
|
mu_sun <- G*M_S*(86400^2)/1e9 # km^3/days^2
|
|
|
|
# t0: jan 1, 1999 00:00:00am
|
|
x0_earth <- c(-27115219762.4, 132888652547.0, 57651255508.0)/1e3 # km
|
|
v0_earth <- c(-29794.2199907, -4924.33333333,-2135.59540741)*(86400)/1e3
|
|
# m/s -> km/day
|
|
|
|
x0_moon <- c(-27083318944, 133232649728, 57770257344)/1e3 # km
|
|
v0_moon <- c(-30864.2207031, -4835.03349304, -2042.89546204)*(86400)/1e3
|
|
# m/s -> km/day
|
|
|
|
# b. Find eclipses
|
|
|
|
# c. keep orbiter at L2 Lagrange point for a year, ignoring the moon's effect
|