git clone https://djg11@bitbucket.org/djg11/numerical-methods-demos.git
And any updates can then be fetched using git pull.
Alphabetical Listing:
A matlab/octave program where the supposed inverse to a singular matrix has been computed. But we do see a large condition number (18 or so) which reveals there to be no significant figures of accuracy.
We convert from binary to ASCII characters, showing the mechanics of the printing a number in human-readable form. This is essentially a long division where the denominator is the largest power of ten required so that the leading digit is in the range 1 to 9.
We demonstrate power series economy using a Chebychev basis. You can uncomment various examples and also a direct Taylor implementation for comparison. We should see that, for the same number of terms, the Chebychev basis gives more even error, which means less error at the larger values of the input range. After all, we are generally just interested in the worst case error, so solving the worst case situation (perhaps at the expense of other less troublesome situations) is a saving.
We implement the lectured equations for the Cholesky decomposition. We use a little of our SimuSolve code to print the matrix and to print itself multiplied by its transpose. Since Cholesky is essentially a square-root function in matrix terms, we should get back the original when we multiply it by its transpose.
Not lectured - reveals a compiler bug.
Not lectured: This demonstrates what happens when we try to keep incrementing a 'float' forever. After a certain point, we see no further increase owing to underflow. We should expect the complete underflow to occur at about 10^7.
It will print this output 1.67772160e+07 (0x4b800000 or 151:000000)
Monte-Carlol Simulation to find Pi - Dart Throwing. We choose uniform random (X, Y) co-ordinates in the unit square. The percentage inside the unit circle will be proportional to Pi. What is the rate of convergence: very poor indeed. The program also include an implementation of antithetic variants. These were not lectured, and the implementation here has not been checked by DJG.
We have a simple set of flow equations, as in the lectured slides, where one of the components is non-linear. The non-linear component is replaced with a fantom flow and a fantom conductance. We see that we have to iterate to get convergence of the flow matrix solution.
C and ML versions of the floating-point to ASCII conversion function. Their complexity shows that saving data in human-readable form is expensive as well as being often unnecessary.
The C version implements the specific tests for NaN and Inf and so on. It does not use any floating-point operators.
The ML version is easier to understand than the C version. Since ML is strongly-typed, the code contains a call to floor. The floor function, in general, would do a lot of the work that is required for float-to-ASCII, and so its use should be classed as cheating. But owing to the work we do before and after, the floor calls used here will end up themselves implementing nothing more than a bit shift. It is worth noting that calls to floor (which themselves are implied in many programming languages where float and int can be tacitly mixed) often involve a lot of work and should be avoided for high performance.
The int-to-float folder, not lectured, contains further implementations, including many in Verilog RTL, which is reading for Part Ib CST students ECAD+ARCH course.
Two implementations of a simple iteration to find a root. We learn that where one direction expands the relative error on each iteration, the other direction will generally converge.
Variable-latency, 'school-method' long division of a pair of natural numbers. We need to recall that division is an expensive operation, as well as how to do it!
The mechanics of long multiplication. We see examples using Booth, Vehdic and school methods.
Loss of significance demo - A program that executes what we did by hand on our calculators in lectures. We take a recurring fraction and subtract off some digits at the start so that further digits of the mantissa are revealed. Beyond a certain point we see nonsense.
A program that forms L/U decomposition of a matrix using partial pivoting. Also code for forwards and backwards substitution.
Not lectured: This program shows how you can convert a +ve 'int' to a double only using bit operations. We abuse the C language type system using a union. The union type in C enables a piece of physical memory to be accessed with various types and no checking.
Two examples where the NaN IEEE feature is generated. The printed result is -nan and -nan.
This exercises the arithmetic system and shows the range expressible as 'float' numbers, including dnorm/infinity. Not mentioned in lectures, but worthy of private study.
Demonstration of noise caused by rounding error during different numerical evaluations of (x+1)^7. We see that one implementation is monotonic and the other is widely not. From this we should learn a useful principle when implementing library functions that are required to be monotonic: avoid subtraction of variable quantities.
Investigation of printing significant figures. How many are there in a single-precision number: 7, 8 or 9?
Quadratic equation solving, precisely as lectured, using an IF statement to avoid cancellation and loss of significance.
This folder contains two linear FDTD simulations. They are linear in terms of the equations of motion and also in the sense that the state vector is a linear array of elements. The difference is that one only contains one state variable per element and the other has two (position and velocity).
First-order system: FDTD simulation of heat conduction or RC network (no inductance).
Second-order system: FDTD simulation of an elastic material such as a violin string - this has velocity and position in state vector. The acceleration is not part of the second system's state since it can be computed statically from the positions. Note that elastic materials follow Hooke's law: the restoring force is proportional to the strain (ie. to the amount of deformation).
Simpson's Rule for Quadrature - We explore the best setting of h. We find it is at roughly the fifth root of macheps.
Not lectured in detail. A gif demo.
A tabulation (for plotting) of the error in Euler's method. For exponential decay, we get a modified rate of decay, but still a perfect exponential.
We explore the Taylor Series for the sine function. Although this converges for all x, it is not sensible to use it for large values of x owing to large amounts of cancelling with large terms of alternating signs: the diplodocus effect. Hence we should always apply range reduction first.
We see the unexpected results that can arise with equality comparisons on floats and doubles.
A hand-crafted finite difference simulation of the world's most successful rocket launcher.
We explore whether summing a series in the forwards direction is the same as the backwards direction. We find floating-point addition is non-associative. We also glance at Kahan's method for combining a pair of floats for greater precision.
We have a simple set of flow equations, as in the lectured slides, where one of the components has properties that vary over time. Hence we need to maintain a state vector and perform FDTD simulation, solving the flow equations for each time step. The component in question is a capacitor or bucket and this is replaced with a fantom flow and a fantom conductance.
A pair of second-order, linear FDTD simulations.
FDTD simulation of a two-pole ring oscillator. If we have a pair of cascaded integrators and feed some of the output back to the input with inversion of sign, we have the classic conditions for SHM.
Not lectured: The same goes for a pair of differentiators, but this does not make for a stable simulation. Moral: before starting a simulation, integrate your equations as needed to avoid differentiators.
A demonstration that a chaotic system need only have a single scalar in its state. Verhulst's logistic map is an iteration that just jumps around a lot. It has two meta-stable roots but we typically find neither of them.
A manual coding of the FDTD system containing several pipes and two buckets and one flow generator. Clearly we can write ad-hoc code to simulate any such system, but we are motivated to phrase our problem in matrix form, going forward.