Programs
All programs were created by David Moore unless otherwise stated.
Simulation of an electron gun. I used numeric.js to solve Laplace's equation on the domain, and then calculated electron trajectories from an equation relying only on the direction vector $\hat{s}$ and potential difference $\varphi(x)-\varphi_0$.
This is a fun exercise in linearizing systems. A system of ideal springs in more than one dimension is a nonlinear system because of the square root taken in finding the distance. I wrote a Mathematica program which, given a mesh of springs, finds the eigenmodes of the system under small oscillations.
A disk rolling without slipping on a plane is an example of a nonholonomic system. This program solves for the dynamics of a massive disk, using the tools one normally uses for nonholonomic systems. The algebra was done in Mathematica and the plotting is done in three.js.
This is a solution of Kepler's problem. Most books on mechanics will give the full solution of the problem, but won't match the solution to initial position/velocity conditions. That turns out to have a few differences, so it's a good exercise to implement a program like this.
You can pass complex values into the logistic map and investigate its limiting properties. What is plotted is the real part of the limiting value.
I wrote this while thinking about ways to simulate the effects of the Earth's tidal bulges on the moon. What's plotted is how much the gravitational field diverges in the perpendicular/angular direction from a perfectly radial field.
I'm still on a quest to better visualize Riemannian manifolds. OK, sitting somewhere and trying to draw functions including their branch cuts is informative, but it leaves one with no idea how, for example, the exponential function looks anything like the logarithmic function. Usually, you say, "log is shaped like a spiral, or a depressed sheet" and "exp is shaped like a sheet that's flat on one side and rumpled on the other." But they're different views of the same manifold! This program lets you walk around on a Riemannian manifold and view the area close to the plane tangent to the manifold wherever you're at. I don't think it turned out too well, and I'm still working on different ways to view Riemannian manifolds.
To calculate the linked images, I essentially used the shooting method to calculate the eigensystem of the situation and then approximate initial waveforms. I did this because I was tired of the type of problems I kept running into when doing a general numerical PDE solution.
I made this for fun while answering a question on stackexchange. The most interesting thing about this one is that I calculated all the intersections using vector math. Sometimes it's easier to think in terms of angles, but you end up having to deal with a lot of special cases. Vector math is more well behaved.
The idea here is that you replace the 1/r^2 attraction law which governs the acceleration of two gravitationally interacting objects with a law governing the velocities of the objects.
As above, but with three bodies. I have a feeling this must be explicitly solvable, but I haven't tried my hand at it yet.
A three-body gravitational system is chaotic. This interactive program (javascript+processing.js) draws the path of particles over all of time, and varies the initial velocity of a single particle based on the mouse position. The whole simulation is in a zero-momentum frame, centered on the center of mass.
An interactive de Jong attractor drawer. A de Jong attractor is [the limit of a] chaotic trigonometric map in two dimensions.
Described in this blog post, I wrote a program which analyzes the dynamics of an otherwise simple conservation of momentum problem.
This is a noninteractive small program, which, for fun, simulates the simple wave equation where a wave pulse is emitted from a moving source. The whole simulation is then given a Lorentz boost. What is a standard reflection (frequency'=frequency) in one frame, is a relativistic doppler shift in another!
It's always fun to mess around with iterated function systems. This is a basic app which I'll expand on in the future.
It's easy to show that the cantor set has the same cardinality of the reals, so of course this means there is some function which maps the cantor set to the filled in unit square. Still, actually constructing such a function is amazing! The program shows one such construction.
Potential energy curves are useful for analyzing kinematics. The effective potential energy of an object in orbit isn't a true potential energy curve, because it consists partly of kinetic energy. By taking into account the conservation of angular momentum, you can rewrite energy conservation to depend on only the radius and radial velocity, and then you can do everything you can with a normal potential energy curve. Khan academy program at: http://www.khanacademy.org/cs/effective-potential-energy-perturbation-orbital-mechanics/1193669002.
This is an awesome-looking, animated program, that's easy to write. It was one of the first things I did on Khan Academy's CS website..
This is a simulation of a bunch of particles undergoing gravitational attraction. Unfortunately it's very slow, and I don't think there's much I can do to speed it up with javascript. Khan academy program at: http://www.khanacademy.org/cs/challenge-modeling-accretion-disks/1188226284.
Quadtrees can be used to speed up gravitational simulations. Unfortunately, in javascript building a quadtree appears to take a long time. Traversing it seems fast enough. I haven't done much cost-analysis on this program though.
This program compares simulated results with the results of an exact solution of a damped driven oscillator. The general/exact solution for a damped driven oscillator is a bit long, and I haven't implemented solutions for critical damping. Khan academy program at: http://www.khanacademy.org/cs/driven-oscillator/1115275237.
The book I was studying, A. P. French's Vibrations and Waves, discussed the normal modes of a chain of ideal springs (with small displacements). One practice problem was to find the steady-state solution of a chain with a driven boundary condition instead of a stationary one. This is the solution to that problem, with the exact solution compared to a simulated one with exact initial conditions.
This program finds the normal modes of a 2-mass spring system, given initial conditions. It visually shows the addition of modes, to get the final real result. Khan academy program at: http://www.khanacademy.org/cs/coupled-oscillators/1120086006.
Nothing too special, just a 3D projection! I was trying to find a parametric equation for a 3D knot, and I didn't use mathematica or google, for whatever reason.
This program takes a bunch of particles and uses quadtree optimization to approximate the forces between every particle. It's pretty efficient, written in Java, and you can find an applet online here: /archive/programs/javagrav/.
Visually exciting? Not really. Technically exciting? Absolutely! This program takes an initial position and velocity function, takes its Fourier transform, and then uses that to find the exact solution to the wave equation with those initial conditions and stationary boundary conditions (to some number of terms).
I always had a particularly difficult time solving the 2D wave equation using standard techniques, so I took a look at the code at http://www.falstad.com/ripple/ and implemented the same thing. It rotates numbers in the complex plane, and it seems to run fast and stable! I did manage to show that this type of integration is the same as the standard wave equation (up to a constant factor) in the limit as the timestep goes to zero, but it still seems a bit magical. Cool stuff! (oh, yeah, and the program has free-boundary edges with two spots driven periodically)
This is a very simple inverse kinematics problem. It looks great when animated.
In the Khan Academy intro CS programs there's one that draws a cool-lookin' curve. I wanted to draw the curve differently, and to do so I needed to find a parametric equation for the curve. I assumed the original program was really just approximating another curve, and used a bit of calculus to get the exact equation (the original program used a "rate of change" variable, so I had to integrate). The new one is cool because you can generalize it by changing the constants of the curve.
This is a great program showing what happens to a wave when it meets a different medium. It's a basic simulation, but it shows the relationship between the magnitudes and sizes of the reflected/transmitted waves.
Algebraic Numbers
A C++ program that generates, solves, and plots the roots of polynomials.
This program is based off of a program on wikipedia, but was re-implemented with some additions in C++ with SDL. The code was updated to be object-oriented, with the ability to drag and zoom in real time.
This program generated the image at the ams visualinsight blog, and I've put the source code on sourceforge.
Multithreaded Mandelbrot Renderer
A multithreaded C++/SDL program that plots a supersampled black and white image of the Mandelbrot set.
Wolfram 1D Cellular Automata
A C++/SDL program that plots a supersampled black and white simulation of the Wolfram 1-dimensional cellular automata. Among these images are rules 30, 169, and 110. Two simulations were performed for the colored ones, one with random initial conditions, and another with the same initial conditions, with the middle cell flipped. The difference between the two simulations is colored, and was done with GIMP.
A bit of Mathematica code that draws the result of a hyperbola projected on a sphere, Riemann-sphere style.
There are .gif animations of this on the gif page.
It turns out that, due to the mass of the rope, a bungee jumper will accelerate slightly faster than gravity. This simulation has a bungee jumper on a heavily damped chain with about the same mass as the jumper.
A vector/velocity field generated based on iterations of the Julia set of a quadratic polynomial.
There are .gif animations of this on the gif page.
This is a really cool project and it's a shame that more people aren't aware of it. Newton's method, given a function and its derivative, defines a map from the reals to the reals, which hopefully causes points to converge to roots as you apply this map again and again. So, if you generate a polynomial with integer coefficients, you can start at some point, apply Newton's method, and then when the calculation converges (and you've verified it's a root), you can color the starting point of the calculation based on which root it converges to, and how long it took to converge. This is exactly that, but with complex-coefficient polynomials, and complex starting points.
Lorentz Attractor
A simulation of the classic Lorentz attractor my friend wrote in Java.
Lindenmayer System
A Lindenmayer-System generator/drawer I wrote in high school computer science, using Java.
Base-Motif Fractal and Fractal Sweeps
Base-motif fractals are a type of recursive replacement fractal described in Mandelbrot's "The Fractal Geometry of Nature". He also defines fractal sweeps, which allow you to flip patterns horizontally or vertically. This is a cool program, the original one I wrote in DBPro at http://forum.thegamecreators.com/?m=forum_view&t=180312&b=6, and I also made a similar one (which doesn't handle fractal sweeps) on the Khan Academy CS website http://www.khanacademy.org/cs/linedrawer1/1074070969.
There are .gif animations of this on the gif page.
Squig Generator
A Squig is defined in Mandelbrot's Fractal Geometry of Nature from subdividing triangles. The program was written with C++/SDL.
There are .gif animations of this on the gif page.
Minutephysics posted a video about why the night sky is dark, so I thought I'd go through some related math to generate my own night sky. I figured out what the probability was of finding a star in a given volume was, found the volume element for a standard screen projection, and then used the resulting probability distribution to find out how far away a star should be for each given pixel. I was disappointed in the render at first, but now I think it's pretty awesome!
Efficient C++ N-Body Simulation
This does essentially the same thing as my previous n-body code, using object oriented quadtrees. However, this simulation is written in C++, and it handles repulsion (and collision damping, so that the collisions aren't elastic!) as well as attraction. You can see a video of this simulation with 1,000,000 particles here: http://www.youtube.com/watch?v=nWJLG7s3Cfc. Note that, in the last part of the video, with 20x zoom, (most visible if you watch the video in high quality) the small clumps of particles (clumps on the order of tens or hundreds of particles) behave more like solid, rigid bodies.
There are still noticeable inaccuracies in the simulation, and I plan on correcting these, as well as adding new features. Switching from per-pixel drawing to OpenGL is one goal. Improving force-calculation methods is another. Simple CPU multithreading is also possible. I'd love to start doing massively-parallel work but I don't plan on studying this in the near future.
Mass on a Rail
This is from a bit of Mathematica code I wrote to handle a mass sliding on a frictionless rail, under the influence of gravity. I was practicing Mathematica and was also was refreshing my memory on some vector stuff, so the program actually uses NDSolve[]'s magic on an ugly vector differential equation (well, the equation is a scalar equation, but it's the dot of the derivatives of two vectors), even though energy considerations would make things WAY easier. This wasn't a waste of time though, because Mathematica did the work for me, so it was just as easy. But hey, this makes it easy to, say, vary the force direction over time, or to add friction.
There are .gif animations of this on the gif page.
Logistic Map - Cobweb Diagram & Lyapunov Exponent
This was written with Mathematica. It draws cobweb diagrams, and also calculates the Lyapunov exponent, of the logistic map, for different values of r. It's actually written in a way that it can calculate the Lyapunov exponent and cobweb diagram of ANY real map, and it wouldn't take many modifications to calculate maps of vectors. (The map must be differentiable in order to calculate the Lyapunov exponent. the Lyapunov exponent is just a sum of logarithms of absolute values, and you could just replace "abs" with "norm" to get an equation for vectors. Cobweb diagrams would need changes too, but it could work.)
This program was based off a section in "Introduction to Computational Physics" by DeVries. There are .gif animations of this on the gif page.
Schrodinger Equation with Multiple Wells
This was written with Mathematica. It calculates solutions to the time independent schrodinger equation, for a single particle, with multiple potential wells, in 1D. The challenge is calculating an allowed energy state, that leads to a solution which decays exponentially on either end. When it works, it works perfectly! It works best for low numbers of wells, and low separation between the wells.
I said, "when it works". It turns out that this can get difficult. With 40 wells, this program comes up with an equation that is the result of a matrix to the 40th power, takes that matrix and contracts it with a row and column vector to get a scalar, and then finds the roots of that scalar equation. The power of identical matrices can be sped up through diagonalization, but the root finding presents another challenge. The matrices contain sine and cosine terms, but they also contain sinh and cosh terms! These are exponential, so they can get very large very fast. Combine this with a 40th power something or other, and you have problems. At x=1, f(x) could be 10^300. At x=1.1, f(x) might be the zero you're looking for, and at x=1.2 it might be -10^300. How do you find the roots to such an equation numerically? I have no clue. And whenever Mathematica attempts to, it crashes. Not gracefully either, despite my best attempts, when it crashes it freezes and I need to alt-f4.
This program was based off a section in "Introduction to Computational Physics" by DeVries.
I'm saving work related to "Introduction to Computational Physics" by DeVries. There is too much content to show on here, in PDF form the mathematica notebooks are 83 pages total. These include the above two projects and many others. The mathematica notebooks are on sourceforge.
Chapter PDFs:
- ch1 Introduction.pdf
- ch2 Functions and Roots.pdf
- ch3 Interpolation and Approximation.pdf
- ch4 Numerical Integration.pdf
- ch5 Ordinary Differential Equations.pdf
Relativistic Transformations
This just contains some miscellaneous work with Lorentz transformations. A bunch of equations are re-derived using mathematica, and some things like the optical effects of relativistic speeds (pictured) and the properties [and optical effects of] throwing a ball on a spaceship as observed from earth are explored.
See also: "Bungee Jumper Acceleration"
A Lagrangian for a hanging massive chain predicting that the chain accelerates more quickly than gravity alone accounts for. Here are some good questions: What if the chains were unlinked and fell into place by some tiny but strong threads? Clearly the links can't accelerate faster than gravity if they can't exert forces on each other. So where in calculating the original equations of motion did we assume that the chain was linked in the first place?
So, Purcell's book on electricity and magnetism goes into a little bit of depth about modelling a neutral conductor as a bath of positive ions and electrons, which are treated as classical billiard-ball particles. Normally you might think to model this as a stationary grid through which electrons are bouncing. I chose not to model it this way because a stationary grid would ignore the forces on the positive ions, which - even though Purcell's book talks all about transporting electrons, forces of electrons, accelerations of electrons and so on - is just as important for momentum conservation! If this sounds dumb, I think of it this way: there's a huge total force on the positive ions in a large electric field. What's stopping the conductor from leaping out of the laboratory is the force of the humble electron slamming in to the rest of the conductor, leaving on average a zero net change in momentum. The current is averaged over time and I do have to artificially damp the system so that the energy doesn't run away to infinity.
This is from the book by Landau and Lifshitz. So far I haven't written too detailed a simulation.
