January 10, 2017
Paths of a projectile.
I had my Numerical Methods student calculate the angle that would give a ballistic projectile its maximum range, then I had them write a program that did the the same by just trying a bunch of different angles. The diagram above is what they came up with.
It made an interesting pattern that I converted into a face-plate cover for a light switch that I made using the laser at the TechShop.
Face plate cover.
Citing this post: Urbano, L., 2017. Projectile Paths, Retrieved January 20th, 2017, from Montessori Muddle: http://MontessoriMuddle.org/ .
Attribution (Curator's Code ): Via: ᔥ Montessori Muddle; Hat tip: ↬ Montessori Muddle.
Posted in Physics, ProgrammingNo Comments » - Tags: ballistics, math, numerical methods, programming, TechShop
December 18, 2016
One of the middle schoolers built a potato gun for his math class. He was looking a the mathematical relationship between the amount of fuel (hair spray) and the hang-time of the potato. To augment this work, I had my Numerical Methods class do the math and create analytical and numerical models of the projectile motion.
One of the things my students had to figure out was what angle would give the maximum range of the projectile? You can figure this out analytically by finding the function for how the horizontal distance (x) changes as the angle (theta) changes (i.e. x(theta)) and then finding the maximum of the function.
Initial velocity vector (v) and its component vectors in the x and y directions for a given angle.
Distance as a function of the angle
In a nutshell, to find the distance traveled by the potato we break its initial velocity into its x and y components (v_{x} and v_{y}), use the y component to find the flight time of the projectile (t_{f}), and then use the v_{x} component to find the distance traveled over the flight time.
Starting with the diagram above we can separate the initial velocity of the potato into its two components using basic trigonometry:
,
so,
,
Now we know that the height of a projectile (y) is given by the function:
(you can figure this out by assuming that the acceleration due to gravity (a) is constant and acceleration is the second differential of position with respect to time.)
To find the flight time we assume we’re starting with an initial height of zero (y_{0} = 0), and that the flight ends when the potato hits the ground which is also at zero ((y_{t} = 0), so:
Factoring out t gives:
Looking at the two factors, we can now see that there are two solutions to this problem, which should not be too much of a surprise since the height equation is parabolic (a second order polynomial). The solutions are when:
The first solution is obviously the initial launch time, while the second is going to be the flight time (t_{f}).
You might think it’s odd to have a negative in the equation, but remember, the acceleration is negative so it’ll cancel out.
Now since we’re working with the y component of the velocity vector, the initial velocity in this equation (v_{0}) is really just v_{y}:
so we can substitute in the trig function for v_{y} to get:
Our horizontal distance is simply given by the velocity in the x direction (v_{x}) times the flight time:
which becomes:
and substituting in the trig function for v_{x} (just to make things look more complicated):
and factoring out some of the constants gives:
Now we have distance as a function of the launch angle.
We can simplify this a little by using the double-angle formula:
to get:
Finding the maximum distance
How do we find the maxima for this function. Sketching the curve should be easy enough, but because we know a little calculus we know that the maximum will occur when the first differential is equal to zero. So we differentiate with respect to the angle to get:
and set the differential equal to zero:
and solve to get:
Since we remember that the arccosine of 0 is 90 degrees:
And thus we’ve found the angle that gives the maximum launch distance for a potato gun.
Citing this post: Urbano, L., 2016. Maximum Range of a Potato Gun, Retrieved January 20th, 2017, from Montessori Muddle: http://MontessoriMuddle.org/ .
Attribution (Curator's Code ): Via: ᔥ Montessori Muddle; Hat tip: ↬ Montessori Muddle.
Posted in Mathematics, ProgrammingNo Comments » - Tags: math, math and programming, numerical methods, physics, projectile motion
November 3, 2016
Previously, I showed how to solve a simple problem of motion at a constant velocity analytically and numerically. Because of the nature of the problem both solutions gave the same result. Now we’ll try a constant acceleration problem which should highlight some of the key differences between the two approaches, particularly the tradeoffs you must make when using numerical approaches.
The Problem
- A ball starts at the origin and moves horizontally with an acceleration of 0.2 m/s^{2}. Print out a table of the ball’s position (in x) with time (every second) for the first 20 seconds.
Analytical Solution
We know that acceleration (a) is the change in velocity with time (t):
so if we integrate acceleration we can find the velocity. Then, as we saw before, velocity (v) is the change in position with time:
which can be integrated to find the position (x) as a function of time.
So, to summarize, to find position as a function of time given only an acceleration, we need to integrate twice: first to get velocity then to get x.
For this problem where the acceleration is a constant 0.2 m/s^{2} we start with acceleration:
which integrates to give the general solution,
To find the constant of integration we refer to the original question which does not say anything about velocity, so we assume that the initial velocity was 0: i.e.:
at t = 0 we have v = 0;
which we can substitute into the velocity equation to find that, for this problem, c is zero:
making the specific velocity equation:
we replace v with dx/dt and integrate:
This constant of integration can be found since we know that the ball starts at the origin so
at t = 0 we have x = 0, so;
Therefore our final equation for x is:
Summarizing the Analytical
To summarize the analytical solution:
These are all a function of time so it might be more proper to write them as:
Velocity and acceleration represent rates of change which so we could also write these equations as:
or we could even write acceleration as the second differential of the position:
or, if we preferred, we could even write it in prime notation for the differentials:
The Numerical Solution
As we saw before we can determine the position of a moving object if we know its old position (x_{old}) and how much that position has changed (dx).
where the change in position is determined from the fact that velocity (v) is the change in position with time (dx/dt):
which rearranges to:
So to find the new position of an object across a timestep we need two equations:
In this problem we don’t yet have the velocity because it changes with time, but we could use the exact same logic to find velocity since acceleration (a) is the change in velocity with time (dv/dt):
which rearranges to:
and knowing the change in velocity (dv) we can find the velocity using:
Therefore, we have four equations to find the position of an accelerating object (note that in the third equation I’ve replaced v with v_{new} which is calculated in the second equation):
These we can plug into a python program just so:
motion-01-both.py
from visual import *
# Initialize
x = 0.0
v = 0.0
a = 0.2
dt = 1.0
# Time loop
for t in arange(dt, 20+dt, dt):
# Analytical solution
x_a = 0.1 * t**2
# Numerical solution
dv = a * dt
v = v + dv
dx = v * dt
x = x + dx
# Output
print t, x_a, x
which give output of:
>>>
1.0 0.1 0.2
2.0 0.4 0.6
3.0 0.9 1.2
4.0 1.6 2.0
5.0 2.5 3.0
6.0 3.6 4.2
7.0 4.9 5.6
8.0 6.4 7.2
9.0 8.1 9.0
10.0 10.0 11.0
11.0 12.1 13.2
12.0 14.4 15.6
13.0 16.9 18.2
14.0 19.6 21.0
15.0 22.5 24.0
16.0 25.6 27.2
17.0 28.9 30.6
18.0 32.4 34.2
19.0 36.1 38.0
20.0 40.0 42.0
Here, unlike the case with constant velocity, the two methods give slightly different results. The analytical solution is the correct one, so we’ll use it for reference. The numerical solution is off because it does not fully account for the continuous nature of the acceleration: we update the velocity ever timestep (every 1 second), so the velocity changes in chunks.
To get a better result we can reduce the timestep. Using dt = 0.1 gives final results of:
18.8 35.344 35.532
18.9 35.721 35.91
19.0 36.1 36.29
19.1 36.481 36.672
19.2 36.864 37.056
19.3 37.249 37.442
19.4 37.636 37.83
19.5 38.025 38.22
19.6 38.416 38.612
19.7 38.809 39.006
19.8 39.204 39.402
19.9 39.601 39.8
20.0 40.0 40.2
which is much closer, but requires a bit more runtime on the computer. And this is the key tradeoff with numerical solutions: greater accuracy requires smaller timesteps which results in longer runtimes on the computer.
Post Script
To generate a graph of the data use the code:
from visual import *
from visual.graph import *
# Initialize
x = 0.0
v = 0.0
a = 0.2
dt = 1.0
analyticCurve = gcurve(color=color.red)
numericCurve = gcurve(color=color.yellow)
# Time loop
for t in arange(dt, 20+dt, dt):
# Analytical solution
x_a = 0.1 * t**2
# Numerical solution
dv = a * dt
v = v + dv
dx = v * dt
x = x + dx
# Output
print t, x_a, x
analyticCurve.plot(pos=(t, x_a))
numericCurve.plot(pos=(t,x))
which gives:
Comparison of numerical and analytical solutions using a timestep (dt) of 1.0 seconds.
Citing this post: Urbano, L., 2016. Numerical and Analytical Solutions 2: Constant Acceleration, Retrieved January 20th, 2017, from Montessori Muddle: http://MontessoriMuddle.org/ .
Attribution (Curator's Code ): Via: ᔥ Montessori Muddle; Hat tip: ↬ Montessori Muddle.
Posted in Calculus, Mathematics, ProgrammingNo Comments » - Tags: motion, numerical methods, physics, physics and calculus, programming
November 3, 2016
We’ve started working on the physics of motion in my programming class, and really it boils down to solving differential equations using numerical methods. Since the class has a calculus co-requisite I thought a good way to approach teaching this would be to first have the solve the basic equations for motion (velocity and acceleration) analytically–using calculus–before we took the numerical approach.
Constant velocity
- Question 1. A ball starts at the origin and moves horizontally at a speed of 0.5 m/s. Print out a table of the ball’s position (in x) with time (t) (every second) for the first 20 seconds.
Analytical Solution:
Well, we know that speed is the change in position (in the x direction in this case) with time, so a constant velocity of 0.5 m/s can be written as the differential equation:
To get the ball’s position at a given time we need to integrate this differential equation. It turns out that my calculus students had not gotten to integration yet. So I gave them the 5 minute version, which they were able to pick up pretty quickly since integration’s just the reverse of differentiation, and we were able to move on.
Integrating gives:
which includes a constant of integration (c). This is the general solution to the differential equation. It’s called the general solution because we still can’t use it since we don’t know what c is. We need to find the specific solution for this particular problem.
In order to find c we need to know the actual position of the ball is at one point in time. Fortunately, the problem states that the ball starts at the origin where x=0 so we know that:
So we plug these values into the general solution to get:
solving for c gives:
Therefore our specific solution is simply:
And we can write a simple python program to print out the position of the ball every second for 20 seconds:
motion-01-analytic.py
for t in range(21):
x = 0.5 * t
print t, x
which gives the result:
>>>
0 0.0
1 0.5
2 1.0
3 1.5
4 2.0
5 2.5
6 3.0
7 3.5
8 4.0
9 4.5
10 5.0
11 5.5
12 6.0
13 6.5
14 7.0
15 7.5
16 8.0
17 8.5
18 9.0
19 9.5
20 10.0
Numerical Solution:
Finding the numerical solution to the differential equation involves not integrating, which is particularly good if the differential equation can’t be integrated.
We start with the same differential equation for velocity:
but instead of trying to solve it we’ll just approximate a solution by recognizing that we use dx/dy to represent when the change in x and t are really, really small. If we were to assume they weren’t infinitesimally small we would rewrite the equations using deltas instead of d’s:
now we can manipulate this equation using algebra to show that:
so the change in the position at any given moment is just the velocity (0.5 m/s) times the timestep. Therefore, to keep track of the position of the ball we need to just add the change in position to the old position of the ball:
Now we can write a program to calculate the position of the ball using this numerical approximation.
motion-01-numeric.py
from visual import *
# Initialize
x = 0.0
dt = 1.0
# Time loop
for t in arange(dt, 21, dt):
v = 0.5
dx = v * dt
x = x + dx
print t, x
I’m sure you’ve noticed a couple inefficiencies in this program. Primarily, that the velocity v, which is a constant, is set inside the loop, which just means it’s reset to the same value every time the loop loops. However, I’m putting it in there because when we get working on acceleration the velocity will change with time.
I also import the visual library (vpython.org) because it imports the numpy library and we’ll be creating and moving 3d balls in a little bit as well.
Finally, the two statements for calculating dx and x could easily be combined into one. I’m only keeping them separate to be consistent with the math described above.
A Program with both Analytical and Numerical Solutions
For constant velocity problems the numerical approach gives the same results as the analytical solution, but that’s most definitely not going to be the case in the future, so to compare the two results more easily we can combine the two programs into one:
motion-01.py
from visual import *
# Initialize
x = 0.0
dt = 1.0
# Time loop
for t in arange(dt, 21, dt):
v = 0.5
# Analytical solution
x_a = v * t
# Numerical solution
dx = v * dt
x = x + dx
# Output
print t, x_a, x
which outputs:
>>>
1.0 0.5 0.5
2.0 1.0 1.0
3.0 1.5 1.5
4.0 2.0 2.0
5.0 2.5 2.5
6.0 3.0 3.0
7.0 3.5 3.5
8.0 4.0 4.0
9.0 4.5 4.5
10.0 5.0 5.0
11.0 5.5 5.5
12.0 6.0 6.0
13.0 6.5 6.5
14.0 7.0 7.0
15.0 7.5 7.5
16.0 8.0 8.0
17.0 8.5 8.5
18.0 9.0 9.0
19.0 9.5 9.5
20.0 10.0 10.0
Solving a problem involving acceleration comes next.
Citing this post: Urbano, L., 2016. Numerical versus Analytical Solutions, Retrieved January 20th, 2017, from Montessori Muddle: http://MontessoriMuddle.org/ .
Attribution (Curator's Code ): Via: ᔥ Montessori Muddle; Hat tip: ↬ Montessori Muddle.
Posted in Calculus, Mathematics, Physics, ProgrammingNo Comments » - Tags: calculus, differential equations, motion, numerical methods, programming
September 13, 2016
Genetic algorithm trying to find a series of four mathematical operations (e.g. -3*4/7+9) that would result in the number 42.
I’m teaching a numerical methods class that’s partly an introduction to programming, and partly a survey of numerical solutions to different types of problems students might encounter in the wild. I thought I’d look into doing a session on genetic algorithms, which are an important precursor to things like networks that have been found to be useful in a wide variety of fields including image and character recognition, stock market prediction and medical diagnostics.
The ai-junkie, bare-essentials page on genetic algorithms seemed a reasonable place to start. The site is definitely readable and I was able to put together a code to try to solve its example problem: to figure out what series of four mathematical operations using only single digits (e.g. +5*3/2-7) would give target number (42 in this example).
The procedure is as follows:
- Initialize: Generate several random sets of four operations,
- Test for fitness: Check which ones come closest to the target number,
- Select: Select the two best options (which is not quite what the ai-junkie says to do, but it worked better for me),
- Mate: Combine the two best options semi-randomly (i.e. exchange some percentage of the operations) to produce a new set of operations
- Mutate: swap out some small percentage of the operations randomly,
- Repeat: Go back to the second step (and repeat until you hit the target).
And this is the code I came up with:
genetic_algorithm2.py
''' Write a program to combine the sequence of numbers 0123456789 and
the operators */+- to get the target value (42 (as an integer))
'''
'''
Procedure:
1. Randomly generate a few sequences (ns=10) where each sequence is 8
charaters long (ng=8).
2. Check how close the sequence's value is to the target value.
The closer the sequence the higher the weight it will get so use:
w = 1/(value - target)
3. Chose two of the sequences in a way that gives preference to higher
weights.
4. Randomly combine the successful sequences to create new sequences (ns=10)
5. Repeat until target is achieved.
'''
from visual import *
from visual.graph import *
from random import *
import operator
# MODEL PARAMETERS
ns = 100
target_val = 42 #the value the program is trying to achieve
sequence_length = 4 # the number of operators in the sequence
crossover_rate = 0.3
mutation_rate = 0.1
max_itterations = 400
class operation:
def __init__(self, operator = None, number = None, nmin = 0, nmax = 9, type="int"):
if operator == None:
n = randrange(1,5)
if n == 1:
self.operator = "+"
elif n == 2:
self.operator = "-"
elif n == 3:
self.operator = "/"
else:
self.operator = "*"
else:
self.operator = operator
if number == None:
#generate random number from 0-9
self.number = 0
if self.operator == "/":
while self.number == 0:
self.number = randrange(nmin, nmax)
else:
self.number = randrange(nmin, nmax)
else:
self.number = number
self.number = float(self.number)
def calc(self, val=0):
# perform operation given the input value
if self.operator == "+":
val += self.number
elif self.operator == "-":
val -= self.number
elif self.operator == "*":
val *= self.number
elif self.operator == "/":
val /= self.number
return val
class gene:
def __init__(self, n_operations = 5, seq = None):
#seq is a sequence of operations (see class above)
#initalize
self.n_operations = n_operations
#generate sequence
if seq == None:
#print "Generating sequence"
self.seq = []
self.seq.append(operation(operator="+")) # the default operation is + some number
for i in range(n_operations-1):
#generate random number
self.seq.append(operation())
else:
self.seq = seq
self.calc_seq()
#print "Sequence: ", self.seq
def stringify(self):
seq = ""
for i in self.seq:
seq = seq + i.operator + str(i.number)
return seq
def calc_seq(self):
self.val = 0
for i in self.seq:
#print i.calc(self.val)
self.val = i.calc(self.val)
return self.val
def crossover(self, ingene, rate):
# combine this gene with the ingene at the given rate (between 0 and 1)
# of mixing to create two new genes
#print "In 1: ", self.stringify()
#print "In 2: ", ingene.stringify()
new_seq_a = []
new_seq_b = []
for i in range(len(self.seq)):
if (random() < rate): # swap
new_seq_a.append(ingene.seq[i])
new_seq_b.append(self.seq[i])
else:
new_seq_b.append(ingene.seq[i])
new_seq_a.append(self.seq[i])
new_gene_a = gene(seq = new_seq_a)
new_gene_b = gene(seq = new_seq_b)
#print "Out 1:", new_gene_a.stringify()
#print "Out 2:", new_gene_b.stringify()
return (new_gene_a, new_gene_b)
def mutate(self, mutation_rate):
for i in range(1, len(self.seq)):
if random() < mutation_rate:
self.seq[i] = operation()
def weight(target, val):
if val <> None:
#print abs(target - val)
if abs(target - val) <> 0:
w = (1. / abs(target - val))
else:
w = "Bingo"
print "Bingo: target, val = ", target, val
else:
w = 0.
return w
def pick_value(weights):
#given a series of weights randomly pick one of the sequence accounting for
# the values of the weights
# sum all the weights (for normalization)
total = 0
for i in weights:
total += i
# make an array of the normalized cumulative totals of the weights.
cum_wts = []
ctot = 0.0
cum_wts.append(ctot)
for i in range(len(weights)):
ctot += weights[i]/total
cum_wts.append(ctot)
#print cum_wts
# get random number and find where it occurs in array
n = random()
index = randrange(0, len(weights)-1)
for i in range(len(cum_wts)-1):
#print i, cum_wts[i], n, cum_wts[i+1]
if n >= cum_wts[i] and n < cum_wts[i+1]:
index = i
#print "Picked", i
break
return index
def pick_best(weights):
# pick the top two values from the sequences
i1 = -1
i2 = -1
max1 = 0.
max2 = 0.
for i in range(len(weights)):
if weights[i] > max1:
max2 = max1
max1 = weights[i]
i2 = i1
i1 = i
elif weights[i] > max2:
max2 = weights[i]
i2 = i
return (i1, i2)
# Main loop
l_loop = True
loop_num = 0
best_gene = None
##test = gene()
##test.print_seq()
##print test.calc_seq()
# initialize
genes = []
for i in range(ns):
genes.append(gene(n_operations=sequence_length))
#print genes[-1].stringify(), genes[-1].val
f1 = gcurve(color=color.cyan)
while (l_loop and loop_num < max_itterations):
loop_num += 1
if (loop_num%10 == 0):
print "Loop: ", loop_num
# Calculate weights
weights = []
for i in range(ns):
weights.append(weight(target_val, genes[i].val))
# check for hit on target
if weights[-1] == "Bingo":
print "Bingo", genes[i].stringify(), genes[i].val
l_loop = False
best_gene = genes[i]
break
#print weights
if l_loop:
# indicate which was the best fit option (highest weight)
max_w = 0.0
max_i = -1
for i in range(len(weights)):
#print max_w, weights[i]
if weights[i] > max_w:
max_w = weights[i]
max_i = i
best_gene = genes[max_i]
## print "Best operation:", max_i, genes[max_i].stringify(), \
## genes[max_i].val, max_w
f1.plot(pos=(loop_num, best_gene.val))
# Pick parent gene sequences for next generation
# pick first of the genes using weigths for preference
## index = pick_value(weights)
## print "Picked operation: ", index, genes[index].stringify(), \
## genes[index].val, weights[index]
##
## # pick second gene
## index2 = index
## while index2 == index:
## index2 = pick_value(weights)
## print "Picked operation 2:", index2, genes[index2].stringify(), \
## genes[index2].val, weights[index2]
##
(index, index2) = pick_best(weights)
# Crossover: combine genes to get the new population
new_genes = []
for i in range(ns/2):
(a,b) = genes[index].crossover(genes[index2], crossover_rate)
new_genes.append(a)
new_genes.append(b)
# Mutate
for i in new_genes:
i.mutate(mutation_rate)
# update genes array
genes = []
for i in new_genes:
genes.append(i)
print
print "Best Gene:", best_gene.stringify(), best_gene.val
print "Number of iterations:", loop_num
##
When run, the code usually gets a valid answer, but does not always converge: The figure at the top of this post shows it finding a solution after 142 iterations (the solution it found was: +8.0 +8.0 *3.0 -6.0). The code is rough, but is all I have time for at the moment. However, it should be a reasonable starting point if I should decide to discuss these in class.
Citing this post: Urbano, L., 2016. Experimenting with Genetic Algorithms, Retrieved January 20th, 2017, from Montessori Muddle: http://MontessoriMuddle.org/ .
Attribution (Curator's Code ): Via: ᔥ Montessori Muddle; Hat tip: ↬ Montessori Muddle.
Posted in Environmental Science, Natural World, ProgrammingNo Comments » - Tags: genetics, numerical methods, programming, programming with VPython, vpython
April 20, 2013
Numerically integrating the area under the curve using four trapezoids.
I gave a quick introduction to programming for my calculus class, which has been working on numerical integration.
Numerical integration is usually used for functions that can’t be integrated (or not easily integrated) but for this example we’ll use a simple parabolic function so we can compare the numerical results to the analytical solution (as seen here).
With the equation:
To find the area under the curve between x = 1 and x = 5 we’d find the definite integral:
which gives the result:
For numerical integration, we break the area of concern into a number of trapezoids, find the areas of all the trapezoids and add them up.
We’ll define the left and right boundaries of the area as a and b, and we can write the integral as:
The left and right boundaries of the area we’re interested in are defined as a and b respectively. The area of each trapezoid is defined as A_{n}.
We also have to choose a number of trapezoids (n) or the width of each trapezoid (dx). Here we choose four trapezoids (n = 4), which gives a trapezoid width of one (dx = 1).
The area of the first trapezoid can be calculated from its width (dx) and the height of the two upper ends of the trapezoid (f(x_{0}) and f(x_{1}).
So if we define the x values of the left and right sides of the first trapezoids as x_{0} and x_{1}, the area of the first trapezoid is:
For this program, we’ll set the trapezoid width (dx) and then calculate the number of trapezoids (n) based on the width and the locations of the end boundaries a and b. So:
and the sum of all the areas will be:
We can also figure out that since the x values change by the same value (dx) for every trapezoid, it’s an arithmetic progression, so:
and,
so our summation becomes:
Which we can program with:
numerical_integration.py
# the function to be integrated
def func(x):
return -0.25*x**2 + x + 4
# define variables
a = 1. # left boundary of area
b = 5. # right boundary of area
dx = 1 # width of the trapezoids
# calculate the number of trapezoids
n = int((b - a) / dx)
# define the variable for area
Area = 0
# loop to calculate the area of each trapezoid and sum.
for i in range(1, n+1):
#the x locations of the left and right side of each trapezpoid
x0 = a+(i-1)*dx
x1 = a+i*dx
#the area of each trapezoid
Ai = dx * (func(x0) + func(x1))/ 2.
# cumulatively sum the areas
Area = Area + Ai
#print out the result.
print "Area = ", Area
And the output looks like
>>>
Area = 17.5
>>>
While the programming is pretty straightforward, it was a bit of a pain getting Python to work for one of my students who is running Windows 8. I still have not figured out a way to get it to work properly, so I’m considering trying to do it using Javascript.
Update
The javascript functions for numerical integration:
function numerically_integrate(a, b, dx, f) {
// calculate the number of trapezoids
n = (b - a) / dx;
// define the variable for area
Area = 0;
//loop to calculate the area of each trapezoid and sum.
for (i = 1; i <= n; i++) {
//the x locations of the left and right side of each trapezpoid
x0 = a + (i-1)*dx;
x1 = a + i*dx;
// the area of each trapezoid
Ai = dx * (f(x0) + f(x1))/ 2.;
// cumulatively sum the areas
Area = Area + Ai
}
return Area;
}
//define function to be integrated
function f(x){
return -0.25*Math.pow(x,2) + x + 4;
}
// define variables
a = 1; // left boundary of area
b = 5; // right boundary of area
dx = 1; // width of the trapezoids
// print out output
alert("Area = "+ numerically_integrate(a, b, dx, f));
This is a demonstration of a full html file that uses the function, and should work in any modern browser (download files: numerical-integration.zip).
Update 2
I've added the above javascript code to the embeddable graphs to allow it to calculate and display numerical integrals: you can change the values in the interactive graph below.
Citing this post: Urbano, L., 2013. Programming Numerical Integration with Python (and Javascript), Retrieved January 20th, 2017, from Montessori Muddle: http://MontessoriMuddle.org/ .
Attribution (Curator's Code ): Via: ᔥ Montessori Muddle; Hat tip: ↬ Montessori Muddle.
Posted in Calculus, Mathematics, ProgrammingNo Comments » - Tags: calculus, embeddable graphs, integration, javascript, math and programming, numerical methods, programming, python
July 24, 2012
The video shows 300 seconds of purely exponential growth (uninhibited), captured from the exponential growth VAMP scenario. Like the exponential growth function itself, the video starts off slowly then gets a lot more exciting (for a given value of exciting).
The modeled growth is based on the exponential growth function:
(1)
where:
- N = number of cells (or concentration of biomass);
- N_{0} = the starting number of cells;
- r = the rate constant, which determines how fast growth occurs; and
- t = time.
Finding the Rate Constant/Doubling Time (r)
You can enter either the rate constant (r) or the doubling time of the particular organism into the model. Determining the doubling time from the exponential growth equation is a nice exercise for pre-calculus students.
Let’s call the doubling time, t_{d}. When the organism doubles from it’s initial concentration the growth equation becomes:
divide through by N_{0}:
take the natural logs of both sides:
bring the exponent down (that’s one of the rules of logarithms);
remember that ln(e) = 1:
and solve for the doubling time:
Decay
A nice follow up would be to solve for the half life given the exponential decay function, which differs from the exponential growth function only by the negative in the exponent:
The UCSD math website has more details about Exponential Growth and Decay.
Finding the Growth Rate
A useful calculus assignment would be to determine the growth rate at any point in time, because that’s what the model actually uses to calculate the growth in cells from timestep to timestep.
The growth rate would be the change in the number of cells with time:
starting with the exponential growth equation:
since we have a natural exponent term, we’ll use the rule for differentiating natural exponents:
So to make this work we’ll have to define:
which can be differentiated to give:
and since N_{0} is a constant:
substituting in for u and du/dt gives:
rearranging (to make it look prettier (and clearer)):
(2)
Numerical Methods: Euler’s method
With this formula, the model could use linear approximations — like in Euler’s method — to simulate the growth of the biomass.
First we can discretize the differential so that the change in N and the change in time (t$) take on discrete values:
Now the change in N is the difference between the current value N^{t} and the new value N^{t+1}:
Now using this in our differentiated equation (Eq. 2) gives:
Which we can solve for the new biomass (N^^{t+1}):
to get:
This linear approximation, however, does introduce some error.
The approximated exponential growth curve (blue line) deviates from the analytical equation. The deviation compounds itself, getting worse exponentially, as time goes on.
Excel file for graphed data: exponential_growth.xls
VAMP
This is the first, basic but useful product of my summer work on the IMPS website, which is centered on the VAMP biochemical model. The VAMP model is, as of this moment, still in it’s alpha stage of development — it’s not terribly user-friendly and is fairly limited in scope — but is improving rapidly.
Citing this post: Urbano, L., 2012. Exponential Cell Growth, Retrieved January 20th, 2017, from Montessori Muddle: http://MontessoriMuddle.org/ .
Attribution (Curator's Code ): Via: ᔥ Montessori Muddle; Hat tip: ↬ Montessori Muddle.
Posted in Biology, Calculus, Mathematics, Natural World, videoNo Comments » - Tags: biology, calculus, imps, math, math and programming, microbiology, numerical methods, population growth, pre-Calculus
March 3, 2012
Based on
Euler's Method, this interactive graph illustrates a numerical method for solving differential equations. This approach is at the core of many sophisticated computer models of physical phenomena (like climate and weather).
If you know the equation for the slope of a curve (the red line for example),
and a point that the curve passes through, such as
, you can integrate to find the equation of the curve:
If you don't have a starting point (initial condition), you can draw a
slope field to see what the general pattern of all the possible solutions.
Even with a starting point, however, there are just times when you can't integrate the slope equation -- it's either too difficult or even impossible.
Then, what you can do is come up with an approximation of what the curve looks like by projecting along the slope from the starting point.
The program above demonstrates how it's done. This approach is called
Euler's Method, and is gives a numerical approximation rather than finding the exact, analytical solution using calculus (integration).
So why use an approximation when you can find the exact solution? Because, there are quite a number of problems that are impossible or extremely difficult to solve analytically, things like: the diffusion of pollution in a lake; how changing temperature in the atmosphere gives you weather and climate; the flow of groundwater in aquifers; stresses on structural members of buildings; and the list goes on and on.
As with most types of numerical approximations, you get better results if you can reduce the step size between projections of the slope. Try changing the numbers and see.
A more detailed version, with solutions, is here:
Euler's Method.
A good reference:
Euler's Method by Paul Dawkins.
Citing this post: Urbano, L., 2012. Euler's Method for Approximating a Solution to a Differential Equation, Retrieved January 20th, 2017, from Montessori Muddle: http://MontessoriMuddle.org/ .
Attribution (Curator's Code ): Via: ᔥ Montessori Muddle; Hat tip: ↬ Montessori Muddle.
Posted in Calculus, Mathematics, TechnologyNo Comments » - Tags: calculus, differential equations, graphing interactives, interactive post, math, numerical methods, online applications, parabolas