Loading [MathJax]/jax/input/TeX/config.js

Radiation Exposure (Riemann Integral) in R & Python

Problem:
In this problem, you are asked to find the amount of radiation a person is exposed to during some period of time. For instance, let's say Sarina has moved into a new apartment. Unbeknownst to her, there is a sample of Cobalt-60 inside one of the walls of the apartment. Initially that sample had 10 MBq of activity, but she moves in after the sample has been there for 5 years. She lives in the apartment for 6 years, then leaves. How much radiation was she exposed to?
We can actually figure this out using the radioactive decay curve from above. What we want to know is her total radiation exposure from year 5 to year 11.



Total radiation exposure corresponds to the area between the two green lines at time = 5 and time = 11, and under the blue radioactive decay curve.


Let's use an approximation algorithm to estimate the area under this curve. We'll do so by first splitting up the area into equally-sized rectangles (in this case, six of them, one rectangle per year)


Once we've done that, we can figure out the area of each rectangle pretty easily. To approximate how much radiation Sarina was exposed to, we next calculate the area of each successive rectangle and then sum up the areas of each rectangle to get the total.
The technique of finding the area under a curve is called integration. This comes to us from calculus. What we're doing in this problem is an approximation of finding the integral under the curve using the summation of rectangular areas known as a Riemann integral. This approximation becomes more and more correct the smaller the width of the rectangles becomes.

Python Code
def radiationExposure(start, stop, step):
'''
Computes and returns the amount of radiation exposed
to between the start and stop times. Calls the
function f (defined for you in the grading script)
to obtain the value of the function at any point.
start: integer, the time at which exposure begins
stop: integer, the time at which exposure ends
step: float, the width of each rectangle. You can assume that
the step size will always partition the space evenly.
returns: float, the amount of radiation exposed to
between start and stop times.
'''
total = 0
n = int((stop - start) / step)
for i in range(0, n):
total += f(start + i * step) * step
return total
## Some results:
# def f(x):
# import math
# return 10*math.e**(math.log(0.5)/5.27 * x)
# radiationExposure(0, 5, 1)
# 39.10318784326239
# def f(x):
# import math
# return 200*math.e**(math.log(0.5)/14.1 * x)
# radiationExposure(72, 96, 0.4)
# 82.61081970598869

R Code
radiationExposure <- function(start, stop, step){
total = 0
n = (stop - start) / step
for (i in 0 : (n-1)){
total = total + f(start + i * step) * step
}
print(total)
}
## Some results
# f <- function(x){
# 200 * exp(log(0.5) / 14.1 * x)
# }
# radiationExposure(72, 96, 0.4)
# [1] 82.61082
# f <- function(x){
# 10 * exp(log(0.5)/5.27 * x)
# }
# radiationExposure(0, 5, 1)
# [1] 39.10319


没有评论:

发表评论