David Bloore's MATLAB Codes
Introduction
There are many codes here arranged roughly in order to coincide with class material. Each code starts with a brief description of its contents. You can download them individually, or download a zip file containing all of them.
I’m using a February 2011 MacBook Pro, with the 2.3 GHz quad i7, 8 GB of RAM, and 1 GB of video RAM. Some codes use all the power my computer has to offer. However, they can easily be modified by looking at the global variables set up in the beginning of the codes. In order to make things run faster, simply reduce parameters, for example the number of x,y pairs used for a function. Often I’ll set this up as a parameter directly, combined with upper and lower bounds for x. In other cases, there’s no way to speed it up, for example performing elimination on K2D9. But for the Fourier codes, the Interpolation codes, and especially the Two and Three Body motion codes, changing the number of evaluation points will make a big difference. However, you may wish to ask yourself if anything else should change, or not, when you do this.
To run the animated codes, first download them into your MATLAB document directory (that you set up during installation to store saved files). Then, execute them in MATLAB. I highlight them in the "Current Folder" window, and press F9 (for PC) or shift-F7 (for Mac). Enjoy!
The animated codes are configured to produce a window that fits on a 1440-by-900 screen. They will be fine on anything larger, but for a smaller screen you may have to make an adjustment. At the command line, type "figure", maximize the blank figure window, then type "get(gcf,'Position')". The vector answer will be what you need to paste into the codes where they set the figure size (set(gcf,'Position',[a b c d]);).
Chapter 1
Basic Boundary Conditions
These codes require an input N, for the number of mesh points. They all assume constant f. Can you change them for variable f? How would you change Fixed/Free to Free/Fixed? What would need to change to look at solutions to a Free/Free system? (And when would a given f for Free/Free be acceptable? Why?)
Chapter 2
Least Squares
This is a fairly straightforward code which shows how to write a code to accept a variety of inputs. The same code works for 2 points or 100 points, for a fit of degree 0 or degree 99.
Two Body Problem
These are simulations of Two Body motion. During each interval, constant acceleration is assumed. Can you explain what the error should be? Are there ways to reduce it? I wrote this code using an approximation even though the analytic solution exists—because I want to simulate Three Body motion, know what the error will be, and how to mitigate it. There exists no analytic solution for Three Body motion, therefore an approximation is required in order to work with such systems.
What growth matrix regime is this operating with? Is it a Backward Euler? Forward? Is it a Trapezoidal Method approximation?
Note: These codes can be computationally intense. Some of the “cases” are much more so than others. Cases 1, 2, 5, 6, 8, and 9 should run well on “most” computers.
Three Body Problem
I invite you to play with sets of initial conditions to produce stable systems. It's not trivial!
Chapter 3
Lagrangian Interpolants and B-Splines
These codes require a scalar value for “m”. Start with m=3. For the 3D codes start with an inital m=3, and a range r=4. That will produce 5 lines or 5 series of lines (4 + the inital).
Do you see any difference in the fit quality as m increases? What about increasing m for the Lagrangian Interpolant compared to the Cubic Spline?
If you’re interested in looking at Lagrangian Interpolant or Spline fits to other functions, try it out. All you need to do is change a couple of lines, and once you know how this process works, you can apply it to any function on any domain!
- Interp1_Interpolant_and_Spline
- Interp2_Interpolants_and_Errors
- Interp3_Interpolant_Errors_3D
- Interp4_Spline_Errors_3D
- Interp5_Interpolant_3D
- Interp6_Spline_3D
Kronecker
Kronecker Products can produce some very interesting effects. There are also rules one can construct for them, for example as in "Products_of_Kroneckers." Also, you can use Kronecker Products of Incidence Matrices to use a graph as a kind of function on another graph, to produce more sophisticated constructs.
Poison Solvers Using K2D
Try both the K2D and K2D sparse codes. Start with a low N, for example 7, and work your way up. If you put a number too large in for N, you may risk not being able to use your computer for a while—possibly a long while. So work your way up with the full matrix-based code, in increments of less than ten, until it starts to take longer than a few seconds. Then switch to the sparse codes and repeat, although the increments will be different.
You may ask, what about even values for N? Well, I incorporated a test to ensure that the code does not accept even N. Why did I do that? When would I not?
Are you getting the same values for u(1/2,1/2) with different values of n? Why do you think that is? Would you want to change that? How, and why?
For the animated codes, start with an initial N of 3, and choose 20 frames. This should calculate pretty quickly, and animate well. Once the system size starts to get large (for your machine) the time between frames will increase, so the above method for testing your computer's ability will give you an idea of where to stop this animation. Remember each frame increases N by 2, not 1!
Note: MATLAB stops displaying solutions to this system beyond a certain n. To overcome this, I scaled the solution up to a displayable size (only in two of the codes, but you can change that if you want). Relative proportions are not affected. Also, in the animations, MATLAB adjusts the limits of the z-axis in response to changes in Z, so please be aware of that when interpreting the meaning of what you see.
- Poisson1_K2D_Full
- Poisson2_K2D_Sparse
- Poisson3_K2D_Sparse_Mesh
- Poisson4_K2D_Animated_Scaled
- Poisson5_K2D_Mesh_Animated
- Poisson6_K2D_Functions_of_N
- Poisson7_K2D_Sparse_f_1
Poison Solvers Using K2D9
(See p.292, problem 5)
Again, start with low values/initial values for N (7 is good), and work your way up (in increments of 10, or 5 once things slow a little).
What do you think the eigenvectors and eigenvalues of K2D9 are? Would you expect a different solution to Poisson’s Equation using K2D9? How would it differ from K2D? Why?
- Poisson8_K2D9_Sparse
- Poisson9_K2D9_Sparse_Mesh
- Poisson10_K2D9_Animated_Scaled
- Poisson11_K2D9_Mesh_Animated
- Poisson12_K2D9_Functions_of_N
- Poisson13_K2D9_Sparse_f_1
Chapter 4
Fourier Series
These are a lot of fun. One question that arose in my mind was, why does the arc length of the Delta Function appear to → ∞ as the number of terms →∞? The delta function is defined as 0 in the interval [π/4,π/2], yet the arc length of the sum of the first n terms of the Fourier Series for the Delta Function is clearly not π/4! In fact, it appears to increase with n, so that the limit as n → ∞ is ∞! Check out case 2 in Delta_Function_2.
If the Fourier codes are running slowly, decrease the variable “s.” That is the number of points in the domain for which the function is evaluated.
Double Fourier Series
These can run slow on any computer. Start with 9 terms evaluated at 150 points per axis. I've run these as high as 300 terms on 2000 points per axis, but that took my computer about 2.5 hours. Can you establish a relation for number of operations as a function of the number of terms and the number of evaluation points?
- Poisson_Double_Sine
- Poisson_Double_Sine_Animated
- Poisson_Double_Sine_Animated_Heat2D
- Double_Delta
- Double_Delta_Animated
Other Codes
This code allows you to look at the accuracy of truncated Taylor Series for f(x)=e^x.
Dynamic Wave Functions
In Physics one often encounters waves, superposition of waves, or other things exhibiting either sinusoidal or periodic oscillation. The following are intended to give you access to these functions and their superpositions, so you can see how frequencies/wavelengths and velocities of superimposed translating waves interact.
One of the codes varies only wavelength, with no translation, yet it still produces quite a pattern!
- Travelling_1Wave
- Travelling_2Waves
- Travelling_3Waves
- Standing_Waves_Dynamic
- Travelling_Waves_Dynamic1
- Travelling_Waves_Dynamic2
- Travelling_Waves_Dynamic3