Exponential Cell Growth

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:

 N = N_0 e^{rt} (1)

where:

  • N = number of cells (or concentration of biomass);
  • N0 = 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, td. When the organism doubles from it’s initial concentration the growth equation becomes:

 2N_0 = N_0 e^{r t_d}

divide through by N0:

 2  =  e^{r t_d}

take the natural logs of both sides:

 \ln 2  =  \ln (e^{r t_d})

bring the exponent down (that’s one of the rules of logarithms);

 \ln 2  =  r t_d \ln (e)

remember that ln(e) = 1:

 \ln 2  =  r t_d

and solve for the doubling time:

 \frac{\ln 2}{r}  =  t_d

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:

 N = N_0 e^{-rt}

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:

 \frac{dN}{dt}

starting with the exponential growth equation:

 N = N_0 e^{rt}

since we have a natural exponent term, we’ll use the rule for differentiating natural exponents:

 \frac{d}{dx}(e^u) = e^u \frac{du}{dx}

So to make this work we’ll have to define:

 u = rt

which can be differentiated to give:

 \frac{du}{dt} = r

and since N0 is a constant:

 N = N_0 e^{u}

 \frac{dN}{dt} = N_0 e^{u} \frac{du}{dt}

substituting in for u and du/dt gives:

 \frac{dN}{dt} = N_0 e^{rt} (r)

rearranging (to make it look prettier (and clearer)):

 \frac{dN}{dt} = N_0 r e^{rt} (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:
 \frac{dN}{dt} = \frac{\Delta N}{\Delta t}

Now the change in N is the difference between the current value Nt and the new value Nt+1:

Now using this in our differentiated equation (Eq. 2) gives:

 \frac{N^{t+1}-N^t}{\Delta t} = N_0 r e^{r\Delta t}

Which we can solve for the new biomass (N^t+1):

  N^{t+1}-N^t = N_0 r e^{r\Delta t} \Delta t

to get:
  N^{t+1}     = N_0 r e^{r\Delta t} \Delta t + N^t

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.
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.

Bending a Soccer Ball

Students from the University of Leicester have published a beautiful short research paper (pdf) on the physics of curving a soccer ball through the air.

It has been found that the amount a football bends depends linearly on the speed of the ball and the amount of spin.

— Sandhu et al., 2011: How to score a goal (pdf) in the University of Leicester’s Journal of Physics Special Topics

They derive the relationship from Bernoulli’s equation using some pretty straightforward algebra. The force (F) perpendicular to the ball’s motion that causes it to curl is:

F = 2 \pi R^3 \rho \omega v

and the distance the ball curls can be calculated from:

D = \frac{\pi R^3 \rho \omega}{ v m } x^2

where:

  • F = force perpendicular to the direction the ball is kicked
  • D = perpendicular distance the ball moves to the direction it is kicked (the amount of curl)
  • R = radius of the ball
  • ρ = density of the air
  • ω = angular velocity of the ball
  • v = velocity of the ball (in the direction it is kicked)
  • m = mass of the ball
  • x = distance traveled in the direction the ball is kicked

The paper itself is an excellent example of what a short, student research paper should look like. And there are number of neat followup projects that advanced, high-school, physics/calculus students could take on, such as: considering the vertical dimension — how much time it take for the ball to rise and fall over the wall; creating a model (VPython) of the motion of the ball; and adding in the slowing of the ball due to air friction.

ScienceDaily

Draining a Bottle Part 2: Linearizing Equations when you have to

Yesterday we used calculus to find the equation for the height of water in a large plastic water bottle as the water drained out of a small hole in the bottom.

Perhaps the most crucial point in the procedure was fitting a curve to the measured reduction of the water’s outflow rate over time. Yesterday, in our initial attempt, we used a straight line for the curve, which produced a very good fit.

Figure 1. The change in the outflow rate over time can be well approximated by a straight line.

The R2 value is a measure of how good a fit the data is to the trendline. The straight line gives an R2 value of 0.9854, which is very close to a perfect fit of 1.0 (the lowest R2 can go is 0.0).

The resulting equation, written in terms of the outflow rate (dV/dt) and time (t), was:

 \frac{dV}{dt} = -0.0035 t + 3.9113

However, if you look carefully at the graph in Figure 1, the last few data points suggest that the outflow does not just linearly decrease to zero, but approaches zero asymptotically. As a result, a different type of curve might be a better trendline.

Types of Equations

So my calculus students and I, with a little help from the pre-Calculus class, tried to figure out what types of curves might work. There are quite a few, but we settled for looking at three: a logarithmic function, a reciprocal function, and a square root function. These are shown in Figure 2.

Figure 2. Example curves that might better describe the relationship between outflow and time.

I steered them toward the square root function because then we’d end up with something akin to Torricelli’s Law (which can be derived from the physics). A basic square root function for outflow would look something like this:

 \frac{dV}{dt} = a \sqrt{t} + b

the a coefficient stretches the equation out, while the b coefficient moves the curve up and down.

Fitting the Curve

Having decided on a square-root type function, the next problem was trying to find the actual equation. Previously, we used Excel to find the best fit straight line. However, while Excel can fit log, exponential and power curves, there’s no option for fitting a square-root function to a graph.

To get around this we linearized the square-root function. The equation, after all, looks a lot like the equation of a straight line already, the only difference is the square root of t, so let’s substitute in:

 x = \sqrt{t}

to get:

 \frac{dV}{dt} = a x + b

Now we can get Excel to fit a straight line to our data, but we have to plot the square-root of time versus temperature instead of the just time versus temperature. So we take the square root of all of our time measurements:

Time Square root of time Outflow rate
t (s) t1/2 = x (s1/2) dV/dt (ml/s)
0.0 0 3.91
45.5 6.75 3.52
97.8 9.89 2.94
140.9 11.87 3.52
197 14.04 3.21
257 16.05 3.01
315.1 17.75 2.81
380.1 19.50 2.53
452.9 21.28 2.23
529.6 23.01 1.92
620.7 24.91 1.69
742.7 27.25 1.45

We can now plot the outflow rate versus the square root of time (Figure 3).

Figure 3. Linear trend relating the outflow rate to the square root of time. The regression coefficient (R2) of 0.9948 is better than the simply linear trend of outlfow rate versus time (which was 0.9854).

The equation Excel gives (Figure 3), is:

 \frac{dV}{dt} = -0.1395 x + 5.21

and we can substitute back in for x=t1/2 to get:

 \frac{dV}{dt} = -0.1395 \sqrt{t} + 5.21

Getting back to the Equation for Height

Now we can do the same procedure we did before to find the equation for height.

First we substitute in V=πr2h:

 \frac{d(\pi r^2 h}{dt} = -0.1395 \sqrt{t} + 5.21

Factor out the πr2 and move it to the other side of the equation to solve for the rate of change of height:

 \frac{dh}{dt} = \frac{-0.1395}{\pi r^2} \sqrt{t} + \frac{5.21}{\pi r^2}

Then integrate to find h(t) (remember \sqrt{t} = t^{1/2} ) :

 \int \frac{dh}{dt} dt = \int \left( \frac{-0.1395}{\pi r^2} t^{1/2} + \frac{5.21}{\pi r^2} \right) dt

gives:

 h =  \frac{-0.1395}{(3/2) \pi r^2} t^{3/2} + \frac{5.21}{\pi r^2} t + c

which might look a bit ugly, but that’s only because I haven’t simplified the fractions. Since the radius (r) is 7.5 cm:

 h =  -0.000526 t^{3/2} + 0.029 t + c

Finally we substitute in the initial value (t=0, h=11) to solve for the coefficient:

 c = 11

giving the equation:

 h =  -0.000526 t^{3/2} + 0.029 t + 11

Plotting the equations shows that it matches the measured data fairly well, although not quite as well as when we used the previous linear function for outflow.

Figure 4. Integrating a square root function for the outflow rate gives a modeled function for the changing height over time that slightly undermatches the measured heights.

Discussion

I’m not sure why the square root function for outflow does not give as good a match of the measurements of height as does the linear function, especially since the former better matches the data (it has a better R2 value).

It could be because of the error in the measurements; the gradations on the water bottle were drawn by hand with a sharpie so the error in the height measurements there alone was probably on the order of 2-3 mm. The measurement of the outflow volume in the beaker was also probably off by about 5%.

I suspect, however, that the relatively short time for the experiment (about 15 minutes) may have a large role in determining which model fit better. If we’d run the experiment for longer, so students could measure the long tail as the water height in the bottle got close to the outlet level and the outflow rate really slowed down, then we’d have found a much better match using the square-root function. The linear match of the outflow data produces a quadratic equation when you integrate it. Quadratic equations will drop to a minimum and then rise again, unlike the square-root function which will just continue to sink.

Conclusions

The linearization of the square-root function worked very nicely. It was a great mathematical example even if it did not produce the better result, it was still close enough to be worth it.

The Draining of a Plastic Bottle: Integrating a Physics Experiment into Calculus

Abstract

I punched a small hole (about 1mm radius) in a one gallon plastic bottle and had my students measure the rate at which water drained. Even though the apparatus and measurement technique was fairly rough, we were able to, with a little calculus, determine the equation for the height of the water in the bottle as a function of time.

Figure 1. A student uses a stopwatch to measure the outflow rate of water from the plastic bottle.

Introduction

Questions about water draining from a tank are a pretty common in calculus textbooks, but there is a significant difference between seeing the problem written down, and having to figure it out from a physical example. The latter is much more challenging because it does not presuppose any relationships for the change in the height of water with time; students must determine the relationship from the data they collect.

The experimental approach mimics the challenges faced by scientists such as Henry Darcy who first determined the formulas for groundwater flow (Darcy, 1733) almost 300 years ago, not long after the development of modern calculus by Newton and Leibniz (O’Conner and Robertson, 1996).

Procedure

Figure 2. The apparatus.

I punched a small hole, about 1mm in radius, in the base of a plastic, one-gallon bottle. I chose this particular type of bottle because the bulk of it was cylindrical in shape.

Students were instructed to figure out how the rate at which water flowed out (outflow rate) changed with time, and how the height of the water (h) in the bottle changed with time. These relationships would allow me to predict the outflow rate at any time, as well as how much water was left in the container at any time.

Data collection:

  1. To measure height versus time, we marked the side of the bottle (within the cylindrical region) in one centimeter increments and recorded the time it took for the water level to drop from one mark to the next. There were a total of 11 marks.
  2. To measure the outflow rate, we intercepted the outflow using a 25 ml beaker (not shown in Figure 2), and measured the time it took to fill to the 25 ml mark.

Results

Time elapsed since last measurement Height of water in container Time to fill 25ml beaker
Δt (s) h (cm) tf (s)
0.0 11 6.4
45.5 10 7.1
52.3 9 8.5
43.1 8 7.1
56.1 7 7.8
60.6 6 8.3
57.6 5 8.9
65.0 4 9.9
72.8 3 11.2
76.7 2 13.0
91.1 1 14.8
122.0 0 17.3

Table 1: Outflow rate, water height change with time.

To analyze the data, we calculated the total time (the cumulative sum of the elapsed time since the previous measurement) and the outflow rate. The outflow rate is the change in volume with time:

 \text{outflow rate} = \frac{volume}{time} = \frac{dV}{dt} = \frac{25 ml}{t_f}

So our data table becomes:

Time Height of water in container Outflow rate
t (s) h (cm) dV/dt (ml/s)
0.0 11 3.91
45.5 10 3.52
97.8 9 2.94
140.9 8 3.52
197 7 3.21
257 6 3.01
315.1 5 2.81
380.1 4 2.53
452.9 3 2.23
529.6 2 1.92
620.7 1 1.69
742.7 0 1.45

Table 2. Height of water and outflow rate of the bottle.

The graph of the height of the water with time shows a curve, although it is difficult to determine precisely what type of curve. My students started by trying to fit a quadratic equation to it, which should work as we’ll see in a minute.

Figure 3. The decrease in height with time is not a straight line (is non-linear).

The plot of the outflow rate versus time, however, shows a pretty good linear trend. (Note that we do not use the first three datapoints (in Table 2), which we believe are erroneous because we were still sorting out the measuring method.)

Although I’ll note here that the data should ideally be modeled using a square root equation (Torricelli’s Law), that is beyond the present scope of the problem (we’ll try that tomorrow as a follow exercise).

Plotting the data in Excel we could add a linear trendline.

Figure 4. The outflow rate decreases linearly with time.

For the linear trendline, Excel gives the equation:

 y = -0.0035 x + 3.9113

The y-axis is outflow rate (the change in volume with time), and the x axis is time (t), so the linear equation becomes:

 \frac{dV}{dt} = -0.0035 t + 3.9113

Notice that this is a differential equation.

To determine the rate of at which the height of water in the container is changing, we need to recognize that the container is cylindrical in shape, and the volume (V) of a cylinder depends on its radius (r) and height (h):

 V = \pi r^2 h

which can be substituted into differential in the rate equation:

 \frac{d(\pi r^2 h)}{dt} = -0.0035 t + 3.9113

since π and the radius (r) are constants (since the jug’s shape is a cylinder), they can be pulled out of the differential:

 \pi r^2 \frac{dh}{dt} = -0.0035 t + 3.9113

Dividing through by πr2 solves for the rate of change of height:

  \frac{dh}{dt} = \frac{-0.0035 t + 3.9113}{\pi r^2}

Isolating the coefficients gives:

  \frac{dh}{dt} = \frac{-0.0035}{\pi r^2} t + \frac{3.9113}{\pi r^2}

This equation should give the rate at which the height changes with time, however, if you look at it carefully you’ll realize that for the time range we’re using (less than 800 seconds) the value of dh/dt will always be positive. We correct this by recognizing that the outflow rate is a loss of water, so a positive outflow should result in a negative change in height, therefore we rewrite the equation as:

  -\frac{dh}{dt} = \frac{-0.0035}{\pi r^2} t + \frac{3.9113}{\pi r^2}

which gives:

  \frac{dh}{dt} = \frac{0.0035}{\pi r^2} t - \frac{3.9113}{\pi r^2}

Now comes the calculus

Now, we can use this rate equation to find the equation for height versus time by integrating with respect to time:

  \int \frac{dh}{dt} dt = \int \left( \frac{0.0035}{\pi r^2} t - \frac{3.9113}{\pi r^2} \right) dt

to get:

 h = \frac{0.0035}{2 \pi r^2} t^2 - \frac{3.9113}{\pi r^2} t + c

And all we have to do to the find the constant of integration is substitute in a known point, an initial value. As is often the case, the best point to use is the starting point where t=0 makes the rest of the calculations easier. In our case, when t=0, h=11:

 11 = \frac{0.0035}{2 \pi r^2} (0)^2 - \frac{3.9113}{\pi r^2} (0) + c

so:

 c = 11

And our final equation becomes:

 h = \frac{0.0035}{2 \pi r^2} t^2 - \frac{3.9113}{\pi r^2} t + 11

which is a quadratic equation as my students guessed before we did the calculus.

Does it work?

You will notice that in the math above, we never use the height data in determining equation for height versus time; all the calculations are based on the trendline for the outflow rate versus time.

As a result, we can compare the results of our equation to the actual measurements to see if our calculations are even close. Remarkably, they are.

Figure 5. The results from our equation (modeled) match the measured heights so well, the data points are difficult to distinguish on the graph.

Discussion

Despite all the potential for error, particularly, our relatively crude measurement techniques, and the imperfect cylindrical shape of our plastic bottle, the experiment went remarkably well.

Students found it quite challenging, and required some assistance even though this is a problem they have seen before in their textbook.

The problems in the textbook use Torricelli’s Law, which should much better describe the draining of a tank than the linear equation we find for dV/dt.

Torricelli’s Law:

 \frac{dV}{dt} = a \sqrt{2gh}

where a is the area of the outlet hole, and g is the acceleration due to gravity.

In our actual experiment it is difficult to tell that a square root function would work better. Excel does not have an option for matching a square root function, so the calculations would become more involved (although it could be set up using Excel’s iterative solver or Goal Seek function).

Conclusion

Our experiment to use calculus to determine the rate of change of the height of water in a leaking plastic water bottle was a successful exercise even though the roughness of the data collection did not permit identification of the square root law for leakage.

Euler’s Method for Approximating a Solution to a Differential Equation

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).

[inline]

  • Starting point: (x,y) = (   ,   )
  • Step size:
  • Direction:


Your browser does not support the canvas element.

  • Slope equation: dy/dx =   x +   
  • Show analytical solution:

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:

[script type=”text/javascript”]

var width=500;
var height=500;
var xrange=10;
var yrange=10;

mx = width/(2.0*xrange);
bx = width/2.0;
my = -height/(2.0*yrange);
by = height/2.0;

function upause () {
ct = 1;
}

function draw_9110(ctx, polys) {
t_9110=t_9110+dt_9110;
//ctx.fillText (“t=”+t, xp(5), yp(5));
ctx.clearRect(0,0,width,height);

polys[0].drawAxes(ctx);
polys[0].slopeField(ctx, 1.0, 1.0, 0.5);
ctx.lineWidth=2;
polys[0].draw(ctx);
polys[0].write_eqn(ctx, “slope = dy/dx”);

eu = eu_9110;
//override nsteps values
//max_nsteps_9110 = eu.nsteps;
//eu.nsteps = nsteps_9110;

if (nsteps_9110 <= max_nsteps_9110) { nsteps_9110 = nsteps_9110 + 1;} else { nsteps_9110 = 0; } eu.nsteps = nsteps_9110; if (show_parabola_9110 == 1) { polys[0].Integrate(ctx, eu.x1, eu.y1); polys[0].integral.color = '#3c3'; polys[0].integral.draw(ctx); } pos = polys[0].EulerApprox(ctx, eu); //document.getElementById('comment_9110').innerHTML = eu.nsteps+ " " + eu.x1 + " " + eu.y1+ " " + eu.dx+ " " + eu.dir+ " " + max_nsteps_9110; //starting point label ctx.textAlign="right"; ctx.textBaseline="alphabetic"; ctx.font = "14pt Calibri"; ctx.fillStyle = "#00f"; ctx.fillText ('start = ('+eu.x1+','+eu.y1+")", xp(eu.x1-0.25), yp(eu.y1-0.5)); //pos = polys[0].EulerApprox(ctx, st_pt_x_9110, st_pt_y_9110, 0.5,20,-1); //pos = polys[0].EulerApprox(ctx, st_pt_x_9110, st_pt_y_9110, 1, 6,1); //draw starting point ctx.strokeStyle="#00f"; ctx.lineWidth=2; ctx.beginPath(); ctx.arc(xp(eu.x1),yp(eu.y1),dxp(0.1),0,Math.PI*2); ctx.closePath(); ctx.stroke(); //ctx.fillText (' c='+move_dir_9110+' '+polys[1].c, xp(5), yp(5)); } //init_mouse(); var c_9110=document.getElementById("myCanvas_9110"); var ctx_9110=c_9110.getContext("2d"); var change = 0.0001; function create_lines_9110 () { //draw line var polys = []; polys.push(addPoly(0, 0.5, -1)); return polys; } function set_max_nsteps_9110 (eu) { max_nsteps_9110 = 1.5*(xrange - eu.x1)/eu.dx; } var polys_9110 = create_lines_9110(); var x1=xp(-10); var y1=yp(1); var x2=xp(10); var y2=yp(1); var dc_9110=0.05; var t_9110 = 0; var dt_9110 = 500; var st_dx_9110 = 1.0; var st_nsteps_9110 = 0; var nsteps_9110 = 0; var max_nsteps_9110 = 20; var dir_9110 = 1; var show_parabola_9110 = 1; var st_pt_x_9110 = -2; var st_pt_y_9110 = -3; var move_dir_9110 = 1.0; // 1 for up //Form interaction function update_form_9110 (x, y, dx, nsteps, dir, b, c, show_para) { x_9110.value = st_pt_x_9110+""; y_9110.value = st_pt_y_9110+""; dx_9110.value = st_dx_9110.toPrecision(2); //c_max_nsteps_9110.value = max_nsteps_9110+""; //document.getElementById('comment_9110').innerHTML = c_dir_9110; if (dir == 1) { c_dir_9110.selectedIndex = 0; } else {c_dir_9110.selectedIndex = 1;} b_9110.value = b+""; c_9110.value = c+""; if (show_para == 1) { ctrl_show.checked = 1; } else {ctrl_show.checked = 0;} } function get_euler_form_9110 () { pos = {x : parseFloat(x_9110.value), y : parseFloat(y_9110.value)} if (c_dir_9110.selectedIndex == 0) { dir = 1;} else {dir = -1;} eu = addEuler(pos.x, pos.y, parseFloat(dx_9110.value), max_nsteps_9110, dir); return eu } function update_solution_9110 (ctx, eu) { //update equations and solution document.getElementById('slope_eqn_9110').innerHTML = polys_9110[0].get_eqn("dy/dx"); //document.getElementById('eqn_9110').innerHTML = "hey"+eu.y1; polys_9110[0].Integrate(ctx, eu.x1, eu.y1); document.getElementById('eqn_9110').innerHTML = polys_9110[0].integral.get_eqn(); document.getElementById('init_pt_9110').innerHTML = "("+eu.x1+","+eu.y1+")"; } function add_euler_line() { eu_9110 = get_euler_form_9110 (); max_nsteps_9110 = eu_9110.nsteps; nsteps_9110 = 0; } var x_9110 = document.getElementById("x_9110"); var y_9110 = document.getElementById("y_9110"); var dx_9110 = document.getElementById("dx_9110"); //var c_max_nsteps_9110 = document.getElementById("nsteps_9110"); var c_dir_9110 = document.getElementById("dir_9110"); var b_9110 = document.getElementById("b_9110"); var c_9110 = document.getElementById("c_9110"); var ctrl_show = document.getElementById("chk_show_para_9110"); update_form_9110(x_9110, y_9110, dx_9110, nsteps_9110, dir_9110, polys_9110[0].b, polys_9110[0].c, show_parabola_9110); var eu_9110 = get_euler_form_9110 (); set_max_nsteps_9110(eu_9110); update_solution_9110(ctx_9110, eu_9110); setInterval("draw_9110(ctx_9110, polys_9110)", dt_9110); x_9110.onchange = function() { eu_9110 = get_euler_form_9110 (); max_nsteps_9110 = eu_9110.nsteps; nsteps_9110 = 0; } y_9110.onchange = function() { eu_9110 = get_euler_form_9110 (); max_nsteps_9110 = eu_9110.nsteps; nsteps_9110 = 0; } c_dir_9110.onchange = function() { eu_9110 = get_euler_form_9110 (); max_nsteps_9110 = eu_9110.nsteps; nsteps_9110 = 0; } dx_9110.onchange = function() { eu_9110 = get_euler_form_9110 (); max_nsteps_9110 = eu_9110.nsteps; //max_nsteps_9110 = 1.4 * (10 - eu_9110.x1)/eu_9110.dx; set_max_nsteps_9110(eu_9110); nsteps_9110 = 0; } b_9110.onchange = function() { polys_9110[0].set_b(parseFloat(this.value)); eu_9110 = get_euler_form_9110 (); max_nsteps_9110 = eu_9110.nsteps; nsteps_9110 = 0; } c_9110.onchange = function() { polys_9110[0].set_c(parseFloat(this.value)); eu_9110 = get_euler_form_9110 (); max_nsteps_9110 = eu_9110.nsteps; nsteps_9110 = 0; } ctrl_show.onchange = function() { if (this.checked == 1) { show_parabola_9110 = 1;} else {show_parabola_9110 = 0} } [/script] [/inline]

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.

Slope Fields

[inline]

Your browser does not support the canvas element.

[script type=”text/javascript”]
var width=500;
var height=500;
var xrange=10;
var yrange=10;

mx = width/(2.0*xrange);
bx = width/2.0;
my = -height/(2.0*yrange);
by = height/2.0;

function draw9083(ctx, polys) {
t=t+dt;
//ctx.fillText (“t=”+t, xp(5), yp(5));
ctx.clearRect(0,0,width,height);

polys[0].drawAxes(ctx);
polys[0].slopeField(ctx, 1.0, 1.0, 0.5);
ctx.lineWidth=2;
polys[0].draw(ctx);
polys[0].write_eqn(ctx, “dy/dx “);

}

//init_mouse();

var c9083=document.getElementById(“myCanvas9083”);
var ctx9083=c9083.getContext(“2d”);

var change = 0.0001;

function create_lines9083 () {
//draw line
//document.write(“hello world! “);
var polys = [];
polys.push(addPoly(0, 0.5, 1));

return polys;
}

var polys9083 = create_lines9083();

var x1=xp(-10);
var y1=yp(1);
var x2=xp(10);
var y2=yp(1);
var dy=0.01;

var t = 0;
var dt = 50;
end_ct = 0;

var move_dir = 1; // 1 for up

//draw9083();
//document.write(“x”+x2+”x”);
//ctx9083.fillText (“n=”, xp(5), yp(5));
setInterval(“draw9083(ctx9083, polys9083)”, dt);

[/script]
[/inline]

Say you have the equation that gives you the slope of a curve (let  \frac{dy}{dx} ) be the slope):
! \frac{dy}{dx} = \frac{1}{2} x + 1

When you use integration to solve the equation, there are quite the number of possible solutions (infinite actually), because when you integrate:

! y = \int (\frac{1}{2} x + 1 ) dx

you get:
! y = \frac{1}{4} x^2 + x + c

where c is a constant. Unfortunately, you don’t know what c is without more information; it could be anything.

However, even without integrating, we can get a feel for what the curve will look like by plotting what the slope will look like at a bunch of different points in space. This comes in really handy when you end up with a equation for slope that is really hard — or even impossible — to solve.

The graph below show a curve of possible solutions to the slope equation. You should be able to see, as the graph slowly moves up and down, how the slope of the graph corresponds to the slope field.

[inline]

Your browser does not support the canvas element.

[script type=”text/javascript”]
var width=500;
var height=500;
var xrange=10;
var yrange=10;

mx = width/(2.0*xrange);
bx = width/2.0;
my = -height/(2.0*yrange);
by = height/2.0;

function draw9083b(ctx, polys) {
t9083b=t9083b+dt9083b;
//ctx.fillText (“t=”+t, xp(5), yp(5));
ctx.clearRect(0,0,width,height);

polys[0].drawAxes(ctx);
polys[0].slopeField(ctx, 1.0, 1.0, 0.5);
ctx.lineWidth=2;
polys[0].draw(ctx);
polys[0].write_eqn(ctx, “dy/dx “);

if (polys[1].c > yrange) {
move_dir9083b = -1.0;
polys[1].c = yrange;
}
else if (polys[1].c < -yrange) { move_dir9083b = 1.0; polys[1].c = -yrange; } polys[1].c = polys[1].c + dc9083b*move_dir9083b; polys[1].draw(ctx); polys[1].write_eqn(ctx); //ctx.fillText (' c='+move_dir9083b+' '+polys[1].c, xp(5), yp(5)); } //init_mouse(); var c9083b=document.getElementById("myCanvas9083b"); var ctx9083b=c9083b.getContext("2d"); var change = 0.0001; function create_lines9083b () { //draw line //document.write("hello world! "); var polys = []; polys.push(addPoly(0, 0.5, 1)); polys.push(addPoly(0.25, 1, 0)); polys[1].color = '#8C8'; return polys; } var polys9083b = create_lines9083b(); var x1=xp(-10); var y1=yp(1); var x2=xp(10); var y2=yp(1); var dc9083b=0.05; var t9083b = 0; var dt9083b = 100; //end_ct = 0; var move_dir9083b = 1.0; // 1 for up //draw9083b(); //document.write("x"+x2+"x"); //ctx9083b.fillText ("n=", xp(5), yp(5)); setInterval("draw9083b(ctx9083b, polys9083b)", dt9083b); [/script] [/inline]

Numerical Integration

This is an attempt to illustrate numerical integration by animating an HTML5’s canvas.

We’re trying to find the area between x = 1 and x = 5, beneath the parabola:

 y = -\frac{1}{4} x^2 + x + 4

By integrating, the area under the curve can be calculated as being 17 ⅔ (see below for the analytical solution). For numerical integration, however, the area under the curve is filled with trapezoids and the total area is calculated from the sum of all the areas. As you increase the number of trapezoids, the approximation becomes more accurate. The reduction in the error can be seen on the graph: with 1 trapezoid there is a large gap between the shaded area and the curve; more trapezoids fill in the gap better and better.

The table below show how the error (defined as the difference between the calculation using trapezoids and the analytic solution) gets smaller with increasing numbers of trapezoids (n).

Number of trapezoids Area (units2) Error (difference from 17.66)
1 15.00 2.66
2 17.00 0.66
3 17.37 0.29
4 17.50 0.16
5 17.56 0.10
6 17.59 0.07
7 17.61 0.05
8 17.63 0.03
9 17.63 0.03
10 17.64 0.02

Analytic solution

The area under the curve, between x = 1 and x = 5 can be figured out analytically by integrating between these limits.

 Area = \int_{_1}^{^5} \left(-\frac{x^2}{4}  + x + 4 \right) \,dx

 Area = \left[-\frac{x^3}{12}  + \frac{x^2}{2} + 4x \right]_1^5

 Area = \left[-\frac{(5)^3}{12}  + \frac{(5)^2}{2} + 4(5) \right] - \left[-\frac{(1)^3}{12}  + \frac{(1)^2}{2} + 4(1) \right]

 Area = \left[-\frac{125}{12}  + \frac{25}{2} + 16 \right] - \left[-\frac{1}{12}  + \frac{1}{2} + 4 \right]

 Area = \left[ 22 \frac{1}{12} \right]  - \left[4 \frac{5}{12} \right]

 Area = 17 \frac{2}{3}  = 17.6\bar{6}

(See also WolframAlpha’s solution).

Update

Demonstration of this numerical integration using four trapezoids in embeddable graphs.

A Study in Linear Equations

The effects of changing the constants in the equation of a line (y=mx+b). Image by Tess R.

My high school pre-Calculus class is studying the subject using a graphical approach. Since we’re half-way through the year I thought it would be useful to introduce some programming by building their own graphical calculators using Vpython.

Now, they all have graphical calculators, and Vpython does have its own graphing capabilities, but they’re fairly simple, only 2-dimensional, and way too automatic, so I prefer to have students program the calculators in full 3d space.

My approach to graphing is fairly simple too, but its nice because it introduces:

  • Co-ordinates: Primarily in 2d (co-ordinate plane), but 3d is easy;
  • Lists: in this case its a list of coordinates on a line;
  • Loops: (specifically for loops) to repeat actions and produce a sequence of numbers (with range); and
  • A sideways glance at matrix-like operations with arrays: A list of numbers can be treated like a matrix in some relatively simple circumstances. However, it’s not real matrix operations: multiplying a scalar by a list works like real matrix multiplication, but multiplying two lists multiplies the corresponding elements in the list.

A Simple Graphing Program

Start the program with the standard vpython header:

from visual import *

x and y axes: curves and lists

Next create the x and y axes. This introduces the curve object and lists, because Vpython draws its curve from a series of points held in a list.

To keep things simple, we’re letting the graph go from -10 to positive 10 along both axes, which makes the x-axis a line segment with only two points:

line_segment = [(-10,0), (10,0)]

The square brackets say that what’s inside as a list. In this case it’s a list of two coordinate pairs, (-10,0) and (10,0).

Now we create a line using Vpython’s curve and tell it that the positions of the points on the curve are the ones we just defined:

xaxis = curve(pos=line_segment)

To create the y-axis, we do the same thing but change the coordinate pair to (0,-10) and (0,10).

line_segment = [(0,-10),(0,10)]
yaxis = curve(pos=line_segment)

Which should produce:

Very simple x and y axes.

Tic-marks: loops

In order to be better able to keep track of things, we’ll need some tic-marks on the axes. Ideally we’d like to label them too, but I think it works well enough to save that for later.

I start by having students create the first few tic-marks and then look for the emerging pattern. Their first attempts usually look something like this:

mark1 = curve(pos=[(-10,0.3),(-10,-0.3)])
mark2 = curve(pos=[(-9,0.3),(-9,-0.3)])
mark3 = curve(pos=[(-8,0.3),(-8,-0.3)])
mark4 = curve(pos=[(-7,0.3),(-7,-0.3)])
A few tic marks.

However, instead of tediously writing out these lines we can automate it by noticing that the only things that change are the x-coordinate of the coordinate pairs: they go from -10, to -9, to -8 etc.

So we want to produce a set of numbers that go from -10 to 10, in increments of 1, and use those number to make the tic-marks. The range function will do just that: specifically, range(-10,10,1). Actually, this list only goes up to 9, but that’s okay for now.

We tell the program to go through each item in the list and give its value to the variable i using a for loop:

for i in range(-10,10,1):
    mark = curve(pos=[(i,0.3),(i,-0.3)])

In python, everything indented after the for statement is within the loop.

Tic marks on the x-axis.

The y-axis’ tic-marks are similar, and its a nice little challenge for students to figure them out. They usually come up with a separate loop, eventually, that looks something like:

for i in range(-10,10,1):
    mark = curve(pos=[(0.3, i),(-0.3,i)])
Our axes.

The Curve

Now to create a line we really only need two points. However, so that we can make other types of curves later on we’ll create a line with a series of points. We’ll create the x and y values separately:

First we set up the set of x values:

line = curve(x=arange(-10,10,0.1))

Note that I use the arange function which is just like the range function but gives you decimal values (so you can do fractions) instead of just integers.

Next we set the y values that go with the x values for the equation (in this example):
! y = 0.5 x + 2

line.y = 0.5 * line.x + 2

Finally, to make it look better, we change the color of the line to yellow:

line.color = color.yellow

In Summary

The final code looks like:

from visual import *


line_segment = [(-10,0),(10,0)]
xaxis = curve(pos=line_segment)

line_segment = [(0,-10),(0,10)]
yaxis = curve(pos=line_segment)

mark1 = curve(pos=[(-10,0.3),(-10,-0.3)])
mark2 = curve(pos=[(-9,0.3),(-9,-0.3)])
mark3 = curve(pos=[(-8,0.3),(-8,-0.3)])
mark4 = curve(pos=[(-7,0.3),(-7,-0.3)])

for i in range(-10,10,1):
    mark = curve(pos=[(i,0.3),(i,-0.3)])

for i in range(-10,10,1):
    mark = curve(pos=[(0.3, i),(-0.3,i)])


line = curve(x=arange(-10,10,0.1))
line.y = 0.5 * line.x + 2
line.color = color.yellow

which produces:

A first line: y=0.5x+2

Note on lists, arrays and matrices: You’ll notice that we create the curve, give it a list of x values (using arange), and then calculate the corresponding y values using matrix multiplication: 0.5 * line.x. This works because line.x actually stores the values as an an array, not as a list. The key difference between lists and arrays, as far as we’re concerned, is that we can get away with this type of multiplication with an array and not a list. However, an array is not a matrix, as is clearly demonstrated by the second part of the command where 2 is added to the result of the multiplication. In this case, 2 is added to each value in the array; if it were an actual matrix you need to add another matrix of the same shape that’s filled with 2’s. Right now, this is invisible to the students. The line of code makes sense. The concern is that when they do start working with matrices there might be some confusion. So watch out.

And to make any other function you just need to adjust the final line. So a parabola:
! y = x^2
would be:

line.y = line.x**2

(The two stars “**” indicates an exponent).

An Assignment

So, to assess learning, and to review the different functions we’ve learned, I asked students to produce “studies” of the different curves by demonstrating what happens when you change the different constants and coefficients in the equation.

For a straight line the general equation is:
! y = mx + b

you what happens when:

  • m > 1;
  • 0 < m < 1;
  • m < 0

and:

  • b > 1;
  • 0 < b < 1;
  • b < 0

The result is, after you add some labels, looks something like the image at the very top of this post.

This type of exercise can be done for polynomials, exponential, trigonometric, and almost any other type of functions.