Fall 2023: Numerical Analysis (Math 436/565)

Overview

Lecture sketches and homework

To help you study for this class, I am posting lecture sketches. These are not a substitute for your personal lecture notes or coming to class (for instance, lots of details and motivation are not included in the sketches). I hope that they are useful to you for revisiting the material and for preparing for exams.

About the homework

• Homework problems are posted for each unit. Homework is submitted online, and you have an unlimited number of attempts. Only the best score is used for your grade.

  Most problems have a random component (which allows you to continue practicing throughout the semester without putting your scores at risk).

• Aim to complete the problems well before the posted due date.

  A 15% penalty applies if homework is submitted late.

• Collect a bonus point for each mathematical typo you find in the lecture notes (that is not yet fixed online), or by reporting mistakes in the homework system. Each bonus point is worth 1% towards a midterm exam.

  The homework system is written by myself in the hope that you find it beneficial. Please help make it as useful as possible by letting me know about any issues!

Exams and practice material

The following material will help you prepare for the exams.

• Midterm Exam 1:
  midterm01-practice.pdf, midterm01-practice-solution.pdf
• Midterm Exam 2:
  midterm02-practice.pdf, midterm02-practice-solution.pdf
• Final Exam:
  final-practice.pdf, final-practice-solution.pdf

Python

If you just want to run a handful quick computations (without saving your work), you can use the text box below. However, keep in mind that the code entered in this box is run as a Python script. That means that we need to use print(...) if we want to see any output.

One convenient option for running Python code in the cloud is Colab by Google Research. You can sign in with any Google account (since Google is managing university email, you can use your university account).

Perform your first computation by entering something like "1+2" and pressing Shift+Enter.

Here are some other things to try:

• The following code plots sine versus a quadratic interpolation:

  from numpy import linspace, pi, sin
  from scipy import interpolate
  import matplotlib.pyplot as plt
  xpoints = [0, pi/2, pi]
  ypoints = [sin(x) for x in xpoints]
  poly = interpolate.lagrange(xpoints, ypoints)
  xplot = linspace(0, pi, 100)
  plt.plot(xpoints, ypoints, 'o', xplot, sin(xplot), '-', xplot, poly(xplot), ':')
  plt.show()
  

• Interpolation using equally spaced points versus Chebyshev nodes:

  from numpy import linspace, pi, cos
  from scipy import interpolate
  import matplotlib.pyplot as plt
  
  def f(x):
      return 1/(1+25*x**2)
  
  n = 5
  
  xpoints = linspace(-1, 1, n)
  ypoints = [f(x) for x in xpoints]
  poly = interpolate.lagrange(xpoints, ypoints)
  
  xpoints2 = [cos((2*j+1)*pi/(2*n)) for j in range(n)]
  ypoints2 = [f(x) for x in xpoints2]
  poly2 = interpolate.lagrange(xpoints2, ypoints2)
  
  xplot = linspace(-1, 1, 100)
  plt.plot(xpoints, ypoints, 'o', xplot, poly(xplot), ':')
  plt.plot(xpoints2, ypoints2, 'o', xplot, poly2(xplot), ':')
  plt.plot(xplot, f(xplot), '-')
  plt.show()
  

• We can construct a (natural) cubic spline and plot it together with the knots:

  from numpy import linspace
  from scipy import interpolate
  import matplotlib.pyplot as plt
  xpoints = [1, 2, 4, 5, 7]
  ypoints = [2, 1, 4, 3, 2]
  spline = interpolate.CubicSpline(xpoints, ypoints, bc_type='natural')
  xplot = linspace(1, 7, 100)
  plt.plot(xplot, spline(xplot), '-', label='spline (natural)')
  plt.plot(xpoints, ypoints, 'o', label='knots')
  plt.legend()
  plt.show()
  

  Other standard choices for the boundary conditions bc_type include 'not-a-knot' (the default) as well as 'clamped' and 'periodic' (this one requires the first and last point to have the same y-coordinates).
• The following compares forward differences (order 1), central differences (order 2) and the Richardson extrapolation of the latter (order 4) for approximating the derivative of $ f'(1) = 2 \ln(2) $ where $ f(x) = 2^x $. Observe how all approximations turn sour if $ h = 10^{-n} $ is too small but how we are able to get a better "best approximation" using methods of higher order:

  from math import log
  def forward_difference(f, x, h):
      return (f(x+h)-f(x))/h
  def central_difference(f, x, h):
      return (f(x+h)-f(x-h))/(2*h)
  def central_difference_richardson(f, x, h):
      return (-f(x+2*h)+8*f(x+h)-8*f(x-h)+f(x-2*h))/(12*h)
  def f(x):
      return 2**x
  print([forward_difference(f, 1, 10**-n) - 2*log(2) for n in range(12)])
  print([central_difference(f, 1, 10**-n) - 2*log(2) for n in range(12)])
  print([central_difference_richardson(f, 1, 10**-n) - 2*log(2) for n in range(12)])
  

• The following implements the trapezoidal rule and applies it for approximating the integral of $ 1/x $ from $ 1 $ to $ 3 $ (which is $ \log(3) $). Observe how the errors confirm that the trapezoidal rule has order 2.

  from math import log
  def trapezoidal_rule(f, a, b, n):
      h = (b - a) / n
      integral = f(a) + f(b)
      for i in range(1,n):
          integral += 2*f(a + i*h)
      return h/2*integral
  def f(x):
      return 1/x
  print([trapezoidal_rule(f, 1, 3, 10**n) - log(3) for n in range(1,6)])
  

• In the following, we apply a Taylor method of order 3 to solve the differential equation $ y' = \cos(x) y $ with $ y(0) = 1 $. For comparison, we also plot the exact solution $ y(x) = e^{\sin(x)} $.

  from numpy import linspace, exp, cos, sin
  import matplotlib.pyplot as plt
  
  def taylor_3_cosy(x0, y0, xmax, n):
      h = (xmax - x0) / n
      ypoints = [y0]
      for i in range(n):
          y0 = y0 + cos(x0)*y0*h + 1/2*(cos(x0)**2-sin(x0))*y0*h**2 + \
               1/6*(cos(x0)**2-3*sin(x0)-1)*cos(x0)*y0*h**3
          x0 = x0 + h
          ypoints.append(y0)
      return ypoints
  
  xpoints = linspace(0, 2, 100)
  plt.plot(xpoints, exp(sin(xpoints)), '-', label='y(x)')
  
  nn = [2, 4, 8]
  for n in nn:
      xpoints = linspace(0, 2, n+1)
      ypoints = taylor_3_cosy(0, 1, 2, n)
      plt.plot(xpoints, ypoints, ':o', label='n = %s' % n)
  
  plt.legend()
  plt.show()
  

• The following illustrates that the midpoint method converges with a global error of 2:

  from math import e, cos, sin
  def midpoint(f, x0, y0, xmax, n):
      h = (xmax - x0) / n
      ypoints = [y0]
      for i in range(n):
          y0 = y0 + f(x0+h/2,y0+f(x0,y0)*h/2)*h
          x0 = x0 + h
          ypoints.append(y0)
      return ypoints
  def f_cosx_y(x, y):
      return cos(x)*y
  print([midpoint(f_cosx_y, 0, 1, 2, 10**n)[-1] - e**sin(2) for n in range(6)])
  

Projects

If you take this class for graduate credit, you need to complete a project. The idea is to gain additional insight into a topic that you are particularly interested in. Some suggestions for projects are listed further below.

• The outcome of the project should be a short paper (about 5 pages)
  • in which you introduce the topic, and then
  • describe how you explored the topic.
  Here, exploring can mean (but is not limited to)
  • computations or visualizations you did in, preferably, Python,
  • working out representative examples, or
  • combining different sources to get an overall picture.

Each project should have either a computational part (this is a great chance to learn more Python!) or have a more mathematical component. Here are some ideas:

• Read the recent article on Newton's Method Without Division by Jeffrey D. Blanchard and Marc Chamberland. Explain the main result and implement several examples to illustrate it.
• Describe and implement the Illinois method for computing roots.
• Introduce and explore Chebyshev polynomials. For instance, prove that their roots provide the optimal nodes for minimizing the generic error bound for polynomial interpolation.
• Introduce orthogonal polynomials and explore how they are used in solving differential equations or for approximating other functions.
• Explore the discrete Fourier transform.
• Describe and implement a fast multiplication algorithm such as the Karatsuba algorithm.