A Toy Model for Straggling

One of the homework problems in Intermediate Physics for Medicine and Biology (Problem 31 in Chapter 16) introduces a toy model for the Bragg peak. I won’t review that entire problem, but students derive an equation for the stopping power, S, (the energy per unit distance deposited in tissue by a high energy ion) as a function of the depth below the tissue surface, x

where S₀ is the ion’s stopping power at the surface (x = 0) and R is the ion’s range. At a glance you can see how the Bragg peak arises—the denominator goes to zero at x = R so the stopping power goes to infinity. That, in fact, is why proton therapy for cancer is becoming so popular: Energy is deposited primarily at one spot well below the tissue surface where a tumor is located, with only a small dose to upstream healthy tissue. 

One topic that comes up when discussing the Bragg peak is straggling. The idea is that the range is not a single parameter. Instead, protons have a distribution of ranges. When preparing the 6th edition of Intermediate Physics for Medicine and Biology, I thought I would try to develop a toy model in a new homework problem to illustrate straggling. 

Section 16.10 

Problem 31 ½. Consider a beam of protons incident on a tissue. Assume the stopping power S for a single proton as a function of depth x below the tissue surface is

Furthermore assume that instead of all the protons having the same range R, the protons have a uniform distribution of ranges between R – δ/2 and R + δ/2, and no protons have a range outside this interval. Calculate the average stopping power by integrating S(x) over this distribution of ranges. 

This calculation is a little more challenging than I had expected. We have to consider three possibilities for x. 

x < R — δ/2

In this case, all of the protons contribute so the average stopping power is

We need to solve the integral 

First, let

With a little analysis, you can show that

So the integral becomes

This new integral I can look up in my integral table

Finally, after a bit of algebra, I get

Well, that was a lot of work and the result is not very pretty. And we are not even done yet! We still have the other two cases. 

 R — δ/2 <  x <  R + δ/2

In this case, if the range is less than x there is no contribution to the stopping power, but if the range is greater than x there is. So, we must solve the integral

I’m not going to go through all those calculations again (I’ll leave it to you, dear reader, to check). The result is 

x   >  R + δ/2

This is the easy case. None of the protons make it to x, so the stopping power is zero. 

Well, I can’t look at these functions and tell what the plot will look like. All I can do is ask Mr. Mathematica to make the plot (he’s much smarter than I am). Here’s what he said: 

The peak of the “pure” (single value for the range) curve (the red one) goes to infinity at x = R, and is zero for any x greater than R. As you begin averaging, you start getting some stopping power past the original range, out to R + δ/2. To me the most interesting thing is that for x = R – δ/2, the stopping power is larger than for the pure case. The curves all overlap for x  > R + δ/2 (of course, they are all zero), and for fairly small values x (in these cases, about x <  0.5) the curves are all nearly equal (indistinguishable in the plot). Even a small value of δ (in this case, for a spread of ranges equal to one tenth the pure range), the peak of the stopping power curve is suppressed. 

The curves for straggling that you see in most textbooks are much smoother, but that’s because I suspect they assume a smoother distribution of range values, such as a normal distribution. In this example, I wanted something simple enough to get an analytical solution, so I took a uniform distribution over a width δ . 

Will this new homework problem make it into the 6th edition? I’m not sure. It’s definitely a candidate. However, the value of toy models is that they illustrate the physical phenomenon and describe it in simple equations. I found the equations in this example to be complicated and not illuminating. There is still some value, but if you are not gaining a lot of insight from your toy model, it may not be worth doing. I’ll leave the decision of including it in the 6th edition to my new coauthor, Gene Surdutovich. After all, he’s the expert in the interaction of ions with tissue.

Knife Edge Diffraction

Intermediate Physics for Medicine and Biology doesn’t analyze diffraction. It’s mentioned a few times—in the chapters about images (Chap. 12), sound (Chap. 13), and light (Chap. 14)—but it’s never investigated in detail. In this post, I want to take a deeper dive into diffraction. In particular, we will examine a specific example that highlights many features of this effect: diffraction from a knife edge.

Assume a plane wave of light, with intensity I₀ and wavelength λ, is incident on a half-infinite opaque screen (the knife edge, shown below with the incident light coming from the bottom upwards). If optics were entirely geometrical (light traveling in straight lines with no effect of its wavelength) the screen would cast a sharp shadow. All the light for x < 0 would continue propagating upward (in the y direction) while all the light for x > 0 would be blocked. But that’s not what really happens. Instead, light diffracts off the screen, causing fringes to appear in the region x < 0, and some light entering the shadow region of x > 0.

Optics,
by Hecht and Zajac.

Knife-edge diffraction is one of the few diffraction problems that we can solve analytically. I’ll follow the solution given in Hecht and Zajac’s textbook Optics (1979), which I used in my optics class when I was an undergraduate physics major at the University of Kansas. The solution to the knife-edge problem involves Fresnel sine and cosine integrals

I’ve plotted C and S in the top panel of the figure below. Both are odd functions that approach one half as ξ approaches infinity, albeit with many oscillations along the way. There’s lots of interesting mathematics behind these integrals, like the Cornu spiral, but this post is long enough that we don’t have time for that digression. 

If we solve for the intensity distribution beyond the screen (y > 0), we get 

This is an interesting function. When I first saw this solution plotted, I noticed oscillations on the left (x < 0) but none in the shadow region on the right (x > 0). But C and S are both odd, so they oscillate on the right and left. The middle two panels in the figure above show how this happens. Taking one half minus C and one half minus S just flips the two functions and adds a constant, so the functions vary from roughly zero to one instead of minus a half to plus a half. When you square these functions, the oscillations that are nearly equal to zero get really small (a small number like one tenth, when squared, gets very small) while the oscillations that are nearly equal to one are preserved (one squared is just one). There are still some small oscillations in the shadow region (x > 0), but somehow when you add the cosine and sine parts even they go away, and you end up with the classic solution in the bottom panel, which you see in all the textbooks.

I was curious to know is how this function behaves for different wavelengths. Diffraction effects are most important for long wavelengths, but for short wavelengths you expect to recover plain old geometrical optics. Interestingly, the wavelength λ only appears as a scaling factor for x in our solution. So changing λ merely stretches or contracts the function along the x axis. The figure below shows the solution for three different wavelengths. Note that the argument of the Fresnel sine and cosine functions is dimensionless: x has units of distance, but so do the wavelength λ and the distance past the opaque screen y, and they appear as a product under a square root. Therefore, we don’t need to worry about the units of x, y, and λ. As long as we use consistent units we are fine. As the wavelength gets small, the distribution gets crowded together close to x = 0. The red curve is the geometrical optics limit (λ = 0). The case of λ = 0.1 approaches this limit in a funny way. The amplitude of the oscillations does not change, but they fall off more quickly, so they basically only exist very near the knife edge. You wouldn’t notice them unless you looked with very fine spatial resolution. The intensity does seep into the shadow regions, and this is more pronounced at large wavelengths than at small. 

The picture above was plotted for one distance, y = 1. It’s what you would get if you put a projector screen just behind the knife edge so you could observe the intensity pattern there. What happens if you move the screen farther back (increase y). Since y and λ enter our solution only as their product, changing y is much like changing λ. Below is a plot of intensity versus x at three different values of y, for a single wavelength. Making this plot was simple: I just changed the labels on the previous plot from λ to y, and from y to λ. As you get farther way (y gets larger), the distribution spreads out. But the spreading is not linear. Because of that square root, the spreading slows down (but never stops) at large y.

You can see the full pattern of intensity in the color plots below. Remember, x is horizontal and the opaque screen is on the right, y is vertical and opaque screen is at the bottom, and color indicates intensity. Yellow is the incident intensity I₀, and blue is zero intensity (darkness). The geometrical limit would be a strip of blue on the right (the shadow) and yellow on the left (the unobstructed incident wave). The case for λ = 0.01 closely approximates this. The case of a really long wavelength (λ = 100) is interesting. The light spreads out all over, giving more of a uniform distribution. For long wavelengths, the light “bends” around the opaque screen. This is why you can still hear music even if there is an obstacle between you and the band. The sound wave diffracts because of its long wavelength (especially the low pitched notes).

Diffraction is fascinating, but it’s a bit too complicated to be included in IPMB. Nevertheless, it’s a part of the physics behind visual acuity, microscope resolution, ultrasound transducers, and many other applications.

Sine and Cosine Integrals and the Delta Function

Trigger warning: This post is for mature audiences only; it may contain Fourier transforms and Dirac delta functions. 

In Chapter 11 of Intermediate Physics for Medicine and Biology, Russ Hobbie and I examine some properties of Fourier transforms. In particular, we consider three integrals of sines and cosines. After some analysis, we conclude that these integrals are related to the Dirac delta function, δ(ω – ω’), equal to infinity at ω = ω’ and zero everywhere else (it’s a strange function consisting of one really tall, thin spike).

Are these equations correct? I now believe that they’re almost right, but not entirely. I propose that instead they should be 

You’re probably thinking “what a pity, the second three equations looks more complicated than the first three.” I agree. But let me explain why I think they’re better. Hang on, it’s a long story.

Let’s go back to our definition of the Fourier transform in Eq. 11.57 of IPMB

The first thing to note is that y(t) consists of two parts. The first depends on cos(ωt), which is an even function, meaning cos(–ωt) = cos(ωt). There’s an integral over ω, implying that many different frequencies contribute to y(t), weighted by the function C(ω). But one thing we know for sure is that when you add up the contributions from all these many frequencies, the result must be an even function (the sum of even functions is an even function). The second part depends on sin(ωt), which is an odd function, sin(–ωt) = – sin(ωt). Again, when you add up all the contributions from these many frequencies weighted by S(ω), you must get an odd function. So, we can say that we’re writing y(t) as the sum of an even part, y_even(t), and an odd part, y_odd(t). In that case, we can rewrite our Fourier transform expressions as

We should be able to take our expression for y_even(t), put our expression for C(ω) into it, and then—if all works as it should—get back y_even(t). Let’s try it and see if it works. To start I’ll just rewrite the first of the four equations listed above

Now for C(ω) I’ll use the third of the four equations listed above. In that expression, there is an integral over t, but t is a dummy variable (it’s an “internal” variable; after you do the integral, the result does not depend on t), so to keep things from getting confusing we’ll call the dummy variable by another name, t'

Next switch the order of the integrals, so the integral over t' is on the outside and the integral over ω is on the inside

Ha! There, inside the bracket, is one of those integrals were’re talking about. Okay, the variables ω and t are swapped, but otherwise it’s the same. So, let’s put in our new expression for the integral

The 2π’s cancel, and a factor of one half comes out. An integral containing a delta function just picks out the value where the argument of the delta function is zero. We get

But, we know that y_even(t) is an even function, meaning y_even(–t) equals y_even(t). So finally

It works! We go “around the loop” and get back our original function.

You could perform another calculation just like this one but for y_odd(t). Stop reading and do it, to convince yourself that again you get back to where you started from, y_odd(t) = y_odd(t).

Now, you folks who are really on the ball might realize that if you had used the old delta function relationships given in IPMB (the first three equations in this post), they would also work. (Again, try it and see.) So why use my fancy new formulas? Well, if you have an integral that adds up a bunch of cos(ωt), you know you’re gonna get an even function. There’s no way it can be equal to δ(ω – ω’), because that function is neither even nor odd. So, it just doesn’t make sense to say that summing up a bunch of even functions gives you something that isn’t even. In my new formula, that sum of two delta functions is an even function. The same argument holds for the integral with sin(ωt), which must be odd.

Finally (and this is what got me started down this rabbit hole), you often see the delta function written as

Jackson even gives this equation, so it MUST be correct. (For those of who aren’t physicists, John David Jackson wrote the highly regarded graduate textbook Classical Electrodynamics, known by physics graduate students simply as “Jackson.”)

In Jackson’s equation, i is the square root of minus one. So, this representation of the delta function uses complex numbers. You won’t see it in IPMB because Russ and I avoid complex numbers (I hate them).

Let’s use the Euler formula e^iθ = cosθ + i sinθ to change the integral in Jackson’s delta function expression to

Now use a couple trig identities, cos(A–B) = cosA cosB + sinA sinB and sin(A–B) = sinA cosB –cosA sinB, to get

This is really four integrals,

Then, using the relations between these integrals and the delta function given in IPMB (the first three equations at the top of this post), you get that the sum of these integrals is equal to

which is obviously wrong; we started with 2πδ and ended up with 4πδ. Even worse, do the same calculation for δ(ω + ω') with a plus instead of a minus in front of the ω'. I’ll leave it to you to work out the details, but you’ll get zero! Again, nonsense. However, if you use the integral relations I propose above (second set of three integrals at the top of this blog), everything works just fine (try it).

Gene Surdutovich, my new coauthor for the sixth edition of IPMB, and I are still deciding how to discuss this issue in the new edition (which we are hard at work on but I doubt will be out within the next year). I don’t want to get bogged down in mathematical minutia that isn’t essential to our book’s goals, but I want our discussion to be correct. Once the sixth edition is published, you can see how we handle it.

I haven’t seen my new delta function/Fourier integral relationships in any other textbook or math handbook. This makes me nervous. Are they correct? Moreover, Intermediate Physics for Medicine and Biology does not typically contain new mathematical results. Maybe I haven’t looked hard enough to see if someone else published these equations (if you’ve seen them before, let me know where…please!). Maybe I’ll find these results in Morse and Feshbach (another one of those textbooks known to all physics graduate students) or some other mathematical tome. I need to make a trip to the Oakland University library to look through their book collection, but right now its too cold and snowy (we got about four to five inches of the white stuff in the last 48 hours).

A Toy Model of Climate Change

Introduction

A screenshot of the online book
Math for the People.

In Introductory Physics for Medicine and Biology, Russ Hobbie and I make use of toy models. Such mathematical descriptions are not intended to be accurate or realistic. Rather, they‘re simple models that capture the main idea without getting bogged down in the details. Today, I present an example of a toy model. It’s not related to medicine or biology, but instead describes climate change. I didn’t originally derive this model. Much of the analysis below comes from other sources, such as the online book Math for the People published by Mark Branson and Whitney George.

Earth Without an Atmosphere

First, consider the earth with no atmosphere. We will balance the energy coming into the earth from the sun with the energy from the earth that is radiated out into space. Our goal will be to calculate the earth’s temperature, T.

The power density (energy per unit time per unit area, in watts per square meter) emitted by the sun is called the solar constant, S. It depends on how far you are from the sun, but at the earth’s orbit S = 1360 W/m². To get the total power impinging on our planet, we must multiply S by the area subtended by the earth, which is πR², where R is the earth’s radius (R = 6.4 × 10⁶ meters). This gives SπR² = 1.8 × 10¹⁷ W, or nearly 200,000 TW (T, or tera-, means one trillion). That’s a lot of power. The total average power consumption by humanity is only about 20 TW, so there’s plenty of energy from the sun.

We often prefer to talk about the energy loss or gain per unit area of the earth’s surface. The surface area of the earth is 4πR² (the factor of four comes from the total surface area of the spherical earth, in contrast to the area subtended by the earth when viewed from the sun). The power per unit area of the earth’s surface is therefore SπR²/4πR², or S/4.

Not all of this energy is absorbed by the earth; some is reflected back into space. The albedo, a, is a dimensionless number that indicates the fraction of the sun’s energy that is reflected. The power absorbed per unit area is then (1 – a)S/4. About 30% of the sun’s energy is reflected (a = 0.3), so the power of sunlight absorbed by the earth per unit of surface area is 238 W/m².

What happens to that energy? The sun heats the earth to a temperature T. Any hot object radiates energy. Such thermal radiation is analyzed in Section 14.8 of Intermediate Physics for Medicine and Biology. The radiated power per unit area is equal to eσT⁴. The symbol σ is the Stefan-Boltzmann constant, σ = 5.7 × 10^–8 W/(m² K⁴). As stated earlier, T is the earth’s temperature. When raising the temperature to the fourth power, T must be expressed as the absolute temperature measured in kelvin (K). Sometimes it’s convenient at the end of a calculation to convert kelvin to the more familiar degrees Celsius (°C), where 0°C = 273 K. But remember, all calculations of T⁴ must use kelvin. Finally, e is the emissivity of the earth, which is a measure of how well the earth absorbs and emits radiation. The emissivity is another dimensionless number ranging between zero and one. The earth is an excellent emitter and absorber, so e = 1. From now on, I’ll not even bother including e in our equations, in which case the power density emitted is just σT⁴.

Let’s assume the earth is in steady state, meaning the temperature is not increasing or decreasing. Then the power in must equal the power out, so 

(1 – a)S/4 = σT⁴. 

Solving for the temperature gives

T = ∜[(1 – a)S/4σ] .

Because we know a, S, and σ, we can calculate the temperature. It is T = 254 K = –19°C. That’s really cold (remember, in the Celsius scale water freezes at 0°C). Without an atmosphere, the earth would be a frozen wasteland.

Earth With an Atmosphere

Often we can learn much from a toy model by adding in complications, one by one. Now, we’ll include an atmosphere around earth. We must keep track of the power into and out of both the earth and the atmosphere. The earth has temperature T_E and the atmosphere has temperature T_A.

First, let’s analyze the atmosphere. Sunlight passes right through the air without being absorbed because it’s mainly visible light and our atmosphere is transparent in the visible part of the spectrum. The main source of thermal (or infrared) radiation (for which the atmosphere is NOT transparent) is from the earth. We already know how much that is, σT_E⁴. The atmosphere only absorbs a fraction of the earth’s radiation, e_A, so the power per unit area absorbed by the atmosphere is e_AσT_E⁴.

Just like the earth, the atmosphere will heat up to a temperature T_A and emit its own thermal radiation. The emitted power per unit area is e_AσT_A⁴. However, the atmosphere has upper and lower surfaces, and we’ll assume they both emit equally well. So the total power emitted by the atmosphere per unit area is 2e_AσT_A⁴.

If we balance the power in and out of the atmosphere, we get 

e_AσT_E⁴ = 2e_AσT_A⁴. 

Interestingly, the fraction of radiation absorbed by the atmosphere, e_A, cancels out of our equation (a good emitter is also a good absorber). The Stefan-Boltzmann constant σ also cancels, and we just get T_E⁴ = 2T_A⁴. If we take the forth root of each side of the equation, we find that T_A = 0.84 T_E. The atmosphere is somewhat cooler than the earth.

Next, let’s reanalyze the power into and out of the earth when surrounded by an atmosphere. The sunlight power per unit area impinging on earth is still (1 – a)S/4. The radiation emitted by the earth is still σT_E⁴. However, the thermal radiation produced by the atmosphere that is aimed inward toward the earth is all absorbed by the earth (since the emissivity of the earth is one, e_E = 1), so this provides another factor of e_AσT_A⁴. Balancing power in and out gives

(1 – a)S/4 + e_AσT_A⁴ = σT_E⁴ .

Notice that if e_A were zero, this would be the same relationship as we found when there was no atmosphere: (1 – a)S/4 = σT_E⁴. The atmosphere provides additional heating, warming the earth.

We found earlier that T_E⁴ = 2T_A⁴. If we rewrite this as T_A⁴ = T_E⁴/2 and plug that into our energy balance equation, we get

(1 – a)S/4 + e_AσT_E⁴/2 = σT_E⁴ .

With a bit of algebra, we find

(1 – a)S/4 = σT_E⁴ (1 – e_A/2) .

Solving for the earth’s temperature gives

T_E = ∜[(1 – a)S/4σ] ∜[1/(1 – e_A/2) ] .

If e_A were zero, this would be exactly the relationship we had for no atmosphere. The fraction of energy absorbed by the atmosphere is not zero, however, but is approximately e_A = 0.8. The atmosphere provides a dimensionless correction factor of ∜[1/(1 – e_A/2)]. The temperature we found previously, 254 K, is corrected by this factor, 1.136. We get T_E = 288.5 K = 15.5 °C. This is approximately the average temperature of the earth. Our atmosphere raised the earth’s temperature from –19°C to +15.5°C, a change of 34.5°C.

Climate Change

To understand climate change, we need to look more deeply into the meaning of the factor e_A, the fraction of energy absorbed by the atmosphere. The main constituents of the atmosphere—oxygen and nitrogen—are transparent to both visible and thermal radiation, so they don’t contribute to e_A. Thermal energy is primarily absorbed by greenhouse gases. Examples of such gases are water vapor, carbon dioxide, and methane. Methane is an excellent absorber of thermal radiation, but its concentration in the atmosphere is low. Water vapor is a good absorber, but water vapor is in equilibrium with liquid water, so it isn’t changing much. Carbon dioxide is a good absorber, has a relatively high concentration, and is being produced by burning fossil fuels, so a lot of our discussion about climate change focuses on carbon dioxide.

The key to understanding climate change is that greenhouse gasses like carbon dioxide affect the fraction of energy absorbed, e_A. Suppose an increase in the carbon dioxide concentration in the atmosphere increased e_A slightly, from 0.80 to 0.81. The correction factor  ∜(1/(1 – e_A/2) ) would increase from 1.136 to 1.139, changing the temperature from 288.5 K to 289.3 K, implying an increase in temperature of 0.8 K. Because changes in temperature are the same if expressed in kelvin or Celsius, this is a 0.8°C rise. A small change in e_A causes a significant change in the earth’s temperature. The more carbon dioxide in the atmosphere, the greater the temperature rise: Global warming.

Feedback

We have assumed the earth’s albedo, a, is a constant, but that is not strictly true. The albedo depends on how much snow and ice cover the earth. More snow and ice means more reflection, a larger albedo, a smaller amount of sunlight absorbed by the earth, and a lower temperature. But a lower temperature means more snow and ice. We have a viscous cycle: more snow and ice leads to a lower temperature which leads to more snow and ice, which leads to an even lower temperature, and so on. Intermediate Physics for Medicine and Biology dedicates an entire chapter to feedback, but it focuses mainly on negative feedback that tends to maintain a system in equilibrium. A viscous cycle is an example of positive feedback, which can lead to explosive change. An example from biology is the upstroke of a nerve action potential: an increase in the electrical voltage inside a nerve cell leads to an opening of sodium channels in the cell membrane, which lets positively charged sodium ions enter the cell, which causes the voltage inside the cell to increase even more. The earth’s climate has many such feedback loops. They are one of the reasons why climate modeling is so complicated.

Conclusion

Today I presented a simple description of the earth’s temperature and the impact of climate change. Many things were left out of this toy model. I ignored differences in temperature over the earth’s surface and within the atmosphere. I neglected ocean currents and the jet stream that move heat around the globe. I did not account for seasonal variations, or for other greenhouse gasses such as methane and water vapor, or how the amount of water vapor changes with temperature, or how clouds affect the albedo, and a myriad of other factors. Climate modeling is a complex subject. But toy models like I presented today provide insight into the underlying physical mechanisms. For that reason, they are crucial for understanding complex phenomena such as climate change.

Taylor Series

The Taylor series is particularly useful for analyzing how functions behave in limiting cases. This is essential when translating a mathematical expression into physical intuition, and I would argue that the ability to do such translations is one of the most important skills an aspiring physicist needs. Below I give a dozen examples from Intermediate Physics for Medicine and Biology, selected to give you practice with Taylor series. In each case, expand the function in the dimensionless variable that I specify. For every example—and this is crucial—interpret the result physically. Think of this blog post as providing a giant homework problem about Taylor series.

Find the Taylor series of:

1. Eq. 2.26 as a function of bt (this is Problem 26 in Chapter 2). The function is the solution for decay plus input at a constant rate. You will need to look up the Taylor series for an exponential, either in Appendix D or in your favorite math handbook. I suspect you’ll find this example easy.
2. Eq. 4.69 as a function of ξ (this is Problem 47 in Chapter 4). Again, the Taylor series for an exponential is required, but this function—which arises when analyzing drift and diffusion—is more difficult than the last one. You’ll need to use the first four terms of the Taylor expansion.
3. The argument of the inverse sine function in the equation for C(r,z) in Problem 34 of Chapter 4, as a function of z/a (assume r is less than a). This expression arises when calculating the concentration during diffusion from a circular disk. Use your Taylor expansion to show that the concentration is uniform on the disk surface (z = 0). This calculation may be difficult, as it involves two different Taylor series. 
4. Eq. 5.26 as a function of ax. Like the first problem, this one is not difficult and merely requires expanding the exponential. However, there are two equations to analyze, arising from the study of countercurrent transport. 
5. Eq. 6.10 as a function of z/c (assume c is less than b). You will need to look up or calculate the Taylor series for the inverse tangent function. This expression indicates the electric field near a rectangular sheet of charge. For z = 0 the electric field is constant, just as it is for an infinite sheet.
6. Eq. 6.75b as a function of b/a. This equation gives the length constant for a myelinated nerve axon with outer radius b and inner radius a. You will need the Taylor series for ln(1+x). The first term of your expansion should be the same as Eq. 6.75a: the length constant for an unmyelinated nerve with radius a and membrane thickness b.
7. The third displayed equation of Problem 46 in Chapter 7 as a function of t/t_C. This expression is for the strength-duration curve when exciting a neuron. Interestingly, the short-duration behavior is not the same as for the Lapicque strength-duration curve, which is the first displayed equation of Problem 46.
8. Eq. 9.5 as a function of [M']/[K]. Sometimes it is tricky to even see how to express the function in terms of the required dimensionless variable. In this case, divide both sides of Eq. 9.5, to get [K']/[K] in terms of [M']/[K]. This problem arises from analysis of Donnan equilibrium, when a membrane is permeable to potassium and chloride ions but not to large charged molecules represented by M’.
9. The expression inside the brackets in Eq. 12.42 as a function of ξ. The first thing to do is to find the Taylor expansion of sinc(ξ), which is equal to sin(ξ)/ξ. This function arises when solving tomography problems using filtered back projection.
10. Eq. 13.39 as a function of a/z. The problem is a little confusing, because you want the limit of large (not small) z, so that a/z goes to zero. The goal is to show that the intensity falls off as 1/z² for an ultrasonic wave in the Fraunhoffer zone.
11. Eq. 14.33 as a function of λk_BT/hc. This problem really is to determine how the blackbody radiation function behaves as a function of wavelength λ, for short wavelength (high energy) photons. You are showing that Planck's blackbody function does not suffer from the ultraviolet catastrophe.
12. Eq. 15.18 as a function of x. (This is Problem 15 in Chapter 15). This function describes how the Compton cross section depends on photon energy. Good luck! (You’ll need it).

Brook Taylor

Who was Taylor? Brook Taylor (1685-1731) was an English mathematician and a fellow of the Royal Society. He was a champion of Newton’s version of the calculus over Leibniz’s, and he disputed with Johann Bernoulli. He published a book on mathematics in 1715 that contained his series.

Numerical Solution of the Quadratic Formula

Homework Problem 7 in Chapter 11 of Intermediate Physics for Medicine and Biology examines fitting data to a straight line. In that problem, the four data points to be fit are (100, 4004), (101, 4017), (102, 4039), and (103, 4063). The goal is to fit the data to the line y = ax + b. For this data, one must perform the calculation of a and b to high precision or else you get large errors. The solution manual (available to instructors upon request) says that

This problem illustrates how students can run into numerical problems if they are not careful. With modern calculators that carry many significant figures, this may seem like a moot point. But the idea is still important and can creep subtly into computer computations and cause unexpected, difficult-to-debug errors.

Numerical Recipes

Are there other examples of numerical errors creeping into calculations? Yes. You can find one discussed in Numerical Recipes that involves the quadratic formula.

We all know the quadratic formula from high school. If you have a quadratic equation of the form

the solution is

For example,
 

has two solutions
 

so x = 1 or 2.

Now, suppose the coefficient b is larger,

The solution is

so x = 300 or 0.00667.

This calculation is susceptible to numerical error. For instance, suppose all numerical calculations are performed to only four significant figures. Then when you reach the step

you must subtract 8 from 90,000. You get 89992, which to four significant figures becomes 89990, which has a square root of (again to four significant figures) 300.0. The solutions are therefore x = 300 or 0. The large solution (300) is correct, but the small solution (0 instead of 0.00667) is completely wrong. The main reason is that when using the minus sign for ± you must subtract two numbers that are almost the same (in this case, 300 – 299.98667) to get a much smaller number.

You might say “so what! Who uses only four significant figures in their calculations?” Okay, try solving

where I increased b from 300 to 3000. You’ll find that using even six significant figures gives one nonsense solution (try it). As you make b larger and larger, the calculation becomes more and more difficult. The situation can cause unexpected, difficult-to-debug errors.

What’s the moral to this story? Is it simply that you must use high precision when doing calculations? No. We can do better. Notice that the solution is fine when using the plus sign in the quadratic equation. We need make no changes. It’s the negative sign that gives the problem,
 

Let’s try a trick; multiply the expression by a very special form of one:

Simplifying, we get

Voilà! The denominator has the plus sign in front of the square root, so it is not susceptible to numerical error. The numerator is simplicity itself. Try solving x² – 300x + 2 = 0 using math to four significant figures,

 

No error, even with just four sig figs. The problem is fixed!

I should note that the problem is fixed only for negative values of b. If b is positive, you can use an analogous approach to get a slightly different form of the solution (I’ll leave that as an exercise for the reader).

So, the moral of the story is: if you find that your numerical calculation is susceptible to numerical error, fix it! Look for a trick that eliminates the problem. Often you can find one.