Bessel function crossings

The previous looked at the angles that graphs make when they cross. For example, sin(x) and cos(x) always cross with the same angle. The same holds for sin(kx) and cos(kx) since the k simply rescales the x-axis.

The post ended with wondering about functions analogous to sine and cosine, such as Bessel functions. This post will look at that question in more detail. Specifically we’ll look at the functions Jν and Yν.

Because these two Bessel functions satisfy the same second order linear homogeneous differential equation, the Strum separation theorem says that their zeros are interlaced: between each pair of consecutive zeros of Jν is exactly one zero of Yν, and between each pair of consecutive zeros of Yν there is exactly one zero of Jν.

Plotting Bessel functions J_3 and Y_3

In the following Python code, we find zeros of Jν, then look in between for places where Jν and Yν cross. Next we find the angle the two curves make at each intersection and plot the angles.

    from scipy.special import jn_zeros, jv, yv
    from scipy.optimize import bisect
    from numpy import empty, linspace, arccos
    import matplotlib.pyplot as plt
    n = 3 # bessel function order
    N = 100 # number of zeros
    z = jn_zeros(n, N) # Zeros of J_n
    crossings = empty(N-1)
    f = lambda x: jv(n,x) - yv(n,x)    
    for i in range(N-1):
        crossings[i] = bisect(f, z[i], z[i+1])
    def angle(n, x):
        # Derivatives of J_nu and Y_nu
        dj = 0.5*(jv(n-1,x) - jv(n+1,x))
        dy = 0.5*(yv(n-1,x) - yv(n+1,x))
        top = 1 + dj*dy
        bottom = ((1 + dj**2)*(1 + dy**2))**0.5
        return arccos(top/bottom)
    y = angle(n, crossings)
    plt.xlabel("Crossing number")
    plt.ylabel("Angle in radians")

This shows that the angles steadily decrease, apparently quadratically.

Angles of crossing of J_3 and Y_3

This quadratic behavior is what we should expect from the asymptotics of Jν and Yν: For large arguments they act like shifted and rescaled versions of sin(x)/√x. So if we looked at √xJν and √xYν rather than Jν and Yν we’d expect the angles to reach some positive asymptote, and they do, as shown below.

Angles of crossing of √x J_3 and √xY_3

Related posts

Orthogonal graphs

Colin Wright posted a tweet yesterday that said that the plots of cosine and tangent are orthogonal. Here’s a plot so you can see for yourself.

Plot of cosine and tangen

Jim Simons replied with a proof so short it fits in a tweet: The product of the derivatives is -sin(x)sec²(x) = -tan(x)/cos(x), which is -1 if cos(x)=tan(x).

This made me wonder whether sine and cosine are orthogonal in the sense of graphs intersecting at right angles. They are orthogonal in the sense that their product integrates to zero over the interval [0, 2π] is zero, a fact that’s important fact in Fourier analysis, but are they orthogonal in the sense of their graphs intersecting at right angles? A graph makes this look plausible:

Plot of sine and cosine

But the graph is misleading. I made the plot without specifying the aspect ratio, using the default in Mathematica. This makes the angle between the graphs look smaller. Setting the aspect ratio to 1 shows the true picture. The two curves intersect at π/4 and 5π/4, and they intersect at an angle of 2π/3, not π/2.

The product of the slopes is not -1, but it is negative and constant, so you could multiply each function by some constant to make the product of slopes -1. Said another way, you could make the curves perpendicular by adjusting he aspect ratio.

Can you do this for other functions that are orthogonal in the inner product sense? Not easily. For example, sin(2x) and sin(3x) are orthogonal in the inner product (integral) sense, but the angles of intersection are different at different places where the curves cross.

What about other functions that are like sine and cosine? I looked at the Bessel functions J1 and J2, but the angles at the intersections vary. Ditto for Chebyshev polynomials. I suppose the difference is that sine and cosine are translations of each other, whereas that’s not the case for Bessel or Chebyshev functions. But it is true for wavelets, so you could find wavelets that are orthogonal in the sense of inner products and in the sense of perpendicular graphs.

Update: See more on how Bessel functions cross in the next post.

Fascination burnout

Here a little dialog from Anathem by Neal Stephenson that I can relate to:

“… I don’t care …”

Asribalt was horrified. “But how can you not be fascinated by—”

“I am fascinated,” I insisted. “That’s the problem. I’m suffering from fascination burnout. Of all the things that are fascinating, I have to choose just one or two.”

Area and volume of Menger sponge

The Menger sponge is the fractal you get by starting with a cube, dividing each face into a 3 by 3 grid (like a Rubik’s cube) and removing the middle square of each face and everything behind it. That’s M1, the Menger sponge at the 1st stage of its construction. The next stage repeats this process on all the little cubes that make up what’s left. That’s M2. Repeat the process over and over, and in the limit you get Menger’s sponge, a fractal with zero volume and infinite area!

This business of zero volume and infinite area may sound unreal, but the steps along the way to the limit are very tangible. Here’s a video by Aaron Benzel that let’s you fly through M3, and watch what happens if you split M3 apart.

You can compute the volume and area at each stage to show that

\mathrm{volume}(M_n) = \left(\frac{20}{27} \right )^n


\mathrm{area}(M_n) = 2\left(\frac{20}{9} \right )^n + 4 \left(\frac{8}{9} \right )^n

From these equations you can see that you can make the volume as small and you’d like, and the area as large as you like, by taking n big enough. And in case that sounds a little hand-wavey, we can get more concrete. Here’s a little code to find exactly how big a value of n is big enough.

    from math import log, ceil
    def menger_volume(n):
        return (20/27.)**n
    def menger_area(n):
        return 2*(20/9.)**n + 4*(8/9.)**n
    def volume_below(v):
        if v >=1:
            return 1
            n = log(v)/log(20/27.)
            return int(ceil(n)) 
    def area_above(a):
        if a < 2:
            return 1
            n = (log(a) - log(2))/log(20/9.)
            return int(ceil(n))
    n = volume_below(0.001)
    print( n, menger_volume(n) )
    n =  area_above(1000)
    print( n, menger_area(n) )

Related posts

Regular expression for ICD-9 and ICD-10 codes

Suppose you’re searching for medical diagnosis codes in the middle of free text. One way to go about this would be to search for each of the roughly 14,000 ICD-9 codes and each of the roughly 70,000 ICD-10 codes. A simpler approach would be to use regular expressions, though that may not be as precise.

In practice regular expressions may have some false positives or false negatives. The expressions given here have only false positives. That is, no valid ICD-9 or ICD-10 codes will go unmatched, but the regular expressions may match things that are not diagnosis codes. The latter is inevitable anyway since a string of characters could coincide with a diagnosis code but not be used as a diagnosis code. For example 1234 is a valid ICD-9 code, but 1234 in a document could refer to other things, such as a street address.

ICD-9 diagnosis code format

Most ICD-9 diagnosis codes are just numbers, but they may also start with E or V.

Numeric ICD-9 codes are at least three digits. Optionally there may be a decimal followed by one of two more digits.

An E code begins with E and three digits. These may be followed by a decimal and one more digit.

A V code begins with a V followed by two digits. These may be followed by a decimal and one or two more digits.

Sometimes the decimals are left out.

Here are regular expressions that summarize the discussion above.

    N = "\d{3}\.?\d{0,2}"
    E = "E\d{3}\.?\d?"
    V = "V\d{2}\.?\d{0,2}"
    icd9_regex = "|".join([N, E, V])

Usually E and V are capitalized, but they don’t have to be, so it would be best to do a case-insensitive match.

ICD-10 diagnosis code format

ICD-10 diagnosis codes always begin with a letter (except U) followed by a digit. The third character is usually a digit, but could be an A or B [1]. After the first three characters, there may be a decimal point, and up to three more alphanumeric characters. These alphanumeric characters are never U. Sometimes the decimal is left out.

So the following regular expression would match any ICD-10 diagnosis code.


As with ICD-9 codes, the letters are usually capitalized, but not always, so it’s best to do a case-insensitive search.

Testing the regular expressions

As mentioned at the beginning, the regular expressions here may have false positives. However, they don’t let any valid codes slip by. I downloaded lists of ICD-9 and ICD-10 codes from the CDC and tested to make sure the regular expressions here matched every code.

Regular expression features used

Character ranges are supported everywhere, such as [A-TV-Z] for the letters A through T and V through Z.

Not every regular expression implementation supports \d to represent a digit. In Emacs, for example, you would have to use[0-9] instead since it doesn’t support \d.

I’ve used \.? for an optional decimal point. (The . is a special character in regular expressions, so it needs to be escaped to represent a literal period.) Some people wold write [.]? instead on the grounds that it may be more readable. (Periods are not special characters in the context of a character classes.)

I’ve used {m} for a pattern that is repeated exactly m times, and {m,n} for a pattern that is repeated between m and n times. This is supported in Perl and Python, for example, but not everywhere. You could write \d\d\d instead of \d{3} and \d?\d? instead of \d{0,2}.

Related posts

[1] The only ICD-10 codes with a non-digit in the third position are those beginning with C4A, C7A, C7B, D3A, M1A, O9A, and Z3A.

A misunderstanding of complexity

Iterating simple rules can lead to complex behavior. Many examples of this are photogenic, and so they’re good for popular articles. It’s fun to look at fractals and such. I’ve written several articles like that here, such as the post that included the image below.

But there’s something in popular articles on complexity that bothers me, and it’s the following logical fallacy:

Complex systems can arise from iterating simple rules, therefore many complex systems do arise from iterating simple rules.

This may just be a popular misunderstanding of complexity theory, but I also get the impression that some people working in complexity theory implicitly fall for the same fallacy.

What fractals, cellular automata, and such systems show is that it is possible to start with simple rules and create a complex system. It says nothing about how often this happens. It does not follow that a particular complex system does in fact arise from iterating simple rules.

There’s a variation on the fallacy that says that while complex systems may not exactly come from iterating simple rules, it’s possible to approximate complex systems this way. Even that is dubious, as I discuss in The other butterfly effect.

At best I think these popular models serve as metaphors and cautionary tales. Strange attractors and such show that systems can exhibit unexpectedly complex behavior, and that forecasts can become useless when looking ahead too far. That’s certainly true more broadly.

But I’m skeptical of saying “You have a complex system. I have a complex system. Let’s see if my complex system models your complex system.” It’s often possible to draw loose analogies between complex systems, and these may be useful. But it’s not realistic to expect quantitative guidance to come out of this, such as using the Mandelbrot set to pick stocks.

Improving on the sieve of Eratosthenes

Ancient algorithm

Eratosthenes had a good idea for finding all primes less than an upper bound N over 22 centuries ago.

Make a list of the numbers 2 to N. Circle 2, then scratch out all the larger multiples of 2 up to N. Then move on to 3. Circle it, and scratch out all the larger multiples of 3.

Every time you start over and find a number that hasn’t been scratched out before, that means it’s not divisible by any numbers smaller than itself, so it’s prime. Circle it and repeat. This algorithm is known as the sieve of Eratosthenes.

You could turn this into an algorithm for factoring every number less than N by not just scratching off composite numbers but keeping track of what numbers they were divisible by.

Not bad for an algorithm that predates Hindu-Arabic numerals by nine centuries. But as you might imagine, there have been a few improvements over the last couple millennia.

New algorithm

A paper by Helfgott published last week [1] gives a refined version of the sieve that takes less time and less space. Helfgott is not the first to improve on the sieve of Eratosthenes, but the latest.

His paper shows that it is possible to find all primes less than N in time

O(N log N)

and space

O(N1/3 (log N)2/3).

Furthermore, it is possible to factor all integers less than N in time

O(N log N)

and space

O(N1/3 (log N)5/3).

He also addresses finding all primes and factoring all integers in an interval [N – Δ, N + Δ] provided

N1/3 (log N)2/3 ≤ Δ ≤ N.

In the case of such an interval, one can find all the primes in the interval in time O(Δ log N) and space O(Δ). And one can factor all integers in the interval in time and space O(Δ log N).

Helfgott’s paper gives detailed pseudocode for all the algorithms.

Related posts

[1] Harald Andrés Helfgott. An improved sieve of Eratosthenes, Mathematics of Computation, Published online April 23, 2019.

How category theory is applied

Instead of asking whether an area of mathematics can be applied, it’s more useful to as how it can be applied.

Differential equations are directly and commonly applied. Ask yourself what laws govern the motion of some system, write down these laws as differential equations, then solve them. Statistical models are similarly direct: propose a model and feed it data.

Linear algebra is extremely useful in application, but it’s not often applied so directly. Rarely would you look at a system and immediately see a linear system. Instead you might describe a system in terms of differential equations or statistical models, then use linear algebra as a key part of the solution process.

Numerical analysis is useful because, for example, it helps you practically solve large linear systems (which help you solve differential equations, which model physical systems).

Many areas of math are useful in application, but some are applied more directly than others. It’s a mistake to say something is not applied just because it is not applied directly.

The following quote from Colin McLarty describes how some of the most abstract math, including cohomology and category theory, is applied.

Cohomology does not itself solve hard problems in topology or algebra. It clears away tangled multitudes of individually trivial problems. It puts the hard problems in clear relief and makes their solution possible. The same holds for category theory in general.

While McLarty speaks of applications to other areas of math, the same applies to applications to other areas such as software engineering.

I suspect many supposed applications of category theory are post hoc glosses, dressing up ideas in categorical language that were discovered in a more tangible setting. At the same time, the applications of category theory may be understated because the theory works behind the scenes. As I discuss here, category theory may lead to insights that are best communicated in the language of the original problem, not in the language of categories.

Related posts

Rare and strange ICD-10 codes

ICD-10 is a set of around 70,000 diagnosis codes. ICD stands for International Statistical Classification of Diseases and Related Health Problems. The verbosity of the name is foreshadowing.

Some of the ICD-10 codes are awfully specific, and bizarre.

For example,

  • V95.4: Unspecified spacecraft accident injuring occupant
  • V97.33XA: Sucked into jet engine, initial encounter
  • V97.33XD: Sucked into jet engine, subsequent encounter

As I understand it, V97.33XD refers to a subsequent encounter with a health care professional, not a subsequent encounter with a jet engine. But you have to wonder how many people who have been sucked into a jet engine survive to have one, much less two, medical visits.

There is a specific code, Y92.146, for injuries in a prison swimming pool. It seems strange to combine a medical diagnosis and its location into a single code. Is a swimming injury in a prison pool medically different than a swimming injury in a YMCA pool?

I understand that the circumstance of a diagnosis is not recorded strictly for medical reasons. But while 70,000 is an unwieldy large set of codes, it’s kinda small when it has to account for both malady and circumstance. Surely there are 70,000 circumstances alone that are more common than being in a spacecraft, for instance.

Is there a code for being at the opera? Why yes there is: Y92.253. However, there are no codes that are unique to being at a Costco, Walmart, or Jiffy Lube.

Related posts

State privacy laws to watch

US map with states highlighted

A Massachusetts court ruled this week that obtaining real-time cell phone location data requires a warrant.

Utah has passed a law that goes into effect next month that goes further. Police in Utah will need a warrant to obtain location data or to search someone’s electronic files. (Surely electronic files are the contemporary equivalent of one’s “papers” under the Fourth Amendment.)

Vermont passed the nation’s first data broker law. It requires data brokers to register with the state and to implement security measures, but as far as I have read it doesn’t put much restriction what they can do.

Texas law expands HIPAA’s notation of a “covered entity” so that it applies to basically anyone handling PHI (protected health information).

California’s CCPA law goes into effect on January 1. In some ways it’s analogous to GDPR. It will be interesting to see what the law ultimately means in practice. It’s expected that the state legislature will amend the law, and we’ll have to wait on precedents to find out in detail what the law prohibits and allows.

Update: Maine passed a bill May 30, 2019 that prohibits ISPs from selling browsing data without consent.

Related posts