Summarizing multiple-choice questions
Quiz
Exercise 6.1: Characterize a finite difference
Question: We can approximate the derivative at a point using two function values:
What is this type of difference called and how large is the error?
Choice A: This is a centered difference with error proportional to \( \Delta t^2 \).
No, a centered difference would have the two function values equally displaced to either side of the target point \( t_n \).
Choice B: This is a forward difference with error proportional to \( \Delta t \).
The name is forward difference, or Forward Euler difference in the context of differential equations. In such contexts the formula is also known as Euler's method or formula.
Choice C: This is a Forward Euler difference with error proportional to \( \Delta t^3 \).
One may well call this a Forward Euler difference, but the error is not as "good" as \( \Delta t^3 \).
Choice D: This is a Backward Euler finite difference with error proportional to \( \Delta t \).
Since we use the points \( t_{n+1} \) and \( t_n \) when constructing the difference, we go forward in time, not backward. Therefore, this is a forward difference.
Exercise 6.2: Characterize a finite difference
Question: We can approximate the derivative at a point using two function values:
What is this type of difference called and how large is the error?
Choice A: This is a centered difference with error proportional to \( \Delta t^2 \).
Correct!
Choice B: This is a forward difference with error proportional to \( \Delta t^2 \).
A forward difference makes use of the point itself, \( t_n \), and a point forward in time, \( t_{n+1} \). This is not the case here: we use points to the left and right, and the derivative is approximated in the center point. Also, a forward difference would not have an error \( \mathcal{O}(\Delta t^2) \).
Choice C: This is a centered difference with error proportional to \( \Delta t^4 \).
It is centered, but the error is only \( \mathcal{O}(\Delta t^2) \).
Choice D: This is a Backward Euler finite difference with error proportional to \( \Delta t \).
Since we use the points \( t_{n+1} \) and \( t_n \) when constructing the difference at \( t_{n+1/2} \), the derivative is in the center of the two points, and the difference is therefore a centered difference. A backward difference, or Backward Euler difference, would use \( t_{n+1/2} \) in this case and some point \( t_{n-1/2} \) backward in time.
Exercise 6.3: What is the problem with this program?
Question: We want to solve $$ u ' = -au, \quad u(0)=1, $$ by a Forward Euler scheme, $$ u^{n+1} = u^n - a\Delta t\, u^n,$$ and the following program
from numpy import *
u = zeros(10)
u[0] = 1
dt = 0.1
for i in range(10):
u[i+1] = u[i] - dt*a*u[n]
What is the major problem with this program?
Choice A:
from numpy import * is not recommended; one should import explicitly the functions needed or do import numpy as np.
True, these are recommended rules, but it is not a problem to do the "star import": it is legal and convenient.
Choice B:
The program aborts with a NameError.
True, a is not defined. This is the only reason why one cannot execute the program.
Choice C:
The program is "flat". It should be wrapped in a function solver(U0, dt, N, a).
That is definitely a good idea, but it is not a major problem for computing the solution of the differential equation.
Choice D: The scheme is unstable!
Actually, this can be true. Since a is not defined, it may happen that dt=0.1 is a too long time step for the stability restrictions of the Forward Euler method for this differential equation. We need \( \Delta t\leq 1/a \). However, the answer is wrong in the sense that this potential instability is not the major problem with the program - it is a bigger problem that a is not defined.
Exercise 6.4: Is the solution correct?
Question: You have solved \( u'=-au \) by the Crank-Nicolson scheme and tested your program thoroughly. Suddenly you run the program with \( I=1 \), \( a=10 \), \( T=1 \), and \( \Delta t=0.225 \). You get somewhat unexpected results:
The results are unexpected because we know the exact solution should be monotone and decreasing, while this numerical solution also shows an increasing stage. What is the problem?
Choice A: The program is not tested well enough - yet another bug is there.
No, this computation is in fact correct - numerically.
Choice B: The numerical solution method is more general than the analytical one, and this is an example where the solution can increase, but the analytical solution technique is not capable of dealing with this situation.
The analytical solution \( u(t)=Ie^{-at} \) covers all possible cases.
Choice C: The time step is too large and cause an instability in the form of non-physical oscillations.
True. One needs \( \Delta t \leq 2/a = 0.2 \) in this case to avoid oscillations.
Choice D: The Crank-Nicolson scheme is unconditionally stable and can be used with all time-step sizes, but the problem here is that round-off errors due to a "too large" \( a \) (10, not around 1) accumulate to the effect seen in the plot.
1: The first part is true if unconditionally stable means that the solution is bounded and decays with time, but one may also argue that oscillations, which are non-physical, are a kind of instability, and these occur if \( \Delta t >2/a \).
2: Round-off errors are very small in this problem, compared to the discretization errors, and cannot cause an oscillating solution.
Exercise 6.5: Is this a proper test function?
Question: Suppose we have some function compute that we want to test. We construct
a unit test and implement an associated test function (according to the
rules for test functions in the nose or pytest test frameworks):
def test_compute(n):
expected = (n**2)/4.5 - 1
computed = compute(n)
assert expected != computed
The question is if this test function can be used as is or if improvements must be implemented.
Choice A:
One should use special assert functions from nose.tools, here the choice is nose.tools.assert_almost_equal.
It is true that an "almost equal" type of assert function is appropriate, but there is nothing in the rules that requires use of special assert functions. The requirement is that an AssertionError is raised if the test fails. That is done by a plain assert as used here.
Choice B:
One cannot test compute(n) for only one n value. Many are required for good evidence that the function works.
This depends on what is inside compute. If it has several branches depending on the value of n, one must test for a visit to each branch, which requires multiple n values, but if it is a formula (the test might indicate so), one value can be sufficient.
Choice C:
The test function does not test expected != computed with a tolerance
and the n parameter cannot be an argument.
The formula for expected indicates that this is a real number that is subject to potential round-off errors, so one should use a tolerance: abs(expected - computed) < 1E-14. Also, test functions should never take arguments.
Choice D:
The assert statement also needs a message explaining what is wrong when the test fails.
This is always a good idea, but not a requirement.
Exercise 6.6: Rewrite an expression with array arithmetics
Question: A mesh function is initialized with the code segment
from math import exp
N_t = 100
T = 3.0
dt = T/N_t
t = []
u = []
for i in range(len(t)):
t.append(i*dt)
u.append(1 - exp(-t[-1]))
Rewrite this code such that there is no loop and u is computed by
array arithmetics.
from numpy import *
from math import exp
N_t = 100
t = np.linspace(0, 3, N_t+1)
u = 1 - exp(-t)
The exp function is imported from math (since from math import exp
reimports the exp name that was imported by from numpy import *)
and cannot be used with array argument t.
import numpy as np
N_t = 100
t = np.linspace(0, 3, N_t+1)
u = 1 - np.exp(-t)
Correct!
N_t = 100
dt = 3.0/N_t
t = [i*dt for i in range(N_t+1)]
from numpy import exp
u = 1 - exp(-t)
The computation of u applies array arithmetics, but the list
comprehension for t involves a standard (slow) Python loop.
import numpy
t = numpy.mesh([0, 3], 100)
u = numpy.func(numpy.exp, t)
There are no functions numpy.mesh and numpy.func.
Exercise 6.7: What is the truncation error?
Question: What is the truncation error?
Choice A: The difference between the exact solution and the numerical solution at a mesh point.
This is the global error.
Choice B:
The error in the factor that takes u[i] to u[i+1].
This is the amplification factor error or the local error.
Choice C: The difference between the exact solution and the numerical solution truncated to one decimal.
Nobody applies this error measure.
Choice D: The error in the scheme when the exact solution is inserted in the scheme's difference equation.
Correct!
Exercise 6.8: Recognize a programming language
Question: What kind of programming language is this?
function integral = trapezoidal(f, a, b, n)
%% Integrate f from a to b with n intervals
h = (b-a)/n;
result = 0.5*f(a) + 0.5*f(b);
for i = 1:(n-1)
result = result + f(a + i*h);
end
integral = h*result;
end
Wrong!
Correct!
Wrong!
Wrong!
Exercise 6.9: Recognize a programming language
Question: What kind of programming language is this?
def trapezoidal(f, a, b, n):
# Integrate f from a to b with n intervals
h = (b-a)/float(n)
result = 0.5*f(a) + 0.5*f(b)
for i in range(1, n):
result += f(a + i*h)
return h*result
Correct!
Wrong!
Wrong!
Wrong!
Exercise 6.10: Recognize a programming language
Question: What kind of programming language is this?
real*8 function trapezoidal(f, a, b, n)
C Integrate f from a to b with n intervals
real*8 h, result, f
external f
h = (b-a)/n
result = 0.5*f(a) + 0.5*f(b)
do i = 1, n-1
result = result + f(a + i*h)
end do
trapezoidal = h*result
return
end
Wrong!
Wrong!
Wrong!
Correct!
Exercise 6.11: Recognize a programming language
Question: What kind of programming language is this?
double trapezoidal(
double (*f)(double), double a, double b, int n)
{
double h, result;
h = (b-a)/n;
result = 0.5*f(a) + 0.5*f(b);
for (i=1; i++; i <= n-1) {
result += f(a + i*h);
}
return h*result;
}
Correct!
Wrong!
Wrong!
Wrong!
Exercise 6.12: What is SymPy?
Question: What is SymPy?
Choice A: A Python module for computing with symmetric matrices.
SymPy can compute with symmetric matrices, so that is true, but it can do very much more.
Choice B: A Python package for doing symbolic computations (exact/analytical differentiation, integration, equation solves, etc.).
Correct!
Choice C: A Python package for numerical approximations to differentiation, integration, equation solving, etc.
It's the opposite: analytical, not numerical.
Choice D: A free, open source version of Mathematica.
Could be viewed as such, but Mathematica is a much more advanced tool for symbolic computing.
Exercise 6.13: Testing of code
Question: What is an appropriate test for computing the solution of some ODE?
Choice A: We compare the numerical solution with the analytical solution in a plot. The curves are quite close, showing that the computations are correct (here using 31 mesh points).
This comparison says nothing if the discrepancy is the unavoidable
numerical error or if it also contains the effect of bugs in the program.
Choice B: The maximum error between the numerical and exact solution at the mesh points is \( 1.886\cdot 10^{-5} \), which is small. Therefore, the program works.
Nobody knows if this is the numerical approximation error or if it
also contains programming errors.
Choice C: The scheme is \( \Oof{\Delta t^2} \). With 31 points we get a maximum point-wise error of \( 1.243\cdot 10^{-4} \). With 61 points (halving \( \Delta t \)) the corresponding error is \( 3.108\cdot 10^{-5} \). This is a reduction of approximately a factor 4, which is what we expect for such a scheme when halving \( \Delta t \).
This test involves checking of the convergence rate, which is a good
test. The only knowledge we have of the numerical error is its rate!
Choice D: The numerical solution and the exact solutions get closer and closer in a plot as we reduce \( \Delta t \). Therefore, the program works.
This test points to convergence of the method, but the error can
still contain the effect of bugs. One needs to measure how fast
the curves get closer and closer.
Exercise 6.14: What kind of scheme is this?
Question: Given \( u^{\prime}=-au + b \), where \( a \) and \( b \) are functions of time. Is the following scheme correct? $$ \frac{u^{n}-u^{n-1}}{\Delta t} = -a(t_n)\frac{1}{2}(u^n + u^{n-1}) + b(t_n)\tp$$
Choice A: Yes, it is a Crank-Nicolson type of scheme.
No, then \( a \) and \( b \) should have been evaluated at \( t_{n-1/2} \).
Choice B: Yes, if \( a(t_n) \) and \( b(t_n) \) are evaluated at the previous time step, \( a(t_{n-1}) \) and \( a(t_{n-1}) \).
No, this will be wrong since the arithmetic mean on the right-hand
side points sampling the ODE at \( t_{n-1/2} \). Then \( a \) and \( b \) must
be evaluated at this point.
Choice C: No, the coefficients \( a \) and \( b \) are not evaluated correctly.
That is right: they should be evaluated at \( t_{n-1/2} \).
Choice D: No, the right-hand side should be \( a(t_n)u^n + b(t_n) \).
That is right if the scheme should be a Backward Euler scheme, but
the arithmetic mean on the right-hand side points to a Crank-Nicolson
scheme, though with wrong sampling of \( a \) and \( b \).
Exercise 6.15: What kind of scheme is this?
Question: We want to solve \( y'=g(x, y) \) for \( y(x) \) and have the scheme $$ y^{i+1} = y^i + \Delta t\, g(x_{i+1}, y^{i+1})\tp$$
What is this scheme called?
Choice A: The implicit midpoint scheme.
There is no midpoint \( i+1/2 \) involved here.
Choice B: The Backward Euler scheme or just the backward scheme.
Correct!
Choice C: The Forward Euler scheme, Euler's method, or just the forward scheme.
This is true if we had \( g(x_i, y^{i}) \).
Choice D: The implicit Adams scheme of order one.
Wrong!
Exercise 6.16: What kind of scheme is this?
Question: We want to solve \( y'=g(x, y) \) for \( y(x) \) and have the scheme $$ y^{i+1} = y^{i-1} + 2\Delta t\, g(x_{i}, y^{i})\tp$$
What is this scheme called?
Choice A: The implicit midpoint scheme.
It is a midpoint scheme, but it is explicit rather than implicit.
Choice B: The two-step backward scheme.
Wrong!
Choice C: The Crank-Nicolson scheme.
No, that looks quite different and is based on a centered difference
over \( [x_i,x_{i+1}] \), not \( [x_{i-1},x_{i+1}] \).
Choice D: The leapfrog scheme.
Correct!