MATH 247 - Numerical solutions of Differential Equations

January 21

Topics

Floating-point numbers.

• Introduction to the course.
• Sources of approximation in numerical simulations.
• Absolute and relative error.
• Data error and computational error.
• Floating Point numbers
• Truncation and Rounding errors.
• Examples of catastrophic cancellations
• Floating Point Arithmetic.
• Propagation of errors in floating point arithmetic.

Readings

• What is Numerical Analysis?, an interesting a short essay on the essence of this course's topic by one of its greater exponents, Nick Trefethen.
• Tutorials on on the positional notation:
  1. The positional notation from Wikipedia.
  2. The positional notation by D. Knuth.
• Floating point numbers:
  1. Section 1.7 (p. 36), from Chapter 1 of Moler's Numerical Computing with MATLAB beautiful (and free!) book
  2. Floating Point numbers
  3. IEEE 754 Floating Point Representation
  4. My slides.

Links

• Converter from decimal to single precision floating point
• Single precision floating point numbers from Wikipedia.
• Double precision floating point numbers from Wikipedia.
• Fall96Cleve.pdf, a very interesting article on Floating point systems by Cleve Moler
• The binary system from Wikipedia.
• The harmonic series from Wikipedia.
• This link is not strictly related to the lecture but, if you like Math, it will give you a quite wider and richer view over the harmonic series: The Riemann Zeta function from Wikipedia.
• Some disasters caused by numerical errors, a demonstration of how important is to deal carefully with floating point calculations!
• CVBK.pdf, a very interesting study of a remarkable example of catastrophic cancellation. It is a research article and it is not very easy to read but I suggest you to give it a try because it shows you in good detail how different hardwares work (e.g. Intel, Sun & IBM architectures).

Code

• em.m, a simple example of counter-intuitive floating point result.
• discrderiv.m, an example of truncation and roundoff errors: dependence on h of the evaluation of f'(x) through (f(x+h)-f(x))/h.
• catCanc1.m, first simple example of catastrophic cancelation.
• catCanc2.m, second simple example of catastrophic cancelation.
• float.m, visualization of three digits floating point systems via MatLab.

Homework

Due date: Feb. 5th.

1. Consider a Floating point system with base 10 and 3 digits mantissa (namely significant digits) and evaluate 3(4/3-1)-1 in that system.
2. Do the same for single and double precision standard floating point numbers.
3. Consider the problem of computing the Euclidean norm $\sqrt{a^2+b^2}$ of the vector $(a,b)$. What happens when $a$ and $b$ become "too big"? And what does it mean exactly "too big" in this context assuming to use 64bit floating point numbers? Suggest a way to re-write the formula so that the problem disappears. Write a program in MatLab or Python to evaluate this number in the "naif" way and in the way you suggest and verify that your formula work even when the first does not.
4. When is the expression $\sqrt{x^4+1}-x^2$ subject to catastrophic cancellation? Find a way to re-write the expression that avoids this problem. Write a program in MatLab or Python to evaluate this expression and your own one in a case that shows catastrophic cancellation and verify that your one gives, instead, the correct result.
5. Using either Python or MatLab or Octave (or all of the above!) do the following: define a variable $a$ to be equal to 2.6 and print it out to 15 digits after the dot. Then add to it 0.2 and print it again to verify the correctedness of the operation. After adding 0.2 to $a$ three times you will see the appearance of an error, explain why. What happend if, instead, you add directly 0.6 to 2.6? Why?
6. Consider the quadratic equation $x^2+54.32x+0.1=0$ and evaluate its roots to 11 digits via the standard solution formula with your favorite language. Now, evaluate the same roots using a four digits decimal floating point (namely, after every operation you perform, truncate the result to the first four significant digits) and evaluate the relative error for the two roots in this approximation. Find the reason for the result you find. Suggest an alternate formula to use that produces better results.
7. Write for the second derivative a code using the approximation $f''(x)\simeq(f(x+h)-2f(x)+f(x-h))/h^2$. Find theoretically for which order of magnitude in h the approximation will be the best and compare your result with the numerical data obtained by your code.

January 28

Topics

Coding in MATLAB and Python.

• A first example of program on MATLAB: em.m
• Basics of the MATLAB language -- arrays, matrices, loops
• A second example of program on MATLAB: allocation.m
• Using MATLAB interactively from the CLI.
• Conditionals in MATLAB.
• Basics of the python language, cycles in python.
• Two examples of python coding: em.py, emd.py.
• Plotting in MATLAB and Python.
• Conditionals in Python. Example: the Eratosthene's Sieve in python (eratosthenes.py) and in matlab (eratosthenes.m).

Readings

MATLAB/Octave:

• Section 1.7 (p. 36), from Chapter 1 of Moler's Numerical Computing with MATLAB beautiful (and free!) book
• The Basic Intro to MATLAB from U. of Oxford (UK)

Python:

• The chapter "An introduction to Pyhton" from the Python tutorial for a first course in Numerical Differential Equations at Brown U.
• The first 7 chapters (they are not that long!) from Think Python 2e by Allen B. Downey

Links

Just look at the languages page!

Problems

1. Using either Python or MatLab or Octave (or all of the above!) do the following: define a variable $a$ to be equal to 2.6 and print it out to 15 digits after the dot. Then add to it 0.2 and print it again to verify the correctedness of the operation. After adding 0.2 to $a$ three times you will see the appearance of an error, explain why. What happend if, instead, you add directly 0.6 to 2.6? Why?
2. Write code in MatLab or Python (or both!) to compute $r=p^2-2q^2$ with $p=665857$ and $q=470832$ in single and in double precision. Explain the reason for the unexpected result you get. Is any of those two numbers correct?
   Remark: in MatLab you can force a variable to represent a single-precision floating point number using the instruction single, e.g. p=single(665857). In Python there is no such thing and you will have to call NumPy, e.g. with import numpy. After that, you can use instrutions like p=numpy.float32(665857) to force a number to be stored in 32 bits format.
3. Repeat the previous exercise with $p=10864$, $q=18817$ and $r=9p^4-q^4+2q^2$.
4. Repeat the previous exercise with $p=10^{34}$, $q=-2$ and $r=p+q-p$.
5. Use Python to evaluate $e^{-5.5}$ using the Taylor series for $e^x$ in a floating point system with base 10 and 2 significant digits. You will have to consider many terms before the sum stabilizes but it will help you getting a good feeling of how a floating point system works. Then repeat the calculation in case of $\frac{1}{e^{5.5}}$ and compare the two results with each other and with the exact value.
   Remark: an easy way to avoid doing these calculations completely by hand is using Python's 'decimal' package (many thanks to Bishwa Silwal for pointing out this to me!).
   Here is an example of how to use it:
   >>> from decimal import * # Load the decimal package
   >>> getcontext().prec = 2 # Set the number of significant digits to 2
   >>> x=Decimal(-5.5) # Declare x as a 2-digits decimal fp and set it to -5.5
   >>> x*x
   Decimal('30')
   >>> 1+x+x*x/2+x*x*x/6
   Decimal('-17')

February 4

Topics

Embedding MATLAB/Python in a Web page.

Use SageMathCell! You can embed both Octave/MATLAB and Python code.

Coding in MATLAB and Python.

1. Using MATLAB and Python to find prime numbers with Eratosthenes' sieve:
   • eratosthenes.py,
   • eratosthenes.py.
2. Plotting in MATLAB and Python.

Elementary Methods to solve ODE's IVP.

1. One step methods:
   • Explicit Euler method:
     • explicitEuler.py,
     • explicitEuler.m.
   • Implicit Euler method:
     • implicitEuler.py,
     • implicitEuler.m.
2. Local and global error of the methods above.

Readings

MATLAB/Octave:

• Plotting with MATLAB
• Colors in MATLAB plots

Python:

• Graph plotting in Python
from tutorialspoint.
• Graph plotting with Python
from geeksforgeeks.

ODEs:

• Start reading Initial Value Problems. This is the second chapter of a nice set of slides that we will follow closely throughout the semester.
• Start reading Numerical solution of ordinary differential equations. This is a very interesting introductory survey on numerical methods for ODEs by E. Hairer and his former PhD student C. Lubich. Hairer is one of the main names in the theory of Geometric Integrators, namely numerical methods for ODEs.

Links

• If you forgot about Taylor series, read this
• Examples of Python code for the Explicit Euler method

February 11

Topics

1. Writing Python and MATLAB code in class to solve ODEs
2. Verifying numerically the order of convergence of the explicit Euler method

Readings

ODE slides

Links

Order and stability of explicit and implicit Euler methods

Problems

• Task 1
• Task 2
• Task 3
• Task 4

February 18

Topics

Elementary Methods to solve ODE's IVP.

1. Heun method (explicit trapezoidal rule):
   • heun.py,
   • heun.m.
2. Local and global error of the Heun method.
3. Methods for stiff ODEs:
   • explicitEuler-stiff.py,
   • implicitEuler-stiff.py.
4. Rounding errors and Compensated sums

Readings

Links

Examples of Python code for the Heun method

Problems

• Task 1
• Task 2
• Task 3
• Task 4

February 25

Topics

Euler methods in $\Bbb R^n$.

1. Explicit Euler
2. Implicit Euler
3. Symplectic Euler

March 3

Topics

Runge Kutta Methods.

1. Midpoint method
2. Trapezoidal method
3. RK4
Variable step size.

March 10

Postponed to after the break

March 17

Spring break

March 24

Topics

Boundary Value Problems.

1. Shooting Method
2. Finite Difference Method

March 31

Topics

Boundary Value Problems.

1. Shooting Method
2. Finite Difference Method

April 7

Topics

Hyperbolic PDEs.

1. The Advection PDE

Readings

1. The shooting method for solving BVPs by T. Lakoba (U. of Vermont).
2. The Finite Difference method for solving BVPs by G. Solderling (Lund University).
3. Boundary Value Problems by J. Blanchard (U. of Winsconsin)
4. A nice worked out example of the shooting method in python
5. A nice worked out example of the shooting method in MATLAB

References

1. Programming for Computations - A Gentle Introduction to Numerical Simulations with Python by H.P. Langtangen & S. Linge (Oslo U., Norway)

Code

1. shooting-eEuler.py
2. shooting-eEuler-improved.py
3. finiteDifferences.py

April 14

Topics

Hyperbolic PDEs.

1. The Advection PDE
2. The Method of Characteristics
3. The FTBS method
4. The FTCS method
5. The Lax method

Convergence, Stability and Consistence.

1. The CFL Condition
2. The Von Neumann Condition
3. The Lax Theorem

Readings

1. Solving the Advection PDE with Finite Differences in MATLAB (see Chap. 2 of the pdf) by A. Bergara (University of the Basque Country, Spain)
2. The Advection PDE by S. Gurevich (U. of Muenster, Germany).
3. The Wave PDE by S. Gurevich (U. of Muenster, Germany).
4. Consistency, stability and convergence by R. Fedkiw (Stanford U.)

April 21

Topics

1. Animations in MATLAB and Python
2. Solving Parabolic and Elliptic PDEs with the Finite Differences Method
3. Compiling python code with cython
4. The theory behind the Finite Elements Method.
5. Elementary implementation of the Finite Elements Method for ODEs.

Readings

1. A nice discussion on the computational complexity of the Finite Differences method applied to the Heat equation by J. Demmel (Berkeley).
2. A nice discussion on the computational complexity of the Finite Differences method applied to the Poisson equation by J. Demmel (Berkeley).
3. From Finite Differences to Finite Elements by G. Solderling (Lund University).
4. A Finite Elements primer, a theoretical introduction to the FEM by D. Silvester (U. of Manchester).
5. Remarks around 50 lines of Matlab: short finite element implementation.
6. Finite Elements Method.

April 28

Topics

1. Elementary implementation of the Finite Elements Method for PDEs.
2. Elementary Spectral Methods in PDE

Readings

1. Finite Elements in 1D, an elementary and insightful implementation of the FEM for ODEs, by Niall Madden, NUI Galway (Ireland)
2. Finite Elements in 2D, an elementary and insightful implementation of the FEM for Poisson PDEs, by Niall Madden, NUI Galway (Ireland)
3. FEM, a detailed discussion on how to implement the FEM for the Poisson PDE by O. Mali (U. of Jyväskylä)
4. FEM2, implementation in MATLAB of the algorithms discussed above by O. Mali (U. of Jyväskylä)
5. Fourier Transforms and Fast Fourier Transforms
6. Numerical solution of KdV

Code

The following files are FEM's MATLABB implementation by O. Mali -- I modified the main so that it accepts meshes produced by distmesh.

1. main.m
2. basis_linear.m
3. RefElemQuad.m
4. AssembleGlobalMatrices.m
5. fe2d_Madden.m
6. fe2d_Persson.m