MAS.864 Week 4: Numerical Solutions, ODEs
Benjamin Preis
3/5/2020
library(knitr)
library(kableExtra)
library(ggplot2)
library(ggpubr)## Loading required package: magrittr7.1
The Heur method is given as:
\[ y(x+h)=y(x)+\frac{h}{2}\{f(x,y(x))+f[x+h,y(x)+hf(x,y(x))]\} \] \[ =y(x)+\frac{h}{2}f(x,y(x))+\frac{h}{2}f[x+h,y(x)+hf(x,y(x))] \]
Given that the third term contains an Euler step, we can use a Taylor expansion to estimate this term…
\[ \frac{h}{2}f[x+h,y(x)+hf(x,y(x))]=\frac{h}{2}[f(x,y(x))+h\frac{d}{dh}f(x+h,y(x)+hf(x,y(x)))|_{h=0}+\mathcal{O}[h^2]] \] \[ \frac{h}{2}[f(x,y(x))+h[1\frac{df}{dh}+hf(x,y(x))\frac{df}{dy}]+\mathcal{O}[h^2]] \]
Plugging back in…
\[ y(x+h)=y(x)+\frac{h}{2}f(x,y(x))+\frac{h}{2}[f(x,y(x))+h[1\frac{df}{dh}+hf(x,y(x))\frac{df}{dy}]+\mathcal{O}[h^2]] \]
\[ =y(x)+hf(x,y(x))+\frac{h^2}{2}(\frac{df}{dx}+f(x,y(x))\frac{df}{dy})+\mathcal{O}[h^3] \]
Therefore, the error is \(\mathcal{O}[h^3]\), in spite of the fact that we only do a 2nd order Euler approximation.
7.2
Euler Method
#Trying multiple different step sizes...
h<-c(.0001,.001,.01,.1,.25)
num_steps<-round(100*pi/h)
Y1<-list()
Y2<-list()
for(j in 1:length(h)){
num_steps<-round((100*pi)/h[j])
y1<-c()
y2<-c()
y1[1]<-1
y2[1]<-0
for(i in 1:num_steps){
y1[i+1]<-y1[i]+h[j]*y2[i]
y2[i+1]<-y2[i]-h[j]*y1[i]
}
Y1[[j]]<-data.frame("y1"=y1[2:length(y1)],"exact"=cos((0:(num_steps-1))*h[j]))
Y2[[j]]<-data.frame("y2"=y2[2:length(y1)],"exact"=-sin((0:(num_steps-1))*h[j]))
}
Y1<-lapply(Y1,function(x) cbind(x,x$exact-x$y1))
Y2<-lapply(Y2,function(x) cbind(x,x$exact-x$y2))
mean_error_y1<-unlist(lapply(Y1,function(x) mean(abs(x[,3]))))
mean_error_y2<-unlist(lapply(Y2,function(x) mean(abs(x[,3]))))
mean_y1_error<-data.frame("h"=h,"mean_error_y1"=mean_error_y1)
kable(mean_y1_error,col.names = c("Step Size","Mean Error, Y"),title="Y1,Euler")%>%
kable_styling(full_width = F)| Step Size | Mean Error, Y |
|---|---|
| 0.0001 | 5.027700e-03 |
| 0.0010 | 5.273770e-02 |
| 0.0100 | 9.077604e-01 |
| 0.1000 | 2.469614e+05 |
| 0.2500 | 6.083680e+14 |
mean_y2_error<-data.frame("h"=h,"mean_error_y2"=mean_error_y2)
kable(mean_y2_error,col.names = c("Step Size","Mean Error, Y''"),title="Y2, Euler")%>%
kable_styling(full_width = F)| Step Size | Mean Error, Y’’ |
|---|---|
| 0.0001 | 5.027300e-03 |
| 0.0010 | 5.273470e-02 |
| 0.0100 | 9.078058e-01 |
| 0.1000 | 2.547900e+05 |
| 0.2500 | 5.849678e+14 |
final_error_y1<-unlist(lapply(Y1,function(x) x[nrow(x),3]))
final_error_y2<-unlist(lapply(Y2,function(x) x[nrow(x),3]))
final_y1_error<-data.frame("h"=h,"final_error_y1"=final_error_y1)
kable(final_y1_error,col.names = c("Step Size","Final Error, Y"),title="Y1,Euler, Final Error")%>%
kable_styling(full_width = F)| Step Size | Final Error, Y |
|---|---|
| 0.0001 | -1.583200e-02 |
| 0.0010 | -1.700893e-01 |
| 0.0100 | -3.809932e+00 |
| 0.1000 | -3.321023e+06 |
| 0.2500 | -3.522898e+16 |
final_y2_error<-data.frame("h"=h,"final_error_y2"=final_error_y2)
kable(final_y2_error,col.names = c("Step Size","Final Error, Y''"),title="Y2, Euler, Final Error")%>%
kable_styling(full_width = F)| Step Size | Final Error, Y’’ |
|---|---|
| 0.0001 | 9.950000e-05 |
| 0.0010 | 8.323000e-04 |
| 0.0100 | -3.756890e-02 |
| 0.1000 | -5.176277e+06 |
| 0.2500 | 2.190518e+15 |
Runge-Kutta
#Define the Runge-Kutta Function
RK<-function(Y1,Y2,H){
k1.y1<-H*Y2
k1.y2<-(-H)*Y1
k2.y1<-H*(Y2+k1.y2/2)
k2.y2<-(-H)*(Y1+k1.y1/2)
k3.y1<-H*(Y2+k2.y2/2)
k3.y2<-(-H)*(Y1+k2.y1/2)
k4.y1<-H*(Y2+k3.y2)
k4.y2<-(-H)*(Y1+k3.y1)
Y1.RK<-(Y1+k1.y1/6+k2.y1/3+k3.y1/3+k4.y1/6)
Y2.RK<-(Y2+k1.y2/6+k2.y2/3+k3.y2/3+k4.y2/6)
temp<-c(Y1.RK,Y2.RK)
names(temp)<-c("y1","y2")
return(temp)
}
Y1RK<-list()
Y2RK<-list()
h<-c(.0001,.001,.01,.1,.25,.5,1)
for(j in 1:length(h)){
num_steps<-round((100*pi)/h[j])
y1<-c()
y2<-c()
y1[1]<-1
y2[1]<-0
for(i in 1:num_steps){
temp<-RK(y1[i],y2[i],h[j])
y1[i+1]<-temp["y1"]
y2[i+1]<-temp["y2"]
}
Y1RK[[j]]<-data.frame("y1"=y1[1:length(y1)],"exact"=cos((0:(num_steps))*h[j]))
Y2RK[[j]]<-data.frame("y2"=y2[1:length(y1)],"exact"=-sin((0:(num_steps))*h[j]))
}
Y1RK<-lapply(Y1RK,function(x) cbind(x,x$exact-x$y1))
Y2RK<-lapply(Y2RK,function(x) cbind(x,x$exact-x$y2))
mean_error_y1rk<-unlist(lapply(Y1RK,function(x) mean(abs(x[,3]))))
mean_error_y2rk<-unlist(lapply(Y2RK,function(x) mean(abs(x[,3]))))
mean_y1_errorrk<-data.frame("h"=h,"mean_error_y1"=mean_error_y1rk)
kable(mean_y1_errorrk,col.names = c("Step Size","Mean Error, Y1"),title="Y1, Mean Error, RK")%>%
kable_styling(full_width = F)| Step Size | Mean Error, Y1 |
|---|---|
| 0.0001 | 0.0000000 |
| 0.0010 | 0.0000000 |
| 0.0100 | 0.0000000 |
| 0.1000 | 0.0000833 |
| 0.2500 | 0.0032518 |
| 0.5000 | 0.0508693 |
| 1.0000 | 0.4550929 |
mean_y2_errorrk<-data.frame("h"=h,"mean_error_y2"=mean_error_y2rk)
kable(mean_y2_errorrk,col.names = c("Step Size","Mean Error, Y"),title="Y2, Mean Error, RK")%>%
kable_styling(full_width = F)| Step Size | Mean Error, Y |
|---|---|
| 0.0001 | 0.0000000 |
| 0.0010 | 0.0000000 |
| 0.0100 | 0.0000000 |
| 0.1000 | 0.0000833 |
| 0.2500 | 0.0032519 |
| 0.5000 | 0.0507112 |
| 1.0000 | 0.4533847 |
final_error_y1rk<-unlist(lapply(Y1RK,function(x) x[nrow(x),3]))
final_error_y2rk<-unlist(lapply(Y2RK,function(x) x[nrow(x),3]))
final_y1_errorrk<-data.frame("h"=h,"final_error_y1"=final_error_y1rk)
kable(final_y1_errorrk,col.names = c("Step Size","Final Error, Y1"),title="Y1, Final Error, RK")%>%
kable_styling(full_width = F)| Step Size | Final Error, Y1 |
|---|---|
| 0.0001 | 0.0000000 |
| 0.0010 | 0.0000000 |
| 0.0100 | 0.0000000 |
| 0.1000 | 0.0000112 |
| 0.2500 | 0.0012489 |
| 0.5000 | 0.0954130 |
| 1.0000 | 1.0362849 |
final_y2_errorrk<-data.frame("h"=h,"final_error_y2"=final_error_y2rk)
kable(final_y2_errorrk,col.names = c("Step Size","Final Error, Y1"),title="Y1, Final Error, RK")%>%
kable_styling(full_width = F)| Step Size | Final Error, Y1 |
|---|---|
| 0.0001 | 0.0000000 |
| 0.0010 | 0.0000000 |
| 0.0100 | 0.0000000 |
| 0.1000 | -0.0002616 |
| 0.2500 | -0.0101357 |
| 0.5000 | -0.1255992 |
| 1.0000 | 0.0203453 |
7.3
end<-100*pi
adaptiveRK.y1<-list()
adaptiveRK.y2<-list()
desired_error<-c(.000001,.00001,.0001,.001,.01,.1,1)
for(j in 1:length(desired_error)){
i<-1
t<-0
time<-c()
time[1]<-0
actual_h<-c()
h<-.5
actual_h[1]<-h
y1<-c()
y2<-c()
y1[1]<-1
y2[1]<-0
while(t<end){
temp<-RK(y1[i],y2[i],h)
#Take two half steps
temp.5<-RK(y1[i],y2[i],(h/2))
temp.5<-RK(temp.5[1],temp.5[2],(h/2))
#Adapt step size based only on position...
if(abs(temp[1]-temp.5[1]) > desired_error[j]){
h<-h*.9
} else{
h<-h*1.1
y1[i+1]<-temp.5["y1"]+((temp.5["y1"]-temp["y1"])/15)
y2[i+1]<-temp.5["y2"]+((temp.5["y2"]-temp["y2"])/15)
i<-i+1
t<-t+(i*h)
time[i]<-t
actual_h[i]<-h
}
}
adaptiveRK.y1[[j]]<-cbind(y1,actual_h)
adaptiveRK.y2[[j]]<-cbind(y2,actual_h)
}
mean_h<-unlist(lapply(adaptiveRK.y1,function(x) mean(x[,2],na.rm=T)))
kable(data.frame(desired_error,mean_h),
col.names = c("Desired Errors","Mean Step Size"),
title="Relationship Between Step Size & Desired Error, Adaptive RK")%>%
kable_styling(full_width = F)| Desired Errors | Mean Step Size |
|---|---|
| 1e-06 | 0.2023960 |
| 1e-05 | 0.3104956 |
| 1e-04 | 0.4801989 |
| 1e-03 | 0.7299420 |
| 1e-02 | 1.0101682 |
| 1e-01 | 1.2720473 |
| 1e+00 | 1.3462919 |
7.4
RK.pendulum<-function(Y1,Y2,H,A,w,t,l,g){
k1.y1<-H*Y2
k1.y2<-(-H)*((g+A*sin(w*t))*sin(Y1)/l)
k2.y1<-H*(Y2+k1.y2/2)
k2.y2<-(-H)*((g+A*sin(w*(t+h/2)))*sin(Y1+k1.y1/2)/l)
k3.y1<-H*(Y2+k2.y2/2)
k3.y2<-(-H)*((g+A*sin(w*(t+h/2)))*sin(Y1+k2.y1/2)/l)
k4.y1<-H*(Y2+k3.y2)
k4.y2<-(-H)*((g+A*sin(w*(t+h)))*sin(Y1+k3.y1)/l)
Y1.RK<-(Y1+k1.y1/6+k2.y1/3+k3.y1/3+k4.y1/6)
Y2.RK<-(Y2+k1.y2/6+k2.y2/3+k3.y2/3+k4.y2/6)
temp<-c(Y1.RK,Y2.RK)
names(temp)<-c("y1","y2")
return(temp)
}
time<-0
time<-c()
time[1]<-0
actual_h<-c()
pendulum.y1<-list()
pendulum.y2<-list()
desired_error<-.00001
h<-.5
actual_h[1]<-h
i<-1
pendulums<-data.frame("A"=c(.5,1,2.5,5),"W"=c(2,4,6,8))
end<-pi*1000
for(j in 1:nrow(pendulums)){
i<-1
t<-0
time<-c()
time[1]<-0
actual_h<-c()
h<-.5
actual_h[1]<-h
y1<-c()
y2<-c()
y1[1]<-pi/4
y2[1]<-0
while(t<end){
temp<-RK.pendulum(y1[i],y2[i],h,A=pendulums[j,1],w=pendulums[j,2],t,l=1,g=(9.8))
#Take two half steps
temp.5<-RK.pendulum(y1[i],y2[i],(h/2),A=pendulums[j,1],w=pendulums[j,2],t,l=1,g=(9.8))
temp.5<-RK.pendulum(temp.5[1],temp.5[2],(h/2),A=pendulums[j,1],w=pendulums[j,2],t,l=1,g=(9.8))
#Adapt step size based only on position...
if(abs(temp[1]-temp.5[1]) > desired_error){
h<-h*.9
} else{
h<-h*1.1
y1[i+1]<-temp.5["y1"]+((temp.5["y1"]-temp["y1"])/15)
y2[i+1]<-temp.5["y2"]+((temp.5["y2"]-temp["y2"])/15)
i<-i+1
t<-t+(i*h)
time[i]<-t
actual_h[i]<-h
}
}
pendulum.y1[[j]]<-cbind(y1,actual_h,time)
pendulum.y2[[j]]<-cbind(y2,actual_h,time)
}
pendulum.y1<-lapply(pendulum.y1,function(x) cbind(x,1:nrow(x)))
pendulum.y2<-lapply(pendulum.y2,function(x) cbind(x,1:nrow(x)))p1<-ggplot(as.data.frame(pendulum.y1[[1]]),aes(x=time,y=y1),group=1)+geom_line()+ggtitle("A=.5,W=2")+xlim(0,3000)
p2<-ggplot(as.data.frame(pendulum.y1[[2]]),aes(x=time,y=y1),group=1)+geom_line()+ggtitle("A=1,W=4")+xlim(0,3000)
p3<-ggplot(as.data.frame(pendulum.y1[[3]]),aes(x=time,y=y1),group=1)+geom_line()+ggtitle("A=2.5,W=6")+xlim(0,3000)
p4<-ggplot(as.data.frame(pendulum.y1[[4]]),aes(x=time,y=y1),group=1)+geom_line()+ggtitle("A=5,W=8")+xlim(0,3000)
ggarrange(p1,p2,p3,p4,nrow=2,ncol=2)## Warning: Removed 7 rows containing missing values (geom_path).
## Warning: Removed 7 rows containing missing values (geom_path).## Warning: Removed 9 rows containing missing values (geom_path).## Warning: Removed 12 rows containing missing values (geom_path).