Ordinary Differential Equations - R
This is a program written in R to numerically solve an Ordinary Differential Equation (ODE) using Euler's Method and Runge-Kutta 4th Order Method. The values are displayed in a table and a graph to the user.
The given scenario involves a mercury thermometer system with a glass-walled design, where the temperature of the bath changes due to a chemical process. The temperature values are stored and a differential equation that describes the process is as follows: 𝜃(𝑡) = 0.1 cos 4𝑡 + 0.2 sin 4𝑡 + 2.9 𝑒 ^ (−2𝑡). The program computes discrete values of the temperature at different step sizes, calculates the relative error of the estimated solution compared to the exact solution of the observed values, and displays the results in a table and a graph. The user can choose one of the two methods of Euler's Method and Runge-Kutta 4th Order Method, and the step size for analysis.
View the Source Code below, or view the full codebase on GitHub
1 # Mostapha Abdelaziz
2 # ODE
3 #
4
5 # function that prints the first menu and validates input
6 menuOne <- function(){
7 # set variable for loop
8 menuOption <- "-1"
9
10 # print menu
11 writeLines("\nChoose the method for solving the ODE of
12 𝜃(𝑡) = 0.1 cos4𝑡 + 0.2sin4𝑡 + 2.9𝑒 ^(−2𝑡)
13 \n1. Euler's Method
14 \n2. Runge-Kutta 4th Order Method
15 \n3. Quit\n")
16 # take input
17 menuOption <- readline(prompt = "Please enter your choice: ")
18 # if it is not 1 or 2 ask again
19 while (menuOption != "1" && menuOption != "2" && menuOption != "3"){
20 writeLines("You must enter \"1\", \"2\" or \"3\", try again.")
21 menuOption <- readline(prompt = "Please enter your choice: ")
22 }
23
24 # return validated input
25 return(menuOption)
26 }
27
28 # function that prints the second menu and validates input
29 menuTwo <- function(){
30 # set variable for loop
31 menuOption <- "3"
32
33 # print menu
34 writeLines("\nChoose step size \"h\" (0.8, 0.2, 0.05)")
35 # take input
36 menuOption <- readline(prompt = "Please enter your choice: ")
37 # if it is not 1 or 2 ask again
38 while (menuOption != "0.8" && menuOption != "0.2"
39 && menuOption != "0.05" && menuOption != ".8"
40 && menuOption != ".2" && menuOption != ".05"){
41 writeLines("You must enter \"0.8\" or \"0.2\" or \"0.05\", try again.")
42 menuOption <- readline(prompt = "Please enter your choice: ")
43 }
44
45 # return validated input
46 return(menuOption)
47 }
48
49 # function for solving ODE using euler method
50 euler <- function(xVec, ynot, stepSize){
51 # create a vector for y values
52 yVec <- vector(mode = "numeric", length = length(xVec))
53 # define ynot as provided value
54 yVec[1] <- ynot
55
56 # loop through to calculate each y value
57 for(i in 1:length(xVec)){
58 # f(x,y) = cos4x - 2y y(i+1) = yi + f(x,y)h
59 yVec[i+1] <- yVec[i] + ( (cos(4*xVec[i]) - 2*yVec[i]) * stepSize )
60 }
61
62 # return our y values
63 return(yVec)
64 }
65
66 # function for solving ODE using runge-kutta 4th order
67 rungekutta <- function(xVec, ynot, stepSize){
68 # create a vector for y values
69 yVec <- vector(mode = "numeric", length = length(xVec))
70 # define ynot as provided value
71 yVec[1] <- ynot
72
73 # loop through to calculate each y value
74 for(i in 1:length(xVec)){
75 # f(x,y) = cos4x - 2y y(i+1) = yi + f(x,y)h
76 # calculate k1,k2,k3,k4
77 k1 <- (cos(4*xVec[i]) - 2*yVec[i])
78 k2 <- (cos(4*(xVec[i] + 0.5*stepSize)) - 2*(yVec[i] + 0.5*k1*stepSize))
79 k3 <- (cos(4*(xVec[i] + 0.5*stepSize)) - 2*(yVec[i] + 0.5*k2*stepSize))
80 k4 <- (cos(4*(xVec[i] + stepSize)) - 2*(yVec[i] + k3*stepSize))
81
82 # caclulate yi
83 yVec[i+1] <- yVec[i] + ( ((1/6)*(k1 + 2*k2 + 2*k3 + k4)) * stepSize )
84 }
85
86 # return our y values
87 return(yVec)
88 }
89
90 # function for calculating the solution
91 solution <- function(xVec, stepSize){
92 # create a vector for y values
93 yVec <- vector(mode = "numeric", length = length(xVec))
94
95 # loop through to calculate each y value
96 for(i in 1:length(xVec)){
97 # f(x) = 0.1cos4x + 0.2sin4x + 2.9e^(-2x)
98 yVec[i] <- 0.1*cos(4*xVec[i]) + 0.2*sin(4*xVec[i]) + 2.9*exp(-2*xVec[i])
99 }
100
101 # return our y values
102 return(yVec)
103 }
104
105 # function for calculating relative error
106 errors <- function(yEstimate, yExact){
107 yVec <- vector(mode = "numeric", length = length(yExact))
108 # loop to calculate error
109 for(i in 1:length(yExact)){
110 absError <- abs(yEstimate[i] - yExact[i])
111 yVec[i] <- abs(absError / yExact[i]) * 100
112 }
113
114 # return vector
115 return(yVec)
116 }
117
118 calculate <- function(menuOption){
119 start <- 0
120 end <- 2
121 ynot <- 3
122
123 # print the second menu
124 stepSize <- as.numeric(menuTwo())
125 # create x vector
126 xVec <- seq(from = start, to = end, by = stepSize)
127
128 # if they chose eulers or runge-kutta call the relevant function
129 if (menuOption == "1"){
130 title <- "Euler's Method Estimated Temperature vs Exact"
131 yEstimate <- euler(xVec, ynot, stepSize)
132 } else {
133 title <- "Runge Kutta's Method Estimated Temperature vs Exact"
134 yEstimate <- rungekutta(xVec, ynot, stepSize)
135 }
136
137 # calculate the exact values
138 yExact <- solution(xVec, stepSize)
139
140 # caclulate the percent error
141 errors <- errors(yEstimate, yExact)
142
143 # print table
144 writeLines("")
145 table <- matrix(c(xVec[2:length(xVec)], yExact[2:length(yExact)],
146 yEstimate[2:length(yExact)], errors[2:length(yExact)]), ncol = 4)
147 colnames(table) <- c("Time(seconds)", " Exact Temp (C)", " Estimated Temp (C)", " Percentage Error (%)")
148 prmatrix(table, rowlab=rep("",length(yExact)))
149
150 # plot a graph of the exact temperature vs the estimated temperature
151 matplot(table[,1], table[,2:3], type=c("l","p"), lty=1, pch=1 , col=c(1,2),
152 main=title, sub = "Exact Values as a black line, Estimated Values as red points",
153 xlab="Time (Seconds)", ylab="Temperature (C)")
154 }
155
156 # main function to loop through menus
157 main <- function(){
158 menuOption <- "-1"
159
160 # loop through menu indefinitely
161 while(menuOption != 3){
162 # print first menu
163 menuOption <- menuOne()
164 # continue unless user quits
165 if (menuOption != "3"){
166 calculate(menuOption)
167 } else {
168 writeLines("Quitting...")
169 }
170
171 } # end of loop
172 }
173
174 # call function to run program
175 main()
176
177