Coverage for pySDC/tutorial/step_1/D_collocation_accuracy_check.py: 100%
50 statements
« prev ^ index » next coverage.py v7.17.0, created at 2026-11-02 12:30 +0000
1# ---
2# jupyter:
3# jupytext:
4# formats: py:percent
5# kernelspec:
6# display_name: Python 3
7# name: python3
8# ---
10# %% [markdown]
11# # Part D: Collocation accuracy check
12#
13# As for space in [Part B](B_spatial_accuracy_check), we now measure the order of accuracy in time: solve the
14# collocation problem from [Part C](C_collocation_problem_setup) for a sequence of shrinking time steps and watch
15# the error.
16#
17# Collocation on $M$ Gauss-Radau nodes has order $2M - 1$. We solve a single step, so we see the *local* error,
18# which is one order higher: $2M = 6$ for three nodes. Beating sixth order in time with a second-order stencil in
19# space needs a fine mesh, otherwise the spatial error hides everything; hence the 16383 unknowns.
21# %%
22from collections import namedtuple
24import matplotlib.pyplot as plt
25import numpy as np
26import scipy.sparse as sp
28from pySDC.core.collocation import CollBase
29from pySDC.implementations.problem_classes.HeatEquation_ND_FD import heatNd_unforced
31problem_params = {
32 'nu': 0.1, # diffusion coefficient
33 'freq': 4, # frequency for the test value
34 'nvars': 16383, # number of DOFs in space
35 'bc': 'dirichlet-zero', # boundary conditions
36}
37prob = heatNd_unforced(**problem_params)
39# instantiate collocation class, relative to the time interval [0,1]
40coll = CollBase(num_nodes=3, tleft=0, tright=1, node_type='LEGENDRE', quad_type='RADAU-RIGHT')
42# assemble list of dt
43dt_list = [0.1 / 2**p for p in range(0, 5)]
45# %% [markdown]
46# The loop is the one from Part C, once per time step, with the results collected as in Part B.
48# %%
49# setup id for gathering the results (will sort by dt)
50ID = namedtuple('ID', 'dt')
53def run_accuracy_check(prob, coll, dt_list):
54 """
55 Routine to build and solve the linear collocation problem
57 Args:
58 prob: a problem instance
59 coll: a collocation instance
60 dt_list: list of time-step sizes
62 Return:
63 the analytic error of the solved collocation problem
64 """
66 results = {}
67 # loop over all nvars
68 for dt in dt_list:
69 # shrink collocation matrix: first line and column deals with initial value, not needed here
70 Q = coll.Qmat[1:, 1:]
72 # build system matrix M of collocation problem
73 M = sp.eye(prob.nvars[0] * coll.num_nodes) - dt * sp.kron(Q, prob.A)
75 # get initial value at t0 = 0
76 u0 = prob.u_exact(t=0)
77 # fill in u0-vector as right-hand side for the collocation problem
78 u0_coll = np.kron(np.ones(coll.num_nodes), u0)
79 # get exact solution at Tend = dt
80 uend = prob.u_exact(t=dt)
82 # solve collocation problem directly
83 u_coll = sp.linalg.spsolve(M, u0_coll)
85 # compute error
86 err = np.linalg.norm(u_coll[-prob.nvars[0] :] - uend, np.inf)
87 # get id for this dt and store error in results
88 id = ID(dt=dt)
89 results[id] = err
91 # add list of dt to results for easier access
92 results['dt_list'] = dt_list
93 return results
96def get_accuracy_order(results):
97 """
98 Routine to compute the order of accuracy in time
100 Args:
101 results: the dictionary containing the errors
103 Returns:
104 the list of orders
105 """
107 # retrieve the list of dt from results
108 assert 'dt_list' in results, 'ERROR: expecting the list of dt in the results dictionary'
109 dt_list = sorted(results['dt_list'], reverse=True)
111 order = []
112 # loop over two consecutive errors/dt pairs
113 for i in range(1, len(dt_list)):
114 # get ids
115 id = ID(dt=dt_list[i])
116 id_prev = ID(dt=dt_list[i - 1])
118 # compute order as log(prev_error/this_error)/log(this_dt/old_dt) <-- depends on the sorting of the list!
119 tmp = np.log(results[id] / results[id_prev]) / np.log(dt_list[i] / dt_list[i - 1])
120 order.append(tmp)
122 return order
125results = run_accuracy_check(prob=prob, coll=coll, dt_list=dt_list)
126order = get_accuracy_order(results)
128# We solve a single step, so this is the local error, which for a collocation method of order 2M-1 is of order 2M.
129expected_order = 2 * coll.num_nodes
130for dt, p in zip(dt_list[1:], order, strict=True):
131 print(f'dt = {dt:.5f}: computed order {p:4.3f} (expected {expected_order})')
133# %% tags=["hide-input"]
134errors = [results[ID(dt=dt)] for dt in dt_list]
135fig, ax = plt.subplots(figsize=(6, 4))
136ax.loglog(dt_list, errors, 'o', label='experiment')
137ax.loglog(
138 dt_list,
139 [errors[-1] * (dt / dt_list[-1]) ** expected_order for dt in dt_list],
140 'k--',
141 label=f'{expected_order}th order',
142)
143ax.set_xlabel(r'$\Delta t$')
144ax.set_ylabel('abs. error')
145ax.grid(alpha=0.3)
146ax.legend(frameon=False)
147fig.tight_layout()
149# %% [markdown]
150# The orders approach 6 from below. The large time steps are not in the asymptotic regime yet: for the largest,
151# $|\lambda \Delta t| = \nu (4\pi)^2 \Delta t \approx 1.6$ for the sine we start from. This test is also less
152# clean than the spatial one because the error we are after is tiny: we are computing a solution that decays
153# towards zero, very, very thoroughly.
154#
155# :::{admonition} Try it yourself
156# :class: tip
157# Switch to `quad_type='LOBATTO'`. Both ends of the interval are nodes now, so the last node is still the end
158# point and nothing else has to change. Which order do you get?
159# :::
160#
161# :::{dropdown} Answer
162# About 4.9 for the smallest step, approaching 5: Gauss-Lobatto collocation with $M$ nodes has order $2M - 2$,
163# so the local error is of order $2M - 1 = 5$. The check in the last cell therefore fails: it only looks at the
164# last, asymptotic value and demands it to be close to 6, which is exactly what tells the two node types apart.
165# :::
166#
167# ## Summary
168#
169# - Collocation on $M$ Radau nodes is a method of order $2M - 1$; one step shows the local order $2M$.
170# - Measuring temporal orders needs a spatial error well below the temporal one.
171#
172# This concludes Step 1. In [Step 2](../step_2), we stop solving the collocation problem directly and let SDC do
173# it.
174#
175# The check the tests run:
177# %%
178assert np.isclose(order[-1], expected_order, atol=0.3), f"ERROR: did not get order of accuracy as expected, got {order}"