And now it’s all this

I just said what I said and it was wrong
Or was taken wrong

Fourier series in Mathematica

June 8, 2026 at 11:35 AM by Dr. Drang

After my last post—the one about using Fourier series—I started thinking about how to use Mathematica to develop Fourier series.¹ I could, of course, use the Integrate function to determine the Fourier coefficients, but Mathematica has other functions that can do the job directly once you understand how they work.

Mathematica has several Fourier functions, but I’m going to stick with the ones associated with sine series, FourierSinCoefficient and FourierSinSeries. They’re meant to be easy to use, and they are, but you need to know how they’re defined.

The coefficients used in both of these functions are defined this way:

$$f_n=\frac{2}{\pi }{\int }_0^{\pi }f(x)\sin nxdx$$

This doesn’t match up exactly with the definition I used for getting the Fourier coefficients of a loading function, $q(x)$:

$$q_n=\frac{2}{L}{\int }_0^{L}q(x)\sin \frac{n\pi x}{L}dx$$

To use the Mathematica functions to get the $q_n$, we have to do a change of variable. Let

$$z=\frac{\pi x}{L}$$

so

$$x=\frac{Lz}{\pi }\mathrm{\text{and}}dx=\frac{L}{\pi }dz$$

(Here, z is just another variable name; I’m not using it to represent a complex number.)

This changes the expression for $q_n$ to

$$q_n=\frac{2}{\pi }{\int }_0^{\pi }q\left(\frac{Lz}{\pi }\right)\sin (nz)dz$$

which means we can use the Mathematica functions as long as we substitute $Lz/\pi$ in for $x$ in the expression for $q$.

Let’s give it a try on this parabolic loading function:

$$q=4q_{max}x(L-x)$$

It shouldn’t take many Fourier terms to get a good approximation of this.

Here are the Mathematica commands I used:

q = 4 qmax x (L - x)
qn = FourierSinCoefficient[q /. x -> L z/Pi, z, n]

You can see the substitution in the first argument to FourierSinCoefficient.

After simplification, the results were

$$q_n=\{\begin{array}{c}\hfill \frac{32L^2{q}_{m}ax}{n^3{\pi }^3}\mathrm{\text{for odd}}n\\ \hfill 0\mathrm{\text{for even}}n\end{array}$$

We could also get the series (through $n=7$) directly with FourierSinSeries:

qapprox2 = FourierSinSeries[q /. x -> L z/Pi, z, 7] /. z -> Pi x /L

where the inverse substitution comes at the end to put the expression back in the form we want. Here’s a screenshot from the Mathematica notebook:

You can see that the coefficients match what we got earlier.

A quick plot of the difference between this truncated Fourier series and the original parabolic function shows that, as expected, the series does a good job of replicating the original:

In this plot, the horizontal axis is the fraction of $L$ and the vertical axis is the fraction of $q_{max}$.

If you’re interested, here’s the entire Mathematica notebook:

Now that I understand the way these functions work, I can do more complicated Fourier analysis in Mathematica without questioning myself on whether I’m using them correctly.

---

Simply supported beam—Fourier series solution of the ODE

June 6, 2026 at 3:36 PM by Dr. Drang

We’re finally here, at the end of all things. In this post, we’ll use a Fourier series to get the formula for the center deflection of a simply supported beam with a uniformly distributed load. We’ll see some of the same math that we saw in the previous Fourier series solution, but the fundamental approach will be different.

Let’s start with the fourth-order ordinary differential equation for a beam with a general loading function, $q(x)$:

$$EI{y}^{iv}=q$$

The simple supports at both ends give us these boundary conditions:

$$y(0)=y(L)=y″(0)=y″(L)=0$$

Let’s work out the solution for a loading condition we haven’t considered before, a distributed load in the form of a sine wave:

$$q(x)=q_m\sin \frac{m\pi x}{L}$$

where $m$ is a positive integer and $q_m$ is some constant amplitude. Given the form of the governing ODE and the boundary conditions, it seems likely that the solution will look like this:

$$y(x)=y_m\sin \frac{m\pi x}{L}$$

It satisfies the boundary conditions, and if we plug it and the expression for q into the governing ODE, we get

$$EI{\left(\frac{m\pi }{L}\right)}^4{y}_m\sin \frac{m\pi }{L}=q_m\sin \frac{m\pi }{L}$$

which is a solution for all values of x if

$$y_m=\frac{1}{EI}{\left(\frac{L}{m\pi }\right)}^4{q}_m$$

This will work for any value of m that’s a positive integer.

Since the governing ODE is linear, if q is the sum of several sine terms, the solution will be the sum of several sine terms. Returning to our favorite problem,

we now have a path for expressing the solution as a sum of sines, as long as we can express the constant function, $q(x)=w$, as a sum of sines.

And of course we can. The Fourier sine series expression that fits a constant is

$$q(x)=\sum _{m=1}^{\infty }{q}_m\sin \frac{m\pi x}{L}=w$$

where

$$q_m=\frac{2}{L}{\int }_0^{L}w\sin \frac{m\pi x}{L}dx=\{\begin{array}{c}\hfill \frac{4w}{m\pi }\mathrm{\text{for odd}}m\\ \hfill 0\mathrm{\text{for even}}m\end{array}$$

This result comes from the orthogonality of the sine function. The Fourier coefficients are the same as those for a square wave.

Plugging this result into the expression for $y_m$, we get

$$y_m=\frac{1}{EI}{\left(\frac{L}{m\pi }\right)}^4\left(\frac{4w}{m\pi }\right)=\frac{4w{L}^4}{m^5{\pi }^{5}EI}$$

for odd m. Therefore,

$$y\left(\frac{L}{2}\right)=\frac{4w{L}^4}{{\pi }^{5}EI}\sum _{m=1,3,5,\dots }^{\infty }\frac{1}{m^5}\sin \frac{m\pi }{2}$$

You may recall from our earlier post that there’s a closed-form solution for this sum:

The last entry in this table with $x=\pi /2$ gives us

$$\begin{array}{cc}\hfill \sum _{m=1,3,5,\dots }^{\infty }\frac{1}{m^5}\sin \frac{m\pi }{2}& =\left(\frac{{\pi }^4}{96}\right)\left(\frac{\pi }{2}\right)-\left(\frac{{\pi }^2}{48}\right){\left(\frac{\pi }{2}\right)}^3+\left(\frac{\pi }{96}\right){\left(\frac{\pi }{2}\right)}^4\hfill \\ \hfill & =\frac{5{\pi }^5}{1536}\hfill \end{array}$$

Substituting the closed-form solution in for the sum yields our old friend:

$$y\left(\frac{L}{2}\right)=\left(\frac{4w{L}^4}{{\pi }^{5}EI}\right)\left(\frac{5{\pi }^5}{1536}\right)=\frac{5w{L}^4}{384EI}$$

---

The purpose of this series of posts—apart from a bit of showing off—was to demonstrate why I love structural analysis. It has, even in the simplest of problems, a depth that rewards those who study it. And when you’re confronted with problems that aren’t so simple, you can draw on that depth to understand and solve them.

---

Simply supported beam—Newmark’s method

June 6, 2026 at 7:44 AM by Dr. Drang

This would have been the last post in the series if I hadn’t realized recently that another method deserved a post. So this is the penultimate. It was described in this 1943 paper by Nathan Newmark and is known—or was known—as Newmark’s method.

I say was because this method is dead. How dead? Forty-five years ago, when I was an undergraduate taking structural analysis courses, Newmark’s method wasn’t being taught in the department that he had been the head of for decades and whose main building is named after him. That’s dead.

Which isn’t to say it wasn’t a good method. Some of my professors talked about using it when they were young, and it was a practical technique when engineers used slide rules and desk calculators, but it was clear by the 80s that its time had passed. Still, it’s desncribed in the textbook I used as an undergrad, and I did eventually learn it well after I was out of school, so here it is.

Newmark’s method is basically a way of doing double integration, so you can use it to go from loading to bending moments or from bending moments to deflection. As Newmark says in the introduction to his paper,

The numerical procedure is approximate, but leads to exact moments (or deflections) when the loading diagram (or diagram of “angle changes”) is made up of segments that are bounded by straight lines or by arcs of parabolas.

For our simply supported, uniformly loaded beam, the “angle change” diagram—which we’ve been calling a curvature diagram—is a parabola, so we can get an exact solution.

Let’s say we have a beam with any kind of support conditions and some general distributed loading:

The fundamental trick of Newmark’s method is to imagine a set of simply supported beams inserted between the load and the real beam.

The small imaginary beam segments are all of the same length, $\lambda$, and the analysis proceeds by calculating the reaction forces at the ends of these beams and then building a table of loads, shears, and moments based on these reactions (reversed in direction) acting on the large original beam.

Here’s a figure from the paper that explains how to calculate the reaction forces for parabolically distributed loading:

There are simpler formulas if the loading is distributed linearly.

You may have noticed the double curvature in the right drawing, which means that the load shown there is not parabolic. Newmark suggests using the formulas in this figure even if the loading curve is of higher order—the parabolic formulas are close enough if you divide the beam into several segments.

I’ve focused here on the process of going from loads to bending moments because I think the mechanics of that is easier to understand. But the same process applies to going from curvature to deflection. We went through that in the conjugate beam post.

In fact, Newmark does our problem in the paper. He uses the constant q as the intensity of the load and splits the beam into four segments:

This may look impenetrable, but we can simplify it by using only two segments and explaining each step. Here’s the table, starting with the moment diagram:

The row labeled Moment is the value of the bending moment at each point. The Curvature row is just the moment at each point divided by EI and with the sign reversed. Those are the easy ones.

The Angle change row is based on the curvature values and the formulas in Figure 5 above. The value at the left end is

$$\frac{L/2}{24}[7\cdot 0+6\left(-\frac{w{L}^2}{8EI}\right)-0]=-\frac{w{L}^3}{64EI}$$

where we’ve used $\lambda =L/2$. This is also the value at the right end.

The value at the center is

$$\frac{L/2}{12}[0+10\left(-\frac{w{L}^2}{8EI}\right)+0]=-\frac{5w{L}^3}{96EI}$$

We start the Average slope row with the knowledge that the slope at the center of the beam is zero, so if the angle change across the center is

$$-\frac{5w{L}^3}{96EI}$$

then the average slope must go from

$$\frac{5w{L}^3}{192EI}\mathrm{\text{to}}-\frac{5w{L}^3}{192EI}$$

as we move from the left half of the beam to the right half.

(If you’re having trouble with this, remember that angle change is analogous to shear and think of a beam with an antisymmetric shear diagram. If the change in shear across the centerline is $-F$, then the shear must be $F/2$ in the left half of the beam and $-F/2$ in the right half.)

For the Deflection row, we know that the deflection is zero at the left and right ends. Since the average slope in the left half of the beam is

$$\frac{5w{L}^2}{192EI}$$

the deflection at the center must be

$$\left(\frac{5w{L}^3}{192EI}\right)\left(\frac{L}{2}\right)=\frac{5w{L}^4}{384EI}$$

an answer we’ve now seen twelve times.

---

Simply supported beam—dummy unit load method

June 5, 2026 at 7:51 AM by Dr. Drang

We’re in the home stretch of this series, so ANIAT will soon go back to complaining about Apple’s UI choices.¹ Next week’s WWDC keynote should provide some inspiration.

But today’s post covers our eleventh method for deriving the center deflection of a uniformly loaded simply supported beam: the dummy unit load method. When I was an undergraduate, this was the method we used to determine the deflections of truss structures, but it’s more general than that. My favorite explanation of why it works is in Nicholas J. Hoff’s The Analysis of Structures. I have the original Wiley hardcover of this book, but you can apparently get both paperback and hardcover reprints.

Hoff uses the principle of virtual work in his explanation, and I’d like to quote him here, but unfortunately his explanation is split between one section on the analysis of trusses and another on the analysis of beams. I couldn’t figure out a nice way to put them together coherently, so you’re just going to have to trust me.

For a beam, you can get the deflection at a particular point by following these steps:

1. Work out an equation for the bending moment, M, in the beam for the given set of applied loads.
2. Imagine that same beam with those loads removed (this is the dummy structure) and a concentrated unit load placed at the point at which we want the deflection. By definition, the magnitude of a unit load is one. The unit load points in the direction of the deflection we want to calculate.
3. Work out an equation for the bending moment, m, in the dummy beam with the unit load.

4. The desired deflection is

   $$\delta ={\int }_0^L\frac{mM}{EI}dx$$

Let’s apply this simple principle to our problem:

The bending moment is

$$M=\frac{w}{2}(Lx-x^2)$$

where the x coordinate starts at the left end of the beam and increases to the right, as usual.

This same structure with just a unit load at the center is

Its bending moment is

$$m=\frac{1}{2}x\mathrm{\text{for}}0\le x\le \frac{L}{2}$$

and m is symmetric about the center of the beam.

Therefore,

$$\delta =y\left(\frac{L}{2}\right)={\int }_0^L\frac{mM}{EI}dx$$

and we can take advantage of symmetry to say

$$\begin{array}{cc}\hfill y\left(\frac{L}{2}\right)& =\frac{2w}{4EI}{\int }_0^{L/2}(L{x}^2-x^3)dx\hfill \\ \hfill & =\frac{w}{2EI}{[\frac{1}{3}L{x}^3-\frac{1}{4}{x}^4]}_0^{L/2}\hfill \\ \hfill & =\frac{w}{2EI}[\frac{L^4}{24}-\frac{L^4}{64}]\hfill \\ \hfill & =\frac{5w{L}^4}{384EI}\hfill \end{array}$$

You can pretty much do this problem in your head if you’re good at figuring out least common multiples. I’m not, but after writing things out, I did realize that the LCM of $24$ and $64$ is $192$.

---

Older posts

Site search

Meta

• drdrang at leancrew
• Blog archive
• RSS feed
• JSON feed
• Mastodon
• GitHub repositories

Recent posts

Credits

This work is licensed under a Creative Commons Attribution-Share Alike 3.0 Unported License.

© 2005–2026, Dr. Drang