Having laid the groundwork in part 1, it is now time to actually simulate the interior of a white dwarf star. I’ve decided to use C# for this. I really like the C# language and with Mono it’s fairly trivial to make it cross-platform, which is nice.
I covered the basics needed in my post Numerical approximations to differential equations. In this case we have a system of coupled differential equations, meaning that each differential equation depends on the other. However the basic approach we’ll use to solve them are very similar to the way I described in that post.
Integrator framework implementation
In order to make the code usefull for other problems, I decided to try to make it a bit generic. I declared a general interface for ODE integrators. The integrator will maintain a state comprised by the independent variable, say time, and one or more dependent variables, like position and velocity. Of course, the interpretation of the variables is of no importance to the integrator.
After setting the initial state one repeatedly calls the Update method to advance the integrator. In order to perform the update the integrator will require the derivatives of the dependent variables. So, the Update method requires a function pointer (or delegate in C#) which it can use to calculate the derivatives for an arbitrary state. This makes it possible to use integrators which internally perform multiple steps.
1 2 3 4 5 6 7 | public interface IIntegrator { void SetInitialState(double x, double[] y, double h); void GetState(out double x, out double[] y); void Update(DerivativesDelegate derivatives); } |
To make my life easier, I made an abstract base class which implements IIntegrator. It has an internal state and implements the methods for setting and getting the state. Specific integrators can then just override the Update method.
The first integrator to be implemented is, of course, the forwards Euler method. Since the implementation is so simple, it’s nice to have during debugging in order to determine if there’s a problem in the integrator or the calculation of the derivatives. The implementation looks like this.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 | public class EulerIntegrator : IntegratorBase { public EulerIntegrator() { } public override void Update(DerivativesDelegate derivatives) { double[] dydx; // get derivatives derivatives(x, y, out dydx); for (int i = 0; i < n; i++) { // apply forward Euler method y[i] = y[i] + h*dydx[i]; } x += h; } } |
It first gets the derivatives \(\dx{y}\) at the current position \(x\), then uses these to take a step forwards. The length of the step is \(h\), which is set by SetInitialState.
Small note regarding the derivatives delegate. Since it’s not possible to pass a const array in C# the derivatives method can modify the internal state of the integrator. In order to avoid this the integrator should pass a copy of the array or a read-only view of some sort. However I felt this would only clutter the code for this project.
The fourth-order Runge-Kutta integrator is a bit more complicated. It takes four “trial steps” in order to perform a single update. Each of these steps are done using the regular forwards Euler method but with different step lengths and derivatives. The steps are as follows:
- Use the derivatives at \(x\) to perform a half-step to \(x+\frac{h}{2}\). Find the derivatives.
- Use derivatives at \(x+\frac{h}{2}\) to perform another half-step from \(x\). Find the derivatives again.
- Use updated derivatives at \(x+\frac{h}{2}\) to take a full step from \(x\) to \(x+h\) and find the derivatives.
- Compute a weighted average of all four derivatives, use this to step from \(x\) to \(x+h\).
Instead of using the finite difference method to approximate the differential equation, it is derived by using the fundamental theorem of calculus and approximating the integral using Simpson’s rule. In addition the midpoint is found by averaging two approximations.
Enough boring background, here’s my implementation.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 | public class RK4Integrator : IntegratorBase { public RK4Integrator() { } public override void Update(DerivativesDelegate derivatives) { double[] k1, k2, k3, k4; double[] yt = new double[n]; double h2 = h / 2; double h6 = h / 6; double x2 = x + h2; double xh = x + h; // compute intermediary steps // get derivatives at x derivatives(x, y, out k1); // get derivatives at x+h/2 for (int i = 0; i < n; i++) { yt[i] = y[i] + h2*k1[i]; } derivatives(x2, yt, out k2); // find new derivatives at x+h/2 for (int i = 0; i < n; i++) { yt[i] = y[i] + h2*k2[i]; } derivatives(x2, yt, out k3); // find derivatives at x+h for (int i = 0; i < n; i++) { yt[i] = y[i] + h*k3[i]; } derivatives(xh, yt, out k4); // get final state by using weighted // average of derivatives for (int i = 0; i < n; i++) { y[i] = y[i] + h6*(k1[i] + 2*k2[i] + 2*k3[i] + k4[i]); } x = xh; } } |
Since the implementation is more complicated, it’s a good idea to test that it all works. I used both integrators to solve a simple ballistic problem, which I could solve using pen and paper. This allowed me to verify their correctness and accuracy.
Stopping in time
The ballistic problem I used as a test case consists of simulating a canon ball fired upwards and find how long it takes for the canon ball to hit the ground. That means that we ideally would like to stop the integration at exactly the point when the ball hits the ground. One way of doing that would be to stop once the ball has passed beyond the ground level and then estimate the exact time of crossing based on the current and last position. This works fine for the ballistic problem. However, as I mentioned in part 1, in the white dwarf simulation \(\gamma\) isn’t defined when \(\bar{\rho}\) is negative. Ideally we want to stop before that happens. And the problem is even worse when using adaptive methods as the time step can be so large that you overstep by a huge amount.
In order to resolve this problem, I came up with the following scheme. It was inspired by how MATLAB handles this issue. After each update, the integrator checks if the solution has overstepped. If it has, it repeats the last update but with a smaller time step. It repeats the process until either the stopping criteria is fulfilled, the time step becomes very small or it reaches the maximum number of iterations allowed.
In order for the integrator to determine if the stopping criteria has been reached or exceeded, I introduced the following delegate.
1 2 3 4 | public enum StoppingCriteriaDirection { Both, Positive, Negative }; public delegate double StoppingCriteriaDelegate(double x, double[] y, out StoppingCriteriaDirection direction); |
In addition I introduced an overloaded Update method.
1 2 3 4 5 6 7 8 9 10 | public interface IIntegrator { void SetInitialState(double x, double[] y, double h); void GetState(out double x, out double[] y); void Update(DerivativesDelegate derivatives); void Update(DerivativesDelegate derivatives, StoppingCriteriaDelegate stoppingCriteria, out bool stopped); } |
The integrator then tries to locate the roots of the StoppingCriteriaDelegate. The StoppingCriteriaDirection is handy to prevent premature termination of the integration. Consider for instance if we in the ballistic problem had shot the cannon ball from ground level, but wanted to know how long it took to fall down on the roof of a house. Our implementation of the StoppingCriteriaDelegate would then return y - y_roof. In this case the stopping criteria would be met twice: first on the way up, in which case we want to keep going, and then on the way down. So we specify that we’re only interested in the case where the zero is approached from above, that is when the derivative is negative.
The neat thing about this implementation is that it reuses the existing framework in such a way that the new Update method only has to be implemented in the IntegratorBase class and can thus be used transparently on all descendants.
Simulation time
Seems I have to postpone the juicy bits for the next part so I can get this published before the apocalypse. In addition to the ballistic problem I’ve added a second test problem for which there is a simple analytical solution. Unlike the for the ballistic test the Runge-Kutta methods cannot reproduce this solution exactly.
In order to show how the framework can be used I’ll show the main loop of the ballistic simulation.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 | while (!stopped) { integrator.Update( // derivatives delegate delegate(double _x, double[] _y, out double[] dydx) { dydx = new double[2]; dydx[0] = _y[1]; // dx/dt = v dydx[1] = -9.81; // dv/dt = a in ms^-2 }, // stopping criteria delegate delegate(double _x, double[] _y, out StoppingCriteriaDirection direction) { direction = StoppingCriteriaDirection.Negative; return _y[0]; // stop when hitting the ground }, // stopped flag, true if stopping criteria is met out stopped); integrator.GetState(out x, out y); } |
The independent variable is x (time) while the dependent variables are y[0] (position) and y[1] (velocity). In this case I’ve opted to use anonymous methods for the two delegates. As you can see the stopping criteria delegate simply returns the position, with the ground being at \(x = 0\). Since we’re only interested in when the cannon ball hits the ground on the way down, the direction is set to Negative. Just to show how bad the Euler integrator is, here’s the test run from the ballistic simulation with an initial upwards velocity of \(13 m/s\) and a step length of \(0.1 s\):
Testing integrator 'ODE.Integrators.EulerIntegrator' Time: 2.749448 Position: 0.000000 Velocity: -13.972081 Testing integrator 'ODE.Integrators.RK4Integrator' Time: 2.650357 Position: 0.000000 Velocity: -13.000000 |
Since we’re not considering air resistance the cannon ball should have the same velocity when it hits the ground, except in a downwards direction of course. The time taken should be \(\sim2.650356779s\). It’s clear that the Euler integrator is having a difficult time with the “large” step time, however the RK4 integrator is spot on. This is not surprising given that it is a fourth order method and the ballistic problem is quadratic in nature.
If you want to play around you can download the entire source code here: WhiteDwarfSim.zip. I’ve developed it using MonoDevelop, but tested it on Visual Studio as well. So far the name is a bit misleading but I will justify it in part 3…