Runge-Kutta Radioactivity Most of us are familiar with carbon-14, the naturally occurring, radioactive isotope of carbon used in radiocarbon dating, but few know of its less-useful cousin, carbon15. In contrast to the relatively long-lasting carbon-14, which has a half-life of 5,730 years, carbon-15 has a half-life of only 2.45 seconds. The amount of carbon-15 over time is given by the following decay equation: dy/dt = In(2)/t1/2 y
where y is the remaining amount of carbon-15 and t1/2 is the half-life in seconds. At time t = 0, the initial amount of carbon-15, yo, is
y(t) = yo exp (-In(2)/t1/2 y
We will use the exact solution to evaluate and compare the accuracy of the three different numerical methods to capture this rapid decay process.
RK1 Atf(tkyk) = yk +1 - Atf(tk, yk)
RK2 = Atf (tx + At, yk + c) = yk + 2 - Atf(tkyk) - Atf (tx + At, yk + })
RK4 - Atf (tk + At, Yk + }) - Atf (tx + At, yx + ) Jk+1 = k + c + c2 + 3^3 + 3c4 A
s discussed in class, we will compare and contrast first order (RK1, aka explicit Euler), second-order (RK2), and fourth order (RK4) Runge-Kutta methods, which generally rely on braking a single time-step down into multiple, incremental calcula- tions in order to reduce overall error.
(a) Write a single function containing all three methods to advance the discretized solution by one time-step. The function must use exactly the function header shown below, where the function name and number and order of inputs and outputs must be followed exactly.
function [y] = advanceRK( y, dt, method)
where dt is the time-step size and method is either 1, 2, or 4 to choose whether to calculate the value at yk+1 using the first, second, or fourth order approxima- tion, respectively. The current amount of carbon-15, y, is used to calculate the dt RKI RK4 value at the next time-step, y = y + ..., and the update value is returned as the single output
Note: Although all three methods are included in a single function, only the result generated by the method option will be calculated with a single function call. t1/2 is a constant in this problem, not to be confused with the time-step size At, or time t. Establish a clear variable naming convention and use it consistently.
(b) Write a program that repeatedly calls your advanceRk function to solve Equa- tion 2 up to t final = 15 seconds using each of three different numerical meth- ods (you will need to call your function 3 times, specifying a different method each time) using three different time-step values, At = [1, 0.1, 0.01). Gener- ate 3 plots, one for each value of At, containing four curves each, 3 numerical method solutions and one exact solution. Additionally, print the average error associated with each method using exactly the format shown below (values shown are for illustration purposes only):
RK2 1.00: 3.14e-02 3.682-03 1.64e-05 0.10: 3.24e-033.262-05 1.44e-09 0.01: 3.24e-04 3.21e-07 1.42e-13
The values in the table are the mean values of the absolute errors between the numerical approximation and the exact solution given by Equation 2. For ex- ample, to calculate the average RK4 error for a given At, calculate the vector of differences, RK4 - exact, ignore negative signs, and calculate their average.
Note: You can open a new plotting window with the keyword figure to ensure that plots 1, 2, and 3 do not simply overwrite one another.
(c) Using your most accurate numerical method, what percentage of the original amount of carbon-15 remains after just 15 seconds? Just by examining the scaling of the average errors, explain how you can verify that the three different methods incur error that scales like (dt), (dt), and (dt); that is, first, second, and fourth-order, respectively? What happens to all three solutions (particularly RK1 / Euler) when we take a very large time-step, i. e. At = 5s? Why does this occur?