Take a topic from your own degree, find the differential equation inside it, understand it properly, and teach it to your class in 15 minutes.
Every project below starts from something you have studied, or will study, in your own department — a machine, a structure, a process, a network, a classroom — and asks one question: which differential equation governs it, and what does its solution actually tell an engineer? You are not asked to invent new mathematics. You are asked to understand one real model well enough to explain it to your classmates without hiding behind slides.
Everyone in the group must be in your own section, because you present to your own class. Mixed-specialty groups are welcome, as long as every member can explain the whole project.
Each brief names the course in your degree it builds on, the equation, what you must learn and what you must deliver. Read three or four briefs fully before deciding. You may take a project from another specialty if your group is ready to do it justice.
One member signs in and books for the group, giving the group name and every member's university email. The moment it is booked, that project disappears for your section, so no two groups in the same class present the same topic. Each student may belong to one group only.
"What your group must learn" is the minimum, not the whole job. Work through the mathematics by hand, then check it numerically. Divide the reading between you, but make sure every member understands all of it: any member may be asked any question.
Every brief lists five deliverables: a derivation, a solved or simulated case with real numbers, a check that your answer is right, one clear figure, and an honest statement of where the model stops working.
Every member speaks. Time yourselves in advance — running over costs marks. Expect questions from your classmates and from me afterwards.
| Criterion | Weight | What I am looking for |
|---|---|---|
| Understanding the equation | 30% | Correct equation, every term explained, solution or simulation done properly |
| Connection to your field | 25% | A real engineering decision that depends on this model, explained in your own words |
| Depth of study | 20% | Evidence you went past a web search: derivation, verification, honest limits, sources |
| Presentation | 15% | Clear structure, readable figures, 15 minutes kept, questions answered |
| Whole-group ownership | 10% | Every member speaks and can answer questions about any part |
Browse by specialty, or filter the list. A project marked Taken has been booked by another group in your section.
When a fault occurs on a transmission line, protection clears it in a few cycles. Whether the generator returns to synchronism or pole-slips depends on how long the fault lasted, and that question is answered by integrating a nonlinear second-order ODE. Utilities set relay timings from exactly this calculation.
The differential equation: $$\frac{2H}{\omega_s}\frac{d^2\delta}{dt^2}+D\frac{d\delta}{dt}=P_m-\frac{E'V}{X}\sin\delta$$
\(\delta\) is the rotor angle in electrical radians, \(H\) the inertia constant in seconds, \(\omega_s\) the synchronous speed in electrical rad/s, \(D\) the damping coefficient, \(P_m\) the mechanical input power in per unit, \(E'\) the machine internal voltage behind transient reactance, \(V\) the infinite-bus voltage and \(X\) the total reactance between them. All powers are in per unit on the machine base.
If your group wants more: Add a simple first-order exciter (automatic voltage regulator) as a third state and show how it changes the critical clearing time and can introduce a low-frequency oscillatory mode.
A voltage step launched onto a cable does not appear instantly at the far end; it travels, reflects off mismatched terminations and can double at an open circuit. Signal integrity on a PCB trace and switching overvoltages on a feeder are the same mathematics. The model is a pair of coupled first-order PDEs that reduce to a wave equation.
The differential equation: $$\frac{\partial v}{\partial x}=-L\frac{\partial i}{\partial t}-Ri,\qquad \frac{\partial i}{\partial x}=-C\frac{\partial v}{\partial t}-Gv \;\Longrightarrow\; \frac{\partial^2 v}{\partial x^2}=LC\frac{\partial^2 v}{\partial t^2}+(RC+LG)\frac{\partial v}{\partial t}+RGv$$
\(v(x,t)\) and \(i(x,t)\) are the line voltage and current at position \(x\) and time \(t\); \(R\), \(L\), \(G\) and \(C\) are the series resistance, series inductance, shunt conductance and shunt capacitance per unit length. With \(R=G=0\) the second equation becomes the lossless wave equation with speed \(u=1/\sqrt{LC}\).
If your group wants more: Add frequency-dependent series resistance R(omega) proportional to sqrt(omega) to represent skin effect and show, by inverse Fourier transform of the phasor solution, how the step edge is smeared.
Every carrier recovery circuit, frequency synthesiser and clock-and-data recovery block is a PLL. The small-signal transfer function taught in class hides the question that matters in practice: given a frequency offset, will the loop lock at all, and how long will it take? That question is a nonlinear phase-plane problem.
The differential equation: $$\frac{d\phi}{dt}=\Delta\omega-K\left(x+\frac{\tau_2}{\tau_1}\sin\phi\right),\qquad \frac{dx}{dt}=\frac{\sin\phi}{\tau_1} \;\Longleftrightarrow\; \frac{d^2\phi}{dt^2}+\frac{K\tau_2}{\tau_1}\cos\phi\,\frac{d\phi}{dt}+\frac{K}{\tau_1}\sin\phi=0$$
\(\phi\) is the phase error between input and VCO in radians, \(x\) is the integrator state of the active loop filter \(F(s)=(1+\tau_2 s)/(\tau_1 s)\), \(K\) is the loop gain in rad/s (phase-detector gain times VCO gain), \(\Delta\omega\) is the constant offset between the input frequency and the free-running VCO frequency, and \(\tau_1,\tau_2\) are the loop-filter time constants.
If your group wants more: Add band-limited white noise at the phase detector and estimate the mean time to the first cycle slip by Monte Carlo, comparing the trend with the exponential dependence on loop signal-to-noise ratio.
A boost converter is a switched circuit, so it has no single linear model; state-space averaging replaces the switching with a continuous ODE in the duty ratio. The resulting control-to-output transfer function has a zero in the right half plane, which is why a boost converter's output dips before it rises and why its voltage loop must be made slow.
The differential equation: $$L\frac{di_L}{dt}=V_g-(1-d)v_C,\qquad C\frac{dv_C}{dt}=(1-d)i_L-\frac{v_C}{R},\qquad G_{vd}(s)=\frac{\hat v_C(s)}{\hat d(s)}=\frac{V_g}{(1-D)^2}\cdot\frac{1-\dfrac{sL}{(1-D)^2R}}{1+\dfrac{sL}{(1-D)^2R}+\dfrac{s^2LC}{(1-D)^2}}$$
\(i_L\) is the averaged inductor current, \(v_C\) the averaged output capacitor voltage, \(d(t)\) the duty ratio with steady-state value \(D\), \(V_g\) the input voltage, \(L\), \(C\) and \(R\) the inductance, output capacitance and load resistance; hats denote small-signal perturbations about the operating point.
If your group wants more: Derive the discontinuous-conduction-mode averaged model for the same converter and show that the right-half-plane zero moves to a much higher frequency, then check it against a switched simulation at light load.
The steady-state equivalent circuit taught for induction machines cannot predict the torque spike and current surge during direct-on-line starting, because those are transients. Replacing the three phase quantities by two dq components turns time-varying inductances into constants and gives a fifth-order nonlinear ODE system that drives every modern vector-control scheme.
The differential equation: $$v_{qs}=r_si_{qs}+\omega\lambda_{ds}+\frac{d\lambda_{qs}}{dt},\quad v_{ds}=r_si_{ds}-\omega\lambda_{qs}+\frac{d\lambda_{ds}}{dt},\quad 0=r_r'i_{qr}'+(\omega-\omega_r)\lambda_{dr}'+\frac{d\lambda_{qr}'}{dt},\quad 0=r_r'i_{dr}'-(\omega-\omega_r)\lambda_{qr}'+\frac{d\lambda_{dr}'}{dt},\quad \frac{2J}{P}\frac{d\omega_r}{dt}=T_e-T_L,\quad T_e=\frac{3P}{4}L_m\left(i_{qs}i_{dr}'-i_{ds}i_{qr}'\right)$$
Subscripts \(q\) and \(d\) are the quadrature and direct axes, \(s\) and \(r\) the stator and (referred) rotor, \(v\), \(i\) and \(\lambda\) are voltage, current and flux linkage, \(r_s\) and \(r_r'\) the resistances, \(\omega\) the reference-frame speed, \(\omega_r\) the electrical rotor speed, \(P\) the number of poles, \(J\) the inertia, \(T_e\) and \(T_L\) the electromagnetic and load torques; the flux linkages close the system through \(\lambda_{qs}=L_si_{qs}+L_mi_{qr}'\), \(\lambda_{ds}=L_si_{ds}+L_mi_{dr}'\), \(\lambda_{qr}'=L_r'i_{qr}'+L_mi_{qs}\), \(\lambda_{dr}'=L_r'i_{dr}'+L_mi_{ds}\) with \(L_s=L_{ls}+L_m\) and \(L_r'=L_{lr}'+L_m\).
If your group wants more: Implement indirect rotor-flux-oriented control on the same model and compare the starting torque and current with the direct-on-line case at equal run-up time.
Energising an unloaded transformer can draw a current many times rated for several cycles, which trips protection that was not designed to expect it. The cause is a first-order ODE in flux linkage with a strongly nonlinear magnetisation characteristic, and the peak depends on the point on the voltage wave at which the breaker closes and on the residual flux left in the core.
The differential equation: $$\frac{d\lambda}{dt}=\sqrt{2}\,V\sin(\omega t+\alpha)-R\,i(\lambda),\qquad i(\lambda)=\frac{\lambda}{L_m}\;\text{for}\;|\lambda|\le\lambda_s,\qquad i(\lambda)=\operatorname{sgn}(\lambda)\left[\frac{\lambda_s}{L_m}+\frac{|\lambda|-\lambda_s}{L_a}\right]\;\text{for}\;|\lambda|\gt\lambda_s$$
\(\lambda\) is the core flux linkage in Wb-turns, \(V\) the rms supply voltage, \(\omega\) the supply angular frequency, \(\alpha\) the closing angle on the voltage wave, \(R\) the winding plus source resistance, \(L_m\) the unsaturated magnetising inductance, \(\lambda_s\) the saturation flux linkage and \(L_a\) the much smaller air-core inductance above the knee; the initial condition \(\lambda(0)=\lambda_r\) is the residual flux.
If your group wants more: Extend to a three-phase three-limb transformer with the three closing instants separated by 120 degrees and show that the second-harmonic content of the differential current is what makes second-harmonic restraint in differential protection work.
In a microgrid with no rotating machine, inverters must share load without communicating with each other. Frequency droop achieves this, but the sharing is a dynamic equilibrium of a nonlinear system, and choosing the droop gain and the power-measurement filter badly makes the two units oscillate against each other. The stability boundary comes straight out of the eigenvalues.
The differential equation: $$\frac{d\delta_{12}}{dt}=-m_1P_1^m+m_2P_2^m,\qquad \tau\frac{dP_1^m}{dt}=-P_1^m+\tfrac12 P_L+P_c(\delta_{12}),\qquad \tau\frac{dP_2^m}{dt}=-P_2^m+\tfrac12 P_L-P_c(\delta_{12}),\qquad P_c(\delta_{12})=\frac{V_1V_2}{X_1+X_2}\sin\delta_{12}$$
\(\delta_{12}=\delta_1-\delta_2\) is the angle difference between the two inverter internal voltages, \(P_i^m\) is the low-pass filtered measured power of unit \(i\) with filter time constant \(\tau\), \(m_i\) is the frequency droop gain of unit \(i\) in rad/s per watt, \(P_L\) the total load power, \(P_c\) the circulating power, \(V_1,V_2\) the inverter voltage magnitudes and \(X_1,X_2\) the (assumed inductive) output reactances.
If your group wants more: Add a secondary integral control term that restores nominal frequency and show how its gain must be an order of magnitude slower than the droop loop to avoid a new oscillatory mode.
An MPPT algorithm is usually presented as a flowchart, which hides why it oscillates around the peak and why it can run away when irradiance changes quickly. Treating the tracker as a slow gradient law on top of the converter's own ODEs turns the tuning question into a time-scale separation question that can be answered with eigenvalues.
The differential equation: $$C_{pv}\frac{dv}{dt}=i_{pv}(v)-i_L,\qquad L\frac{di_L}{dt}=v-(1-d)V_{dc},\qquad \frac{dv_{ref}}{dt}=\gamma\left(i_{pv}(v_{ref})+v_{ref}\left.\frac{di_{pv}}{dv}\right|_{v_{ref}}\right),\qquad i_{pv}(v)=I_{ph}-I_0\left(e^{v/(nV_T)}-1\right)$$
\(v\) is the PV array terminal voltage, \(i_L\) the boost inductor current, \(d\) the duty ratio produced by an inner regulator that drives \(v\) towards \(v_{ref}\), \(V_{dc}\) the (stiff) output bus voltage, \(C_{pv}\) the input capacitance, \(L\) the boost inductance, \(\gamma\) the tracker gain, \(I_{ph}\) the light-generated current, \(I_0\) the diode saturation current, \(n\) the ideality factor and \(V_T=kT/q\) the thermal voltage; the bracket in the third equation is exactly \(dP/dv\).
If your group wants more: Add a second series-connected module with a different irradiance and a bypass diode to create two local maxima, then show numerically that the gradient tracker locks onto the wrong peak depending on its starting point.
Every state feedback controller needs states it cannot measure, so it estimates them. The estimation error itself obeys a linear ODE whose eigenvalues the designer places, and the optimal choice of gain comes from a nonlinear matrix ODE, the Riccati equation, whose steady state is what a Kalman filter actually implements. The scalar case can be solved by hand with Chapter 1 methods.
The differential equation: $$\dot e=(A-LC)e,\qquad \dot P=AP+PA^{\mathsf T}-PC^{\mathsf T}R^{-1}CP+Q,\qquad L=PC^{\mathsf T}R^{-1}$$
\(e=x-\hat x\) is the estimation error vector, \(A\) and \(C\) are the plant state and output matrices, \(L\) the observer gain, \(P\) the error covariance matrix, \(Q\) the process-noise intensity matrix and \(R\) the measurement-noise intensity; the second equation is the continuous-time matrix Riccati differential equation whose steady state gives the Kalman gain.
If your group wants more: Compare the optimal Kalman gain with a pole-placement observer whose eigenvalues you set five times faster, and show by Monte Carlo which gives the lower mean squared error under the stated noise.
A power MOSFET's on-resistance rises with junction temperature, so conduction loss rises, so the junction gets hotter still. Below a critical current this loop settles; above it the device destroys itself. The whole design rule is the sign of one coefficient in a first-order linear ODE, and it tells you the largest heatsink thermal resistance you may use.
The differential equation: $$C_{th}\frac{dT_j}{dt}=I^2R_{ds}(T_j)-\frac{T_j-T_a}{R_{th}},\qquad R_{ds}(T_j)=R_{25}\left[1+\alpha\left(T_j-25\right)\right]$$
\(T_j\) is the junction temperature in degrees Celsius, \(T_a\) the ambient temperature, \(C_{th}\) the thermal capacitance in J/K, \(R_{th}\) the junction-to-ambient thermal resistance in K/W, \(I\) the (constant) drain current, \(R_{25}\) the on-resistance at 25 C and \(\alpha\) its temperature coefficient per kelvin.
If your group wants more: Replace the single lump with a three-stage Foster network fitted to a datasheet transient thermal impedance curve, integrate it as a third-order system, and show how a pulsed current profile can be safe even though its peak exceeds the DC critical current.
The AC resistance of a conductor is larger than its DC resistance because current crowds towards the surface, which is why busbars are laminated and why transformer windings use Litz wire. The field inside the conductor satisfies a diffusion equation whose steady-state phasor form is Bessel's equation with an imaginary parameter, solved by a power series.
The differential equation: $$\frac{d^2E_z}{dr^2}+\frac1r\frac{dE_z}{dr}-j\omega\mu\sigma E_z=0,\qquad E_z(r)=E_z(a)\,\frac{J_0\!\left(j^{3/2}\sqrt{\omega\mu\sigma}\,r\right)}{J_0\!\left(j^{3/2}\sqrt{\omega\mu\sigma}\,a\right)},\qquad \delta=\sqrt{\frac{2}{\omega\mu\sigma}}$$
\(E_z(r)\) is the phasor axial electric field at radius \(r\) inside a conductor of radius \(a\), \(\omega\) the angular frequency, \(\mu\) the permeability, \(\sigma\) the conductivity, \(J_0\) the Bessel function of the first kind of order zero and \(\delta\) the skin depth; the current density is \(J_z=\sigma E_z\), and \(J_0(j^{3/2}x)=\mathrm{ber}(x)+j\,\mathrm{bei}(x)\) defines the Kelvin functions.
If your group wants more: Solve the same diffusion equation in the time domain for a step of applied surface field, using separation of variables and a Bessel-Fourier series, and show how long it takes the current to redistribute.
A digital filter running on a microcontroller is a numerical scheme applied to an analogue prototype's ODE. Which scheme you choose decides whether the filter stays stable at a given sampling rate and where its cut-off actually lands. Forward Euler, trapezoidal (bilinear) and Runge-Kutta all approximate the same ODE and give measurably different frequency responses.
The differential equation: $$\frac{d^2y}{dt^2}+2\zeta\omega_n\frac{dy}{dt}+\omega_n^2y=\omega_n^2u(t),\qquad H(s)=\frac{\omega_n^2}{s^2+2\zeta\omega_ns+\omega_n^2},\qquad s=\frac{2}{T}\frac{1-z^{-1}}{1+z^{-1}},\qquad \omega_a=\frac{2}{T}\tan\!\left(\frac{\omega_dT}{2}\right)$$
\(u(t)\) is the filter input and \(y(t)\) its output, \(\omega_n\) the undamped natural frequency in rad/s and \(\zeta\) the damping ratio; \(T\) is the sampling period, the third expression is the bilinear (trapezoidal) substitution and the fourth is the frequency pre-warping relation between the analogue frequency \(\omega_a\) and the digital frequency \(\omega_d\).
If your group wants more: Repeat the comparison at a sampling rate only three times the cut-off, add fourth-order Runge-Kutta as a third method, and quantise all coefficients to 16 bits to show which realisation degrades first.
Installing capacitors to correct power factor is routine, but the bank and the supply inductance form a parallel resonant circuit. If the resonant order lands near the 5th or 7th harmonic produced by the plant's drives, the capacitor current can multiply and the bank fails within months. The resonant order follows from a second-order ODE driven by a harmonic current source.
The differential equation: $$L_sC\frac{d^2v}{dt^2}+R_sC\frac{dv}{dt}+v=L_s\frac{di_h}{dt}+R_si_h,\qquad i_h(t)=\sum_h I_h\cos(h\omega_1t+\theta_h),\qquad h_r=\sqrt{\frac{S_{sc}}{Q_c}}$$
\(v\) is the bus voltage, \(L_s\) and \(R_s\) the supply (Thevenin) inductance and resistance, \(C\) the capacitance of the shunt bank, \(i_h\) the harmonic current injected by the nonlinear load, \(\omega_1\) the fundamental angular frequency, \(h\) the harmonic order, \(S_{sc}\) the short-circuit apparent power at the bus and \(Q_c\) the reactive power rating of the bank.
If your group wants more: Add a detuned reactor at 4.7 times the fundamental in series with the bank, recompute the impedance scan and the capacitor current, and check the design against the capacitor overvoltage and overcurrent limits in IEEE Std 18.
An electric arc is not a resistor; its conductance has its own dynamics, with a thermal time constant of microseconds to milliseconds. The same two first-order models explain why a circuit breaker sometimes fails to interrupt at current zero and why an arc furnace makes the lights of a whole district flicker. Both are nonlinear ODEs coupled to the circuit.
The differential equation: $$\frac{1}{g}\frac{dg}{dt}=\frac{1}{\tau}\left(\frac{u\,i}{P_0}-1\right)\;\text{(Mayr)},\qquad \frac{1}{g}\frac{dg}{dt}=\frac{1}{\tau}\left(\frac{u^2}{U_c^2}-1\right)\;\text{(Cassie)},\qquad L\frac{di}{dt}=v_s(t)-Ri-u,\qquad u=\frac{i}{g}$$
\(g\) is the arc conductance in siemens, \(u\) the arc voltage, \(i\) the arc current, \(\tau\) the arc thermal time constant, \(P_0\) the constant cooling power of the Mayr model, \(U_c\) the constant arc voltage of the Cassie model, and \(v_s(t)\), \(R\) and \(L\) the source voltage and the series resistance and inductance of the supply circuit.
If your group wants more: Drive the arc length with a measured or randomly generated signal, compute the modulation spectrum, and evaluate the short-term flicker severity P_st using the weighting curve of IEC 61000-4-15.
A battery management system estimates state of charge from terminal voltage, but after a current pulse the voltage keeps drifting for tens of minutes. That relaxation is lithium diffusing inside the electrode particles, which is a spherical diffusion PDE, not an RC circuit. Knowing which part of the relaxation is diffusion and which is the double layer decides how long a BMS must wait before it can trust an open-circuit voltage reading.
The differential equation: $$\frac{\partial c}{\partial t}=\frac{D_s}{r^2}\frac{\partial}{\partial r}\left(r^2\frac{\partial c}{\partial r}\right),\quad \left.\frac{\partial c}{\partial r}\right|_{r=0}=0,\quad -D_s\left.\frac{\partial c}{\partial r}\right|_{r=R_s}=\frac{I}{FA_s},\qquad \frac{dz}{dt}=-\frac{I}{Q},\qquad \frac{dv_1}{dt}=-\frac{v_1}{R_1C_1}+\frac{I}{C_1}$$
\(c(r,t)\) is the lithium concentration inside a spherical electrode particle of radius \(R_s\) with solid diffusion coefficient \(D_s\), \(I\) the cell current (positive on discharge), \(F\) the Faraday constant, \(A_s\) the total active surface area, \(z\) the state of charge with cell capacity \(Q\), and \(v_1\) the voltage across the \(R_1C_1\) pair representing charge transfer; the terminal voltage is \(v=U\!\left(c(R_s,t)\right)-v_1-IR_0\) with \(U\) the open-circuit potential and \(R_0\) the ohmic resistance.
If your group wants more: Fit D_s and R_1 C_1 to a published or measured pulse-and-rest voltage curve by least squares, and report how sensitive the fitted diffusion coefficient is to the length of the rest period used.
Every pump, turbine and turbocharger shaft carries residual unbalance, and there are shaft speeds at which the deflection it produces grows large enough to destroy the bearings. Designers must know where those critical speeds sit and how much damping is needed to run through them. The whole calculation reduces to one second-order equation written in a complex variable.
The differential equation: $$ m\ddot{r} + c\dot{r} + k r = m e \Omega^{2} e^{i\Omega t}, \qquad r(t) = x(t) + i\,y(t) $$
Here \(m\) is the disc mass, \(c\) the external viscous damping coefficient, \(k\) the lateral stiffness of the shaft at the disc, \(e\) the distance from the disc centre of mass to the shaft centre (eccentricity), \(\Omega\) the spin speed in rad/s, and \(r=x+iy\) the complex position of the shaft centre in the plane perpendicular to the shaft.
If your group wants more: Make the bearing stiffness different in the two directions. The complex equation no longer closes on r alone, because a term in the complex conjugate of r appears; show that the orbit becomes elliptical, that a backward whirl component exists, and find the two separate critical speeds.
Fins on heat sinks, engine cylinders and heat exchangers are almost always analysed at steady state, but a fin that is warming up or responding to a load change behaves differently for the first several time constants. Whether the fin can be treated as a single lumped temperature at all is decided by the Biot number. Getting this wrong produces electronics thermal designs that pass a steady-state check and then fail on a transient.
The differential equation: $$ \frac{\partial T}{\partial t} = \alpha\frac{\partial^{2} T}{\partial x^{2}} - \frac{hP}{\rho c_{p}A_{c}}\bigl(T - T_{\infty}\bigr), \qquad \alpha = \frac{k}{\rho c_{p}}, \qquad m^{2} = \frac{hP}{kA_{c}} $$
\(T(x,t)\) is the fin temperature at distance \(x\) from the base, \(\alpha\) the thermal diffusivity, \(k\) the thermal conductivity, \(\rho\) the density, \(c_p\) the specific heat, \(h\) the convection coefficient, \(P\) the fin perimeter, \(A_c\) the cross-sectional area, \(T_\infty\) the ambient temperature and \(m\) the standard fin parameter.
If your group wants more: Replace the fixed base temperature with a base temperature that oscillates sinusoidally, as for a duty-cycled processor. Find the periodic steady state, derive the thermal penetration depth, and show how far along the fin the oscillation is still detectable.
Drag on an aerofoil, a ship hull or a flat heat-exchanger plate begins with the laminar boundary layer, and the numbers quoted in every fluid mechanics textbook come from one nonlinear ODE that has no closed-form solution. Reproducing those constants yourself, rather than quoting them, is the point of this project. It is also the cleanest example in the whole curriculum of reducing a partial differential equation to an ordinary one.
The differential equation: $$ 2f'''(\eta) + f(\eta)f''(\eta) = 0, \qquad f(0)=0,\; f'(0)=0,\; \lim_{\eta\to\infty} f'(\eta) = 1, \qquad \eta = y\sqrt{\frac{U_{\infty}}{\nu x}} $$
\(\eta\) is the similarity variable combining the wall-normal coordinate \(y\) with the streamwise coordinate \(x\); \(f(\eta)\) is the dimensionless stream function, so that \(f'(\eta)=u/U_\infty\) is the normalised streamwise velocity; \(U_\infty\) is the free-stream speed and \(\nu\) the kinematic viscosity.
If your group wants more: Solve the Falkner-Skan equation, which adds a pressure-gradient parameter, for several values of that parameter covering both favourable and adverse gradients, and find numerically the value at which the wall shear reaches zero, which is the separation condition.
Brake squeal, machine-tool chatter and the creak of a sliding gate are the same phenomenon: a steady sliding motion feeds energy into an oscillation because friction falls as sliding speed rises. There is no external oscillating force, so the usual forced-vibration picture cannot explain it. The equilibrium itself is unstable and the system settles into a limit cycle, which makes this a phase-plane problem rather than a resonance problem.
The differential equation: $$ m\ddot{x} + c\dot{x} + kx = N\,\mu(v_r)\,\mathrm{sgn}(v_r), \qquad v_r = v_b - \dot{x}, \qquad \mu(v_r) = \mu_k + (\mu_s - \mu_k)e^{-|v_r|/v_0} $$
\(x\) is the displacement of the friction pad or mass, \(m\), \(c\) and \(k\) are its mass, structural damping and stiffness, \(N\) the normal contact force, \(v_b\) the belt or disc surface speed, \(v_r\) the relative sliding speed, \(\mu_s\) and \(\mu_k\) the static and kinetic friction coefficients, and \(v_0\) the speed that sets how quickly friction falls as sliding speed rises.
If your group wants more: Replace the static friction law with a rate-and-state law, in which the friction coefficient depends on an additional state variable with its own first-order ODE. Show that this gives a three-dimensional system and investigate whether the limit cycle can period-double as the belt speed is reduced.
Hydraulic actuators move excavator arms, aircraft flight controls and press machines, and their speed of response is limited not by the pump but by the compressibility of the oil trapped between the valve and the piston. That compressibility acts like a spring, giving the system a hydraulic natural frequency that caps the bandwidth of any control loop built around it. Laplace transforms turn three coupled physical relations into one transfer function that predicts this directly.
The differential equation: $$ K_q x_v - K_{ce}P_L = A\dot{x}_p + \frac{V_t}{4\beta_e}\dot{P}_L, \qquad m_t\ddot{x}_p + B_p\dot{x}_p + Kx_p = AP_L - F_L $$
\(x_v\) is the spool displacement, \(K_q\) the valve flow gain, \(K_{ce}\) the combined flow-pressure and leakage coefficient, \(P_L\) the load pressure difference across the piston, \(A\) the piston area, \(V_t\) the total trapped volume, \(\beta_e\) the effective bulk modulus of the oil, \(x_p\) the piston position, \(m_t\) the total mass referred to the piston, \(B_p\) the viscous damping, \(K\) any load spring rate and \(F_L\) the external load force.
If your group wants more: Close a proportional-plus-derivative position loop around the transfer function and use the Routh-Hurwitz criterion to find the largest proportional gain that keeps the closed-loop system stable, then show how that limit depends on the hydraulic damping ratio.
When a road embankment or a storage tank is placed on soft clay, the load is first carried by the pore water and the settlement takes years to finish. The contractor needs to know how long to wait before paving, and whether vertical drains are needed. The waiting time comes straight from the solution of a diffusion equation for excess pore pressure.
The differential equation: $$\frac{\partial u}{\partial t}=c_v\,\frac{\partial^2 u}{\partial z^2},\qquad c_v=\frac{k}{m_v\gamma_w},\qquad u(z,0)=u_0,\quad u(0,t)=u(2H,t)=0$$
u(z,t) is the excess pore water pressure at depth z and time t, c_v the coefficient of consolidation, k the permeability, m_v the coefficient of volume compressibility, gamma_w the unit weight of water, u_0 the initial excess pressure produced by the applied load, and H the longest drainage path (half the layer thickness when the layer drains top and bottom).
If your group wants more: Extend to a clay layer improved with prefabricated vertical drains by solving the axisymmetric radial consolidation equation (Barron's solution) and combining it with the vertical solution through Carrillo's formula.
A continuous footing or a railway rail does not sit on rigid supports; the soil pushes back in proportion to how far the member sinks. The resulting fourth-order equation decides how far a column load spreads along the footing and how much reinforcement is needed away from the column. Getting the characteristic length wrong gives either an over-designed raft or cracking between columns.
The differential equation: $$EI\frac{d^4 w}{dx^4}+k\,w=q(x),\qquad k=k_s b,\qquad \beta=\left(\frac{k}{4EI}\right)^{1/4}$$
w(x) is the downward deflection of the beam at position x, E the elastic modulus, I the second moment of area, k the foundation modulus per unit length obtained from the modulus of subgrade reaction k_s times the footing width b, q(x) the applied load per unit length, and beta the inverse characteristic length that sets how fast the response decays away from a load.
If your group wants more: Repeat the analysis for a finite footing of length L carrying several column loads, solve the two-point boundary value problem with free ends, and compare the design moments with those from the rigid-footing assumption used in routine design.
Base isolation protects hospitals and data centres by placing a soft layer of bearings under the building so that the structure above moves almost as a rigid body. The design questions are quantitative: how much the isolator lengthens the period, how much the roof acceleration drops, and how much displacement the bearing and the seismic gap must accommodate. All of it follows from a two-mass linear system driven by a recorded ground acceleration.
The differential equation: $$\mathbf{M}\ddot{\mathbf{u}}+\mathbf{C}\dot{\mathbf{u}}+\mathbf{K}\mathbf{u}=-\mathbf{M}\mathbf{r}\,\ddot{u}_g(t),\qquad\text{single mode:}\quad \ddot{u}+2\zeta\omega_n\dot{u}+\omega_n^{2}u=-\ddot{u}_g(t)$$
u is the vector of displacements relative to the ground (base slab and superstructure), M, C and K the mass, damping and stiffness matrices, r the influence vector of ones, u_g the ground displacement so that the doubly dotted term is the recorded ground acceleration, and in the single-mode form omega_n is the natural circular frequency and zeta the damping ratio.
If your group wants more: Replace the linear isolator by a bilinear hysteretic model with yield strength Q and post-yield stiffness, integrate the now nonlinear system, and show how the effective period and the equivalent damping change with the amplitude of shaking.
When a lane closes for maintenance a queue forms, and the tail of that queue travels backwards up the carriageway at a speed that can be predicted, sometimes faster than drivers can react. Traffic engineers use that speed to place advance warning signs and to estimate delay. The whole calculation is a first-order conservation law whose characteristics are ordinary differential equations.
The differential equation: $$\frac{\partial\rho}{\partial t}+\frac{\partial q(\rho)}{\partial x}=0,\qquad q(\rho)=\rho\,v_f\!\left(1-\frac{\rho}{\rho_j}\right),\qquad \frac{dx}{dt}=q'(\rho)=v_f\!\left(1-\frac{2\rho}{\rho_j}\right)$$
rho(x,t) is traffic density in vehicles per kilometre per lane, q the flow in vehicles per hour, v_f the free-flow speed, rho_j the jam density, and dx/dt the speed of a characteristic, along which the density stays constant.
If your group wants more: Solve the same lane closure with a triangular fundamental diagram using Newell's cumulative-count formulation, and compare the predicted total delay with the Greenshields result.
When a turbine shuts down, the water in a long headrace tunnel cannot stop instantly, so it surges up a vertical shaft and then oscillates for several minutes. If the shaft is too narrow the oscillation grows instead of dying away and the plant becomes unstable under governor control. The size of that shaft is fixed by the stability of a pair of coupled nonlinear ordinary differential equations.
The differential equation: $$\frac{L}{g}\frac{dV}{dt}=-z-cV|V|,\qquad A_s\frac{dz}{dt}=A_tV-Q_t(t)$$
V(t) is the mean water velocity in the tunnel, z(t) the surge level in the tank measured above the reservoir surface, L the tunnel length, A_t the tunnel cross-sectional area, A_s the surge tank area, g gravity, c the head-loss coefficient defined by h_f = cV|V|, and Q_t the discharge drawn by the turbine.
If your group wants more: Add a throttled orifice at the base of the tank with an asymmetric loss coefficient, or a differential (Johnson) surge tank, and show numerically how much the maximum upsurge and the required tank height are reduced.
The additive-increase multiplicative-decrease rule of TCP is normally taught as a sawtooth picture, but the aggregate behaviour of thousands of flows sharing a router is a coupled nonlinear system whose equilibrium fixes the queue length, the delay and the loss rate a user actually experiences. Router buffer sizing, active queue management and bufferbloat are all decided from this model. A group that can compute the equilibrium window and test its stability can explain why a large buffer makes latency worse without improving throughput.
The differential equation: $$\frac{dW}{dt}=\frac{1}{R(q)}-\frac{W(t)\,W(t-R)}{2R(t-R)}\,p(t-R),\qquad \frac{dq}{dt}=\frac{N\,W(t)}{R(q)}-C,\qquad R(q)=\frac{q}{C}+T_p$$
\(W\) is the average congestion window in packets, \(q\) the bottleneck queue length in packets, \(C\) the link capacity in packets per second, \(N\) the number of identical long-lived flows, \(T_p\) the fixed propagation delay, \(R(q)\) the round-trip time, and \(p\) the probability that a packet is dropped or marked.
If your group wants more: Replace the constant \(p\) by a RED marking law \(p=p(q)\), keep the round-trip delay in the equations, and find the gain at which the equilibrium loses stability through a Hopf bifurcation, producing a sustained queue oscillation.
Training any model is a discrete loop, but shrinking the step size turns gradient descent into an ODE, and momentum turns it into a damped second-order system. That view explains, with nothing more than characteristic roots, why momentum converges faster than plain gradient descent, why too much momentum makes the loss oscillate, and why an ill-conditioned problem trains slowly. Students who see this stop tuning learning rates blindly.
The differential equation: $$f(x)=\tfrac{1}{2}x^{\top}Ax-b^{\top}x,\qquad \dot{x}=-\nabla f(x)=-(Ax-b),\qquad \ddot{x}+\beta\dot{x}+\nabla f(x)=0$$
\(x(t)\in\mathbb{R}^n\) is the parameter vector at continuous time \(t\), \(f\) the loss, \(A\) a symmetric positive definite Hessian with eigenvalues \(\lambda_1\le\cdots\le\lambda_n\), \(b\) the linear term, and \(\beta>0\) the damping coefficient that corresponds to the momentum parameter.
If your group wants more: Study the Nesterov flow \(\ddot{X}+\frac{3}{t}\dot{X}+\nabla f(X)=0\) with its time-varying damping, and show numerically that its loss decays like \(O(1/t^2)\); or derive the adjoint equation \(\dot{a}=-a^{\top}\partial_x g\) used to differentiate through a neural ODE.
A processor's clock frequency is limited by heat, not by logic. Leakage current grows with temperature and leakage itself produces heat, so the lumped thermal model is a nonlinear first-order equation that can have two equilibria, one stable and one not, and can lose both in thermal runaway. Every dynamic voltage and frequency scaling controller and every thermal throttling rule is designed against this equation.
The differential equation: $$C_{th}\frac{dT}{dt}=\underbrace{\alpha C_L V^2 f}_{\text{dynamic}}+\underbrace{V I_0\exp\!\left(\frac{T-T_{\mathrm{ref}}}{T_0}\right)}_{\text{leakage}}-\frac{T-T_a}{R_{th}}$$
\(T\) is die temperature, \(C_{th}\) the thermal capacitance in J/K, \(R_{th}\) the junction-to-ambient thermal resistance in K/W, \(T_a\) ambient temperature, \(\alpha\) the activity factor, \(C_L\) the switched capacitance, \(V\) supply voltage, \(f\) clock frequency, and \(I_0,\,T_0,\,T_{\mathrm{ref}}\) the leakage reference current, its temperature scale and its reference temperature.
If your group wants more: Add a proportional-integral DVFS controller that adjusts \(f\) to hold \(T\) at a set point, form the closed-loop transfer function, and find the controller gains at which the loop becomes oscillatory when a sensor delay of 1 ms is included.
Cloth, hair and soft bodies in games are mass-spring systems whose stiff springs force an explicit integrator down to time steps far smaller than a video frame, so the simulation either explodes or runs too slowly. The standard fix in graphics is implicit Euler, and the reason it works is entirely a statement about the stability region of a numerical method in the complex plane. This project makes the connection between a graphics engine's frame budget and the eigenvalues of a stiffness matrix.
The differential equation: $$M\ddot{\mathbf{x}}+D\dot{\mathbf{x}}+K\mathbf{x}=\mathbf{f}(t)\;\Longrightarrow\;\ddot{z}_i+2\zeta_i\omega_i\dot{z}_i+\omega_i^{2}z_i=g_i(t),\qquad \omega_i^{2}=\lambda_i\!\left(M^{-1}K\right)$$
\(\mathbf{x}\) collects the displacements of all mass points, \(M\) is the diagonal mass matrix, \(K\) the stiffness matrix assembled from the springs, \(D\) the damping matrix, \(\mathbf{f}\) external forces, and \(z_i,\ \omega_i,\ \zeta_i\) are the modal coordinate, natural frequency and damping ratio of the \(i\)th mode.
If your group wants more: Implement one step of implicit Euler with a Newton solve for a nonlinear spring law, or compare against a second-order method such as the trapezoidal rule or BDF2 and show how the artificial damping changes.
Once a clock edge is faster than the time a signal takes to cross a board trace, a wire stops behaving like a node and becomes a transmission line: the edge reflects off the far end and rings, producing false switching and timing violations. The governing model is a pair of coupled PDEs whose lossless limit is the wave equation. Termination resistor values, the rule about trace length against rise time, and eye diagrams all come out of this analysis.
The differential equation: $$\frac{\partial v}{\partial x}=-L\frac{\partial i}{\partial t}-Ri,\quad \frac{\partial i}{\partial x}=-C\frac{\partial v}{\partial t}-Gv\;\Longrightarrow\;\frac{\partial^{2}v}{\partial x^{2}}=LC\frac{\partial^{2}v}{\partial t^{2}}+(RC+LG)\frac{\partial v}{\partial t}+RGv$$
\(v(x,t)\) and \(i(x,t)\) are voltage and current at position \(x\) along the trace at time \(t\), and \(R,L,G,C\) are the resistance, inductance, conductance and capacitance per unit length of the line.
If your group wants more: Add frequency-dependent loss from the skin effect, solve in the frequency domain for a pseudo-random bit sequence, and build the eye diagram to show intersymbol interference; or extend to two coupled traces and compute crosstalk.
Call centres, emergency departments and airport checkpoints all face arrival rates that change hour by hour, so the steady-state Erlang formulas taught for \(M/M/s\) do not apply directly. The correct object is an ODE for the offered load, whose solution lags the arrival rate by the service time, and a fluid equation for the queue itself. Staffing from a lagged offered load plus a square-root safety margin holds the delay probability steady across the day instead of letting it spike after each peak.
The differential equation: $$\frac{dR}{dt}=\lambda(t)-\mu R(t),\qquad \frac{dx}{dt}=\lambda(t)-\mu\min\{x(t),s(t)\},\qquad s(t)=\left\lceil R(t)+\beta\sqrt{R(t)}\,\right\rceil$$
\(\lambda(t)\) is the arrival rate at time \(t\), \(\mu\) the service rate per server, \(R(t)\) the offered load, that is the mean number in an infinite-server system, \(x(t)\) the fluid approximation to the number in the system with \(s(t)\) servers, and \(\beta\) the quality-of-service parameter that sets the safety margin.
If your group wants more: Add the diffusion refinement, that is, an ODE for the variance of the number in the system, and use it to set \(\beta\) so that the probability of waiting stays near a target such as 0.2 throughout the day.
A plant asking whether to buy a second pump, and a data centre asking whether one repair crew is enough, are asking the same question: how much of the time is the system up, given failure and repair rates. The answer is a linear system of ODEs on the states of a Markov chain, whose steady state gives availability and whose eigenvalues say how quickly the system forgets its starting condition. Contracts are written on availability numbers produced this way.
The differential equation: $$\frac{d\mathbf{p}}{dt}=\mathbf{A}\mathbf{p},\qquad \mathbf{A}=\begin{pmatrix}-2\lambda & \mu & 0\\ 2\lambda & -(\lambda+\mu) & \mu\\ 0 & \lambda & -\mu\end{pmatrix},\qquad A(t)=p_0(t)+p_1(t)$$
\(p_k(t)\) is the probability that exactly \(k\) of two identical units have failed at time \(t\), \(\lambda\) is the failure rate of one working unit, \(\mu\) the repair rate of the single repair crew, and \(A(t)\) the instantaneous availability, since the system works while at most one unit is down.
If your group wants more: Compare one repair crew against two by building the alternative generator matrix, and find the failure rate at which the second crew buys more availability than a third redundant unit; or add an imperfect repair that returns the unit to a degraded state.
Classical EOQ assumes stock sits unchanged until it is sold, which is false for food, pharmaceuticals, chemicals and volatile fuels, where a fraction of the stock is lost per unit time. Once deterioration is included, the inventory level obeys a first-order linear ODE rather than a straight line, the optimal cycle shortens, and the EOQ answer becomes an overestimate. Warehouses in a hot climate care about exactly this correction.
The differential equation: $$\frac{dI}{dt}+\theta I(t)=-D,\quad I(T)=0\;\Longrightarrow\;I(t)=\frac{D}{\theta}\left(e^{\theta(T-t)}-1\right),\quad Q=I(0)=\frac{D}{\theta}\left(e^{\theta T}-1\right)$$
\(I(t)\) is the inventory on hand at time \(t\) within a replenishment cycle of length \(T\), \(D\) the constant demand rate, \(\theta\) the deterioration rate as a fraction of stock lost per unit time, and \(Q\) the order quantity placed at the start of each cycle.
If your group wants more: Allow planned shortages with full backlogging, so a second ODE \(dI/dt=-D\) governs the shortage phase, and optimise jointly over the cycle length and the reorder point; or make demand a decreasing function of inventory age.
Equipment rarely fails at random: a bearing, a cutting tool or a corroding pipe wall accumulates damage that fails when it crosses a threshold. Modelling that damage as drift plus noise turns the probability density of wear into a PDE of advection-diffusion type, and failure becomes a first-passage problem with an absorbing boundary. The resulting life distribution then feeds directly into choosing a preventive replacement interval.
The differential equation: $$\frac{\partial p}{\partial t}=-\nu\frac{\partial p}{\partial x}+\frac{\sigma^{2}}{2}\frac{\partial^{2}p}{\partial x^{2}},\quad p(L,t)=0,\qquad f_T(t)=\frac{L}{\sigma\sqrt{2\pi t^{3}}}\exp\!\left(-\frac{(L-\nu t)^{2}}{2\sigma^{2}t}\right)$$
\(p(x,t)\) is the probability density of accumulated wear \(x\) at time \(t\), \(\nu>0\) the mean wear rate, \(\sigma\) the diffusion coefficient measuring variability of wear, \(L\) the failure threshold, and \(f_T\) the density of the first time the wear reaches \(L\).
If your group wants more: Replace the Brownian wear by a gamma process, which is monotone, obtain the life distribution numerically, and compare the optimal replacement interval and cost against the diffusion answer; or add condition-based maintenance with periodic noisy inspections.
Order variability grows as it moves upstream: a mild wobble in retail demand becomes a violent swing at the factory, driving overtime, expediting and then idle capacity. The cause is not irrational behaviour but feedback with delay, and the standard replenishment policy written as differential equations is a third-order linear system whose transfer function can amplify certain demand frequencies. Once the group has the frequency response, they can retune the smoothing constants and show the amplification fall below one.
The differential equation: $$\dot{\hat d}=\frac{d-\hat d}{T_a},\qquad \dot W=o-\frac{W}{T_p},\qquad \dot I=\frac{W}{T_p}-d,\qquad o=\hat d+\frac{I^{*}-I}{T_i}+\frac{T_p\hat d-W}{T_w}$$
\(d(t)\) is customer demand, \(\hat d(t)\) the smoothed forecast with averaging time \(T_a\), \(W(t)\) the work in progress in the pipeline with production lead time \(T_p\), \(I(t)\) the inventory with target \(I^{*}\), \(o(t)\) the order rate, and \(T_i,\ T_w\) the times over which inventory and pipeline discrepancies are corrected.
If your group wants more: Chain two or three such stages in series, retailer to distributor to factory, and show numerically how the amplification compounds; then tune \(T_i\) and \(T_w\) to keep the factory-level amplification below one and report the inventory cost of doing so.
Magnetic bearings and maglev vehicles hold a mass in the air with no contact, but an electromagnet attracting a steel ball is open-loop unstable: move the ball closer and the force grows, pulling it closer still. No amount of mechanical design removes this, only feedback. The instability appears immediately in the eigenvalues of the linearised system, which makes the rig the clearest available demonstration of why engineers linearise nonlinear models at all.
The differential equation: $$ m\ddot{x} = mg - C\left(\frac{i}{x}\right)^{2}, \qquad L\frac{di}{dt} + Ri = v(t) $$
\(x\) is the air gap between the pole face and the ball, measured downwards so that increasing \(x\) means the ball falls, \(m\) the ball mass, \(g\) gravitational acceleration, \(i\) the coil current, \(C\) the magnetic force constant of the coil and ball geometry, \(L\) and \(R\) the coil inductance and resistance, and \(v(t)\) the applied coil voltage.
If your group wants more: Include a gap-dependent inductance in the electrical equation, which introduces an extra back-EMF term coupling gap rate into the current dynamics. Rederive the Jacobian with this coupling and determine whether the extra term helps or hinders stabilisation.
Industrial robot joints use harmonic drives and gearboxes that are not rigid, so the motor encoder reading and the actual link position are two different quantities. The compliance creates a resonance that limits how hard the controller can be pushed, and an antiresonance that the motor cannot see at all. Laplace transforms show both directly in the pole-zero pattern, and explain why collocated and non-collocated feedback behave so differently.
The differential equation: $$ J_m\ddot{\theta}_m + b_m\dot{\theta}_m + k(\theta_m - \theta_\ell) = \tau_m, \qquad J_\ell\ddot{\theta}_\ell + b_\ell\dot{\theta}_\ell + k(\theta_\ell - \theta_m) = -\tau_\ell $$
\(\theta_m\) and \(\theta_\ell\) are the motor-side and link-side angles referred through the gear ratio, \(J_m\) and \(J_\ell\) their inertias, \(b_m\) and \(b_\ell\) their viscous damping coefficients, \(k\) the torsional stiffness of the joint, \(\tau_m\) the motor torque and \(\tau_\ell\) an external load torque applied at the link.
If your group wants more: Add a dead zone of plus or minus 0.2 milliradians to represent backlash, so that the transmitted torque is zero while the two angles lie within the dead zone. Simulate a motion reversal and show the limit cycle that appears once a proportional-derivative controller is tuned aggressively.
Lightweight robot arms, satellite booms and pick-and-place gantries flex while they move, and the tip keeps vibrating after the motor has stopped. Predicting where the tip actually is requires treating the link as a continuous body, which gives a partial differential equation with an infinite set of natural frequencies rather than one. The first of those frequencies usually sets the maximum useful speed of the whole machine.
The differential equation: $$ \rho A\frac{\partial^{2}w}{\partial t^{2}} + EI\frac{\partial^{4}w}{\partial x^{4}} = 0, \qquad w(0,t)=\frac{\partial w}{\partial x}(0,t)=0, \quad \frac{\partial^{2}w}{\partial x^{2}}(L,t)=\frac{\partial^{3}w}{\partial x^{3}}(L,t)=0 $$
\(w(x,t)\) is the transverse deflection at distance \(x\) along the link, \(\rho\) the density, \(A\) the cross-sectional area, \(E\) Young's modulus, \(I\) the second moment of area of the cross-section and \(L\) the link length; the boundary conditions given are clamped at the hub and free at the tip.
If your group wants more: Add a rigid tip mass with its own rotary inertia. This changes two of the boundary conditions and therefore the frequency equation; derive the new equation, solve it numerically, and plot how the first two frequencies fall as the payload rises from zero to the mass of the link itself.
Every drone, phone and self-balancing robot carries a gyroscope that drifts and an accelerometer that is correct on average but ruined by vibration. Fusing them is not a software trick but a small dynamic system whose own differential equations decide how quickly drift is removed and how much vibration gets through. Treating it as an ODE rather than as three lines of code is what allows it to be tuned for a particular vehicle.
The differential equation: $$ \dot{\hat{\theta}} = \omega_m - \hat{b} + k_p(\theta_a - \hat{\theta}), \qquad \dot{\hat{b}} = -k_i(\theta_a - \hat{\theta}) $$
\(\hat\theta\) is the estimated tilt angle and \(\hat b\) the estimated gyroscope bias, \(\omega_m = \dot\theta + b + n_g\) is the measured angular rate (true rate plus bias plus noise), \(\theta_a = \theta + n_a\) is the tilt angle inferred from the accelerometer, and \(k_p\) and \(k_i\) are the proportional and integral filter gains.
If your group wants more: Derive the steady-state scalar Kalman filter for the same two-state problem and show that it has exactly this structure, with the two gains now fixed by the ratio of gyroscope noise to accelerometer noise rather than chosen by hand.
Piezoelectric stacks position atomic force microscope stages and fuel injectors to nanometre precision, but the displacement produced by a given voltage depends on where the voltage has been: the extension and retraction curves differ by 10 to 15 per cent of full stroke. This memory cannot be captured by any algebraic function of voltage, only by an extra state carrying its own differential equation. Modelling it correctly is the difference between nanometre and micrometre open-loop positioning.
The differential equation: $$ m\ddot{x} + c\dot{x} + kx = k\bigl(d_p v(t) - h\bigr), \qquad \dot{h} = \alpha d_p\dot{v} - \beta|\dot{v}|h - \gamma\dot{v}|h| $$
\(x\) is the actuator displacement, \(m\), \(c\) and \(k\) the effective mass, damping and stiffness of the stack with its load, \(v(t)\) the applied voltage, \(d_p\) the piezoelectric coefficient relating voltage to unconstrained extension, \(h\) the internal hysteretic state measured in units of displacement, and \(\alpha\), \(\beta\) and \(\gamma\) the shape parameters of the hysteresis loop.
If your group wants more: Build an inverse model that takes a desired displacement and computes the voltage required, by integrating the Bouc-Wen equations in the reverse direction, then show by simulation how much of the hysteresis it removes in open loop and what happens when the parameter estimates are 10 per cent wrong.
A cooled stirred-tank reactor running an exothermic reaction can have three possible steady states at the same feed conditions, and the middle one is unstable. Plants have been lost because a small change in coolant temperature moved the reactor from the low-conversion state to the ignited one. The whole behaviour is contained in two coupled nonlinear ordinary differential equations.
The differential equation: $$V\frac{dC_A}{dt}=q\left(C_{A0}-C_A\right)-k_0e^{-E/RT}C_AV$$ $$V\rho c_p\frac{dT}{dt}=q\rho c_p\left(T_0-T\right)+\left(-\Delta H_r\right)k_0e^{-E/RT}C_AV-UA\left(T-T_c\right)$$
C_A is the reactant concentration in the tank, T the reactor temperature, V the reactor volume, q the volumetric feed rate, C_A0 and T_0 the feed concentration and temperature, k_0 and E the pre-exponential factor and the activation energy, R the gas constant, rho and c_p the density and specific heat capacity, (-Delta H_r) the heat of reaction, U and A the overall heat transfer coefficient and the jacket area, and T_c the coolant temperature.
If your group wants more: Add the jacket as a third state, or a proportional controller acting on the coolant temperature, and locate the Hopf bifurcation at which the reactor changes from a stable steady state to sustained temperature oscillations.
Industrial catalysts are porous pellets, and in a fast reaction the reactant never reaches the centre, so much of the expensive metal does nothing. The effectiveness factor tells the designer how much of the pellet is actually working, and whether to buy smaller pellets or a shell-impregnated catalyst at the cost of extra pressure drop. It comes from a boundary value problem in spherical or cylindrical coordinates.
The differential equation: $$\frac{D_e}{r^{2}}\frac{d}{dr}\left(r^{2}\frac{dC_A}{dr}\right)-kC_A=0,\qquad \left.\frac{dC_A}{dr}\right|_{r=0}=0,\quad C_A(R)=C_{As}$$ $$\phi=\frac{R}{3}\sqrt{\frac{k}{D_e}},\qquad \eta=\frac{1}{\phi}\left(\frac{1}{\tanh 3\phi}-\frac{1}{3\phi}\right)$$
C_A(r) is the reactant concentration at radius r inside the pellet, D_e the effective diffusivity of the pore structure, k the first-order rate constant per unit pellet volume, R the pellet radius, C_As the concentration at the outer surface, phi the Thiele modulus and eta the internal effectiveness factor.
If your group wants more: Solve the non-isothermal pellet by adding the energy balance with the Prater relation, and identify the parameter range in which the effectiveness factor exceeds unity and multiple pellet steady states appear.
A real packed bed is neither a perfect plug flow reactor nor a stirred tank: backmixing around the packing spreads the residence times and lowers conversion. Engineers measure that spread with a tracer test and turn it into a Peclet number for design. The same partial differential equation describes both the tracer pulse and the reacting species.
The differential equation: $$\frac{\partial C}{\partial t}=D_{ax}\frac{\partial^{2}C}{\partial z^{2}}-u\frac{\partial C}{\partial z}-kC$$ $$uC_{0}=uC(0^{+},t)-D_{ax}\left.\frac{\partial C}{\partial z}\right|_{0^{+}},\qquad \left.\frac{\partial C}{\partial z}\right|_{z=L}=0,\qquad Pe=\frac{uL}{D_{ax}}$$
C(z,t) is the concentration at axial position z and time t, D_ax the axial dispersion coefficient, u the interstitial velocity, k the first-order rate constant, L the bed length and C_0 the feed concentration; the two conditions are the Danckwerts closed-vessel boundary conditions and Pe is the Peclet number measuring how close the bed is to plug flow.
If your group wants more: Add intraparticle mass transfer resistance so that the bed is described by two coupled equations for bulk and particle concentration, and show how the apparent dispersion coefficient inferred from a tracer test is inflated by uptake into the particles.
Steady-state design says how many trays are needed, but not how long the column takes to return to specification after the feed composition shifts, and that is what decides how much off-specification product is made. A column is a chain of interacting tray holdups, so its dynamics is a large system of coupled ordinary differential equations whose slowest eigenvalue sets the settling time.
The differential equation: $$M_n\frac{dx_n}{dt}=L_{n+1}x_{n+1}+V_{n-1}y_{n-1}-L_nx_n-V_ny_n,\qquad y_n=\frac{\alpha x_n}{1+(\alpha-1)x_n}$$
x_n and y_n are the liquid and vapour mole fractions of the light component on tray n, M_n the molar liquid holdup on that tray, L_n and V_n the liquid and vapour molar flows leaving it, and alpha the relative volatility used in the constant-relative-volatility equilibrium relation.
If your group wants more: Add linear liquid hydraulics through the Francis weir relation so that flow and composition dynamics are coupled, then close a single-ended composition loop with a PI controller and find the gain at which the column becomes oscillatory.
Every temperature loop in a plant has a delay: the sensor sits downstream, the jacket takes time to respond and the analyser reports late. Delay is what limits how tightly a loop can be tuned, and a controller tuned as if there were none will oscillate or trip the reactor. The analysis is done with Laplace transforms and complex arithmetic rather than by trial and error on a live plant.
The differential equation: $$\tau\frac{dT}{dt}+T=K\,T_j(t-\theta),\qquad G(s)=\frac{Ke^{-\theta s}}{\tau s+1},\qquad G_c(s)=K_c\left(1+\frac{1}{\tau_Is}+\tau_Ds\right),\qquad 1+G(s)G_c(s)=0$$
T is the measured reactor temperature, T_j the jacket or coolant temperature that the controller manipulates, tau the process time constant, K the steady-state gain, theta the dead time, s the Laplace variable, K_c, tau_I and tau_D the controller gain, integral time and derivative time, and the last expression is the closed-loop characteristic equation.
If your group wants more: Replace the delay by a first- or second-order Pade approximation, show how well the root locus of the approximated system predicts the true stability limit, then design a Smith predictor and quantify how much faster the loop can then be tuned.
Every method in Chapter 1 assumes that a solution exists and that it is the only one. Picard's theorem says exactly when that assumption is safe, and one cube-root equation shows what happens when the Lipschitz condition fails. The group turns an existence theorem into something they can compute by hand.
The differential equation: $$\frac{dy}{dt} = y^{2/3}, \qquad y(0) = 0$$
Here \(y(t)\) is the unknown function of the independent variable \(t\), and the right-hand side \(f(t,y)=y^{2/3}\) is continuous everywhere but is not Lipschitz in \(y\) at \(y=0\).
If your group wants more: Prove Peano's existence theorem using Euler polygons and the Arzela-Ascoli theorem, and explain why it gives existence with no uniqueness at all.
Chapter 6 expands initial data in sines because the boundary conditions happen to be simple. Sturm-Liouville theory explains why this works far more generally: the operator is self-adjoint in a weighted inner product, so its eigenfunctions are orthogonal and complete. The group shows that a Fourier series is an eigenfunction expansion, not a trick.
The differential equation: $$\frac{d}{dx}\!\left(p(x)\frac{dy}{dx}\right) + q(x)\,y + \lambda\, w(x)\, y = 0, \quad a \lt x \lt b; \qquad \text{concrete case } y'' + \lambda y = 0,\ y(0)=0,\ y'(L)=0 \ \Longrightarrow\ \lambda_n = \left(\frac{(2n-1)\pi}{2L}\right)^{2},\quad y_n(x)=\sin\!\left(\frac{(2n-1)\pi x}{2L}\right),\ n=1,2,\dots$$
\(y(x)\) is the unknown function on \([a,b]\), \(\lambda\) is the eigenvalue parameter, \(p>0\) and \(w>0\) are the coefficient and weight functions, \(q\) is the potential term, and \(L>0\) is the length of the interval in the concrete case.
If your group wants more: Treat a singular Sturm-Liouville problem, either Legendre's equation on [-1,1] or Bessel's equation on [0,a], and explain exactly what changes in the theory when p vanishes at an endpoint.
Linearisation classifies an equilibrium but says nothing about how far its influence reaches, and it fails outright in the borderline undamped case. A Lyapunov function, which here is just the physical energy, settles both questions on the whole plane. The group uses the pendulum to see where the basin of attraction actually ends.
The differential equation: $$mL^{2}\ddot{\theta} + cL^{2}\dot{\theta} + mgL\sin\theta = 0, \qquad V(\theta,\omega)=\tfrac{1}{2}mL^{2}\omega^{2}+mgL\,(1-\cos\theta), \qquad \dot{V} = -cL^{2}\omega^{2} \le 0$$
\(\theta(t)\) is the angle from the downward vertical and \(\omega=\dot\theta\) the angular velocity; \(m\) is the bob mass, \(L\) the rod length, \(g\) the gravitational acceleration and \(c\ge 0\) the linear damping coefficient, while \(V\) is the total mechanical energy.
If your group wants more: Estimate the basin of attraction using the largest sub-level set of V contained in a suitable region, then compare that estimate with the true basin computed numerically and quantify how conservative it is.
On some perfectly smooth problems an explicit solver is forced to take steps thousands of times smaller than the accuracy of the answer requires. The reason is stability, not accuracy, and it is decided by where the eigenvalues sit relative to the method's stability region. The group demonstrates the collapse and then removes it with an implicit method.
The differential equation: $$y' = \lambda\,(y-\cos t) - \sin t, \quad y(0)=0, \quad \lambda=-1000 \qquad \Longrightarrow \qquad y(t)=\cos t - e^{\lambda t}$$
\(y(t)\) is the unknown function of time \(t\) and \(\lambda<0\) is a fixed constant that sets the rate of the fast transient; the exact solution is a slow part \(\cos t\) plus a transient \(-e^{\lambda t}\) that is negligible after about \(5/|\lambda|\).
If your group wants more: Show that the trapezoidal rule is A-stable but not L-stable by examining the amplification factor as the real part of h lambda tends to minus infinity, and demonstrate the resulting spurious oscillation on this very problem.
The course inverts Laplace transforms by matching entries in a table. The inversion is really a contour integral in the complex plane, and the residue theorem turns it into a finite calculation that explains the shape of the answer. Resonance, in particular, is simply a pole of order two sitting on the imaginary axis.
The differential equation: $$y''+4y=\sin 2t,\quad y(0)=y'(0)=0 \;\Longrightarrow\; Y(s)=\frac{2}{(s^{2}+4)^{2}}, \qquad y(t)=\frac{1}{2\pi i}\int_{\gamma-i\infty}^{\gamma+i\infty} Y(s)e^{st}\,ds = \frac{1}{8}\left(\sin 2t - 2t\cos 2t\right)$$
\(y(t)\) is the displacement of an undamped oscillator of natural frequency 2 driven at its own frequency, \(s\) is the complex transform variable, \(Y(s)\) is the Laplace transform of \(y\), and \(\gamma\) is any real number to the right of every singularity of \(Y\).
If your group wants more: Invert a transform with a branch point, such as \(e^{-x\sqrt{s}}/s\), using a keyhole contour, and identify the resulting complementary-error-function solution of the heat equation on the half-line.
Learners commonly treat a differential equation as a formula to be hunted down and cannot say what a solution is. A lesson that begins with the slope field forces the idea that the equation prescribes a direction at every point and that a solution is a curve following those directions. The group designs the lesson, teaches it to classmates, and revises it on the evidence.
The differential equation: $$\frac{dP}{dt} = 0.5\,P\left(1-\frac{P}{100}\right), \quad P(0)=10 \qquad \Longrightarrow \qquad P(t)=\frac{100\,e^{t/2}}{9+e^{t/2}}$$
\(P(t)\) is the population at time \(t\) in years, 0.5 per year is the intrinsic growth rate and 100 is the carrying capacity in the same units as \(P\).
If your group wants more: Redesign the same lesson for an equation with no elementary closed form, such as dy/dx = x^2 + y^2, and argue what learners can still establish without solving it.
When the forcing term already solves the homogeneous equation, the usual guess collapses and learners nearly always conclude they have made an arithmetic slip. The difficulty is conceptual, about the solution space and linear independence, not about algebra. The group builds an instrument that separates students who know the rule from students who know why the rule exists.
The differential equation: $$y''-4y'+4y=e^{2t}, \qquad y_h=(C_1+C_2t)e^{2t}, \qquad y_p=\tfrac{1}{2}t^{2}e^{2t}, \qquad y=(C_1+C_2t+\tfrac12 t^{2})e^{2t}$$
\(y(t)\) is the unknown function of \(t\), \(y_h\) is the general solution of the homogeneous equation whose characteristic root \(r=2\) is repeated, \(y_p\) is a particular solution, and \(C_1,C_2\) are arbitrary constants.
If your group wants more: Extend the instrument to y'' + 4y = sin 2t, where the resonance is oscillatory rather than exponential, and test whether students who corrected the first case transfer the reasoning to the second.
A student can run the Laplace table forwards and backwards without being able to say what a switch turning on at t = 2 does to a physical system. A test made only of transform exercises cannot tell those two students apart. The group builds a paired-item instrument in which each procedural question has a conceptual twin resting on the same mathematics.
The differential equation: $$y''+3y'+2y=u(t-2),\quad y(0)=y'(0)=0 \;\Longrightarrow\; Y(s)=\frac{e^{-2s}}{s(s+1)(s+2)}, \qquad y(t)=u(t-2)\left[\tfrac12 - e^{-(t-2)} + \tfrac12 e^{-2(t-2)}\right]$$
\(y(t)\) is the response of a damped system at rest until time 2, \(u(t-2)\) is the unit step that switches on there, \(s\) is the transform variable and \(Y(s)\) is the Laplace transform of \(y\).
If your group wants more: Add an item using the impulse forcing delta(t-2) in place of the step, where y stays continuous but y' jumps by 1, and check whether students who explained the step case can predict the impulse case.
School learners meet exponential decay as a formula with no origin. A cup of hot water, a thermometer and a timer give them the differential equation itself, and Euler's method lets them build the solution using arithmetic they already have. The group runs the experiment, teaches the modelling lesson, and compares the analytic and numerical treatments honestly.
The differential equation: $$\frac{dT}{dt}=-k\,(T-T_a), \quad T(0)=T_0 \qquad \Longrightarrow \qquad T(t)=T_a+(T_0-T_a)e^{-kt}, \qquad T_{n+1}=T_n-h\,k\,(T_n-T_a)$$
\(T(t)\) is the temperature of the water at time \(t\), \(T_a\) is the constant ambient temperature, \(T_0\) the initial temperature, \(k>0\) the cooling constant in units of inverse time, \(h\) the Euler step size and \(T_n\) the approximation to \(T(nh)\).
If your group wants more: Add a second cup in a different container and model the two temperatures as a coupled 2x2 linear system, solving it with eigenvalues and linking the lesson forward to Chapter 4.
Eigenvectors are invisible in the algebra but plain in a phase portrait once the learner can move the initial point and watch what happens. A tool that lets the user drag the starting point and change the matrix entries makes the classification of equilibria something to be noticed rather than memorised. The group builds such a tool and then finds out whether it actually helps.
The differential equation: $$\begin{pmatrix} x' \\ y' \end{pmatrix} = \begin{pmatrix} 1 & 2 \\ 3 & 2 \end{pmatrix}\begin{pmatrix} x \\ y \end{pmatrix}, \qquad \lambda_1=4,\ \mathbf{v}_1=\begin{pmatrix} 2 \\ 3\end{pmatrix}, \qquad \lambda_2=-1,\ \mathbf{v}_2=\begin{pmatrix} 1 \\ -1\end{pmatrix}, \qquad \mathbf{x}(t)=c_1\mathbf{v}_1e^{4t}+c_2\mathbf{v}_2e^{-t}$$
\(x(t)\) and \(y(t)\) are the two unknown functions of time, the matrix entries are the constant coefficients, \(\lambda_1,\lambda_2\) are its eigenvalues with eigenvectors \(\mathbf{v}_1,\mathbf{v}_2\), and \(c_1,c_2\) are constants fixed by the initial condition.
If your group wants more: Extend the applet to complex and to repeated eigenvalues, and animate a parameter that carries the system across a boundary in the trace-determinant plane so that the change of portrait can be watched.