MATH 247
Numerical solutions of
Differential Equations

Instructor: Roberto De Leo.     Term: Spring 2022.

Lecture 1, January 18

Topics

  • Introduction to the course.
  • Numbers in base 2, 10 and 16.

Readings

Homework #1

Due date: Jan. 27
  1. Write the decimal numbers 65,17,1024 in binary (look at Sec. 3.2.2 of Numerical Algorithms and digital representation).
  2. Write in binary the following fractions: $\frac{1}{2},\,\frac{1}{3},\,\frac{1}{4},\,\frac{1}{5},\,\frac{1}{8},\,\frac{1}{10}$ (look at Sec. 3.3.3 of Numerical Algorithms and digital representation).

Lecture 2, January 20

Topics

  • Floating Point systems
  • Single and Double precision floating point numbers

Readings

  1. Section 1.7 (it begins at the end of p. 35), from Chapter 1 of Moler's Numerical Computing with MATLAB beautiful (and free!) book
  2. The first 8 pages from these Floating Point numbers slides
  3. IEEE 754 Floating Point Representation
  4. My slides.

Useful Links

Homework #1

Due date: Jan. 27
  1. Consider a Floating point system F with base 10 and 3 digits mantissa (namely each number in F has 3 digits) and evaluate the expression $$3(4/3-1)-1$$ in that system.
    Remark: when you do this, you need to do operations one by one starting from the innermost and moving outwards. So, first evaluate 4/3, then subtract to it 1 and so on. At every step you need to cut the result to make it a number of F.
  2. Evaluate the same expression for single (32bits) and double precision (64bits) floating point numbers.
  3. Set $a=10^{35}, b=10^{-5}, c=-10^{35}, d=10^{-5}$. Evaluate by hand the results, in double precision, of the operations $a+b+c+d$ (which means $((a+b)+c)+d$) and $(a+b)+(c+d)$. Verify your guess with MATLAB or Python.
    Remark: and easy way to input, say, $10^{35}$ into a variable a in MATLAB and Python is using the syntax
    a=1e35;

Lecture 3, January 25

Topics

  • Floating Point numbers
  • Floating Point Arithmetic.
  • Propagation of errors in floating point arithmetic.
  • Examples of catastrophic cancellations

Readings

In case you want to know more

Links

  • 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).

Homework #2

Due date: Feb. 3

Solve problems 2,3,7 in 1.5.

Lecture 4, January 27

Topics

(purple topics are about coding in MATLAB/Octave and Python)
  • Propagation of errors in floating point arithmetic.
  • An example of catastrophic cancellations: solving $x^2-2x+\delta=0$.
  • The four operations.
  • Loops and Conditionals.
  • What is Numerical Analysis.
  • Numerical evaluation of the first derivative of a function.
  • Truncation and Rounding errors.

Readings

Math MATLAB Python

In case you want to know more

Homework #2

Due date: Feb. 3

  1. When do you expect that the expression $$\sqrt{x^4+1}-x^2$$ is 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 version in a case that shows catastrophic cancellation and verify that your one gives a much more precise result than the original one.
  2. 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 works even when the first does not.
  3. Write MATLAB and Python codes to evaluate the sum of the cubes of the first 20 natural integers and verify the result using the formula $$ \sum_{k=1}^n k^3 = \frac{n^2(n+1)^2}{4}. $$
Extra credit: explain how is the algorithm in sqrt2.html able to find the digits of $\sqrt{2}$. (This problem is not trivial, I will keep assigning it as extra credit in homework until someone will come up with a solution.)

Lecture 5, February 1

Topics

  • Truncation and round-off errors
  • Numerical evaluation of the first derivative of a function: Forward differencing.
  • Numerical evaluation of the first derivative of a function: Centered differencing.
  • Arrays in MATLAB.
  • Plots in MATAB.

Readings

Math MATLAB Python

Homework #3

Due date: Feb. 10

  1. Using either Python or MatLab/Octave (or both!) 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 happens 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. Modify the MATLAB/Octave code in Section 2.2 to evaluate the precision of the Forward differencing method (2.2.1) and of the Centered differencing method (2.2.2) in case of the function $e^x$ at $x_0=2$.
    Remark: in both MatLab and Python the value of $e^x$ is obtain using the instruction exp(x).
Extra credit: explain how is the algorithm in sqrt2.html able to find the digits of $\sqrt{2}$. (This problem is not trivial, I will keep assigning it as extra credit in homework until someone will come up with a solution.)

Lecture 6, February 3

Topics

  • Arrays in Python.
  • Plots in Python.
  • Ordinary Differential Equations.
  • ODEs Initial Value Problems: Existence and Uniqueness of solutions.

Readings

Math Python

Lecture 7, February 8

Topics

  • Examples of IVP with:
    • Unique solutions defined at all times;
    • Unique solutiosn that blow up at finite time;
    • Non-unique solutions.
  • Explicit Euler method:
    • Local error.
    • An example of implementation (from Sec. 9.3).

Readings

Math Python

Lecture 8, February 10

Topics

  • Explicit Euler method:
    • A second example of implementation.
    • Global error.
    • Round-off error.
    • Compensated summation.
    • Stability.
    • Variable step-size.

Readings

Homework set #4

Due date: Feb. 17

  1. Write (by hand, no numerics involved!) any two distinct solutions of $\dot x=x^{1/3}, x(0)=0$.
  2. Find (by hand, no numerics involved!) the unique solution to $\dot x=x^3, x(x_0)=t_0$ and find at which time does the function blows up.
  3. Modify the code at Sec. 1.3 (or write your own!) to evaluate the solution to the IVP $$\dot x = x^2,\;\; x(1)=1$$ at time $t=1.5$ and find which value of $h=10^{-k}$ is needed to find $x(1.5)$ with a precision of $10^{-6}$.
  4. Use some code at Sec. 1.3 (or write your own!) to show numerically that the local error of the Euler method is of order 2.
Extra credit: explain how is the algorithm in sqrt2.html able to find the digits of $\sqrt{2}$. (This problem is not trivial, I will keep assigning it as extra credit in homework until someone will come up with a solution.)

Lecture 9, February 15

Topics

  • Explicit Euler method:
    • Compensated Summation
    • Variable step-size.
  • Explicit Heun method.

Readings

In case you want more readings

Lecture 10, February 17

Topics

  • Taylor methods
  • Multistep methods
  • Runge-Kutta methods

Readings

Homework set #5

Due date: Feb. 24

  1. In class we used all algorithms I covered to solve ODEs whose exact solution is known in order to check their convergence. In real life, of course we have no such luxury. How can we say then that a numerical solution is correct? A criterion often used in practice is the following heuristic principle: the digits that stabilize as $h$ decreases are correct. A second criterion is solving the same ODE with a different method and verify that one gets the same values with both of them.

    Adopt these principles to find the first 6 digits of the value at $t=10$ of the solution of the following IVP, for which the exact expression is not known: $$ \dot x = \sin t-0.1x^3,\;\;x(0)=4. $$ Use the following numerical methods: Euler, RK4, the Taylor method of order 2 and the Adam-Bashcroft multistep method of order 2 (either write your own code or modify the one in the lecture notes). For each method, find the order of magnitude of the smallest h for which the first 6 digits of $x(10)$ are correct. Send me the code you used in a text attachment.
Extra credit: explain how is the algorithm in sqrt2.html able to find the digits of $\sqrt{2}$. (This problem is not trivial, I will keep assigning it as extra credit in homework until someone will come up with a solution.)

Lecture 11, February 22

Topics

  • Convergence, Consistency and Stability of an ODE method.
  • Implicit Methods
  • Stability regions.

Readings

Lecture 12, February 24

Topics

  • Stiff ODEs.
  • Examples of IVP on the plane.

Readings

  • Sections 1.11, 1.12.
  • Sections 1 and 2 from the Mini-course on Geometric Integrators by D. Cohen (Umea U., Sweden)

Homework set #6

Due date: Mar. 3

  1. Prove that the Implicit Euler method is consistent, namely that its local error goes to 0 with h.
  2. Prove the formula for the stability region of the RK4 method.
  3. Find the formula for the stability region of the implicit trapezoidal algorithm.
Extra credit: explain how is the algorithm in sqrt2.html able to find the digits of $\sqrt{2}$. (This problem is not trivial, I will keep assigning it as extra credit in homework until someone will come up with a solution.)

Lecture 13, March 1

Topics

  • Examples of IVP on the plane.
  • Geometrical integrators.
  • SciPy IVP solvers
  • MATLAB IVP solvers

Readings

  • Sections 1.12, 1.13 and 1.14
  • Sections 1 and 2 from the Mini-course on Geometric Integrators by D. Cohen (Umea U., Sweden)

Lecture 14, March 3

Topics

  • Boundary Value problems:
    • Shooting method;
    • Finite Differences method.

Readings

Lecture 15, March 15

Topics

  • SciPy BVP solvers
  • MATLAB BVP solvers
  • General facts about PDEs

Readings

Lecture 16, March 17

Topics

  • Classification of PDEs based on characteristics.
  • The Advection PDE.
  • The FTBS method.

Readings

Homework set #7

Due date: March 24

  1. Recall that the energy of the physical pendulum ODE in Sec. 1.12 is given by $$E(x,y) = \frac{y^2}{2}+1-\cos x.$$ Write code (or modify the existing code) to plot the difference between the exact energy and the numerical energy versus time for the case solved in the lecture notes with the three methods: Explicit Euler, Implicit Euler and Symplectic Euler. Leave all parameter values unchanged. What do these plots say to you?
  2. Solve numerically the linear BVP in Sec. 1.17 for n=1250,2500,5000,10000,20000. For each of these values of n, evaluate the maximum error between the numerical and the exact solution and use these values to estimate the global order of the Finite Difference method.
Extra credit: explain how is the algorithm in sqrt2.html able to find the digits of $\sqrt{2}$. (This problem is not trivial, I will keep assigning it as extra credit in homework until someone will come up with a solution.)

Lecture 17, March 22

Topics

  • FTBS: local error order
  • FTBS: the CFL condition
  • FTBS: von Neumann conditions
  • FTBS: Python implementation

Readings

Lecture 18, March 24

Topics

  • Different types out of ways of displaying the numerical solutions in Python
  • FTFS method
  • FTCS method
  • Lax method
  • Why do we see dispersion in advection's numerical solutions

Readings

Lecture 19, March 29

Topics

  • FTBS global order
  • Code vectorization
  • 1D animations in Python
  • Leapfrog method

Readings

Readings

Lecture 20, March 31

Topics

  • Lax-Wendroff method
  • The 1D Wave PDE
  • FTBS method
  • Lax method
  • Leapfrog method
  • CTCS method
  • The 2D Wave PDE

Readings

Lecture 21, Apr 5

Topics

  • 2D animations in Python
  • Consistency, Stability and Convergence of a Finite Difference method
  • The Lax Theorem: under consistency, Stability and Convergence are equivalent
  • Parabolic PDEs
  • FTCS method

Readings

Lecture 22, Apr 7

Topics

  • Analytical solutions of the Heat PDE
  • Instantaneous propagation of signals under Parabolic PDEs
  • CTCS method
  • DuFort-Frankel method
  • BTCS method
  • Crank-Nicholson method
  • The 2D Heat PDE
  • The 1D Advection-Diffusion PDE

Readings

  • Section 2.4
  • Watch this MIT video lecture about the difference in speed propagation of signals in Wave and Heat PDEs.

Homework set #8

Due date: April 16

  1. Modify the FTBS method in Section 2.2.1.1 to study numerically the Burgers' nonlinear PDE $$ u_t+u\,u_x=0. $$ Send me your code and the graphs you obtain. You will see that, unlike the linear advection PDE, different parts of the graph move with different speeds. Explain why does this happen by comparing Burgers' PDE with the advection PDE.
  2. Do the same as above with the PDE $$ u_t+\frac{u_x}{1+u}=0. $$
Extra credit: explain how is the algorithm in sqrt2.html able to find the digits of $\sqrt{2}$. (This problem is not trivial, I will keep assigning it as extra credit in homework until someone will come up with a solution.)

Lecture 23, Apr 12

Topics

  • Elliptic PDEs

Readings

Extra readings

Lecture 24, Apr 14

Topics

  • The Finite Element Method: 1D

Readings

Lecture 25, Apr 19

Topics

Readings

Lecture 26, Apr 21

Topics

  • The Finite Element Method: 2D
  • Review

Readings