Step 3 : generate \(Q_\Delta\) coefficients

📜 We denote by \(Q_\Delta\) the approximation of a given \(Q\) matrix related to any \(Q\)-coefficients (from collocation, RK, …) :

\[\begin{split}\begin{array} {c|c} \tau & Q \\ \hline & w^\top \end{array}\end{split}\]

There is two approaches to generate \(Q_\Delta\) approximations, from which you can choose in function of your needs and preferences.

Use a generic function

The quick easy way, simply import :

from qmat import genQDeltaCoeffs

Then considering, e.g, the following collocation method :

from qmat.qcoeff.collocation import Collocation
coll = Collocation.getInstance()  # use default parameters : 4 LEGENDRE RADAU-RIGHT nodes

we can generate the \(Q_\Delta\) matrix from a Backward Euler discretization between the nodes (as in the original SDC method) :

QDelta = genQDeltaCoeffs("BE", nodes=coll.nodes)
print(QDelta)

[[0.08858796 0.         0.         0.        ]
 [0.08858796 0.3208789  0.         0.        ]
 [0.08858796 0.3208789  0.3781926  0.        ]
 [0.08858796 0.3208789  0.3781926  0.21234054]]

… or get the LU approximation from [Weiser, 2015] :

QDelta = genQDeltaCoeffs("LU", Q=coll.Q)
print(QDelta)

[[0.11299948 0.         0.         0.        ]
 [0.234384   0.29050213 0.         0.        ]
 [0.21668178 0.48341808 0.30825766 0.        ]
 [0.22046221 0.46683684 0.44141588 0.11764706]]

… or the diagonal approximation from [van der Houwen & Sommeijer, 1991] :

QDelta = genQDeltaCoeffs("VDHS", nNodes=coll.nNodes, nodeType=coll.nodeType, quadType=coll.quadType)
print(QDelta)

[[0.32049937 0.         0.         0.        ]
 [0.         0.08915379 0.         0.        ]
 [0.         0.         0.18173956 0.        ]
 [0.         0.         0.         0.2333628 ]]

… or even some magical diagonal coefficients obtained with a black box optimizer that does not exists anymore from [Speck, 2021]:

QDelta = genQDeltaCoeffs("MIN3", nNodes=coll.nNodes, nodeType=coll.nodeType, quadType=coll.quadType)
print(QDelta)

[[0.31987868 0.         0.         0.        ]
 [0.         0.08887606 0.         0.        ]
 [0.         0.         0.18123663 0.        ]
 [0.         0.         0.         0.23273925]]

🔍 … and yes, the MIN3 coefficients are almost the same as those from VDHS, even if both where obtained independently ! But who really never ever re-invented the wheel in research 😉

Note that depending on the requested approximation, different arguments may be required for the \(Q_\Delta\)-generator called in the background. If not provided, a descriptive exception is raised, for instance :

try:
    QDelta = genQDeltaCoeffs("BE")
except Exception as e:
    print(f"{e.__class__.__name__}: {e}")

TypeError: TimeStepping.__init__() missing 1 required positional argument: 'nodes'

To allow more generic call, it is possible to give keyword arguments that will be used only for some \(Q_\Delta\) matrices :

for qdType in ["BE", "LU", "VDHS", "MIN3"]:
    QDelta = genQDeltaCoeffs(qdType,
        Q=coll.Q, nodes=coll.nodes, nNodes=coll.nNodes, nodeType=coll.nodeType, quadType=coll.quadType)

And to make it more simple, it is also possible to only give a QGenerator object as qGen keyword argument only, and qmat will extract all required arguments from it :

for qdType in ["BE", "LU", "VDHS", "MIN3"]:
    QDelta = genQDeltaCoeffs(qdType, qGen=coll)

🔔 As for \(Q\)-generators, different unique aliases exists for each \(Q_\Delta\)-generators. For instance :

QDelta = genQDeltaCoeffs("IE", nodes=coll.nodes)  # equivalent to BE

Don’t hesitate to look at the API documentation for a complete list of available generators …

Use generator objects

As for \(Q\)-generators, you can also retrieve the \(Q_\Delta\)-generators with their classes, using one of the following approaches :

• import the generator directly from its submodule

from qmat.qdelta.timestepping import BE
approx = BE(nodes=coll.nodes)

• retrieve it with one of its aliases from the QDELTA_GENERATORS dictionary

from qmat import QDELTA_GENERATORS
Generator = QDELTA_GENERATORS["IE"]
approx = Generator(nodes=coll.nodes)

In both case, you’ll instantiate an object that provides \(Q_\Delta\) matrix through a getQDelta method :

print(approx.getQDelta())

[[0.08858796 0.         0.         0.        ]
 [0.08858796 0.3208789  0.         0.        ]
 [0.08858796 0.3208789  0.3781926  0.        ]
 [0.08858796 0.3208789  0.3781926  0.21234054]]

⚠️ While the \(Q_\Delta\) generators do have a QDelta attribute, it is usually zero-initialized before calling the getQDelta method !

The reason behind the use of a method here (and not a property) is that some approximations can vary depending on a given iteration number, e.g :

approx = QDELTA_GENERATORS["MIN-SR-FLEX"](nNodes=coll.nNodes, nodeType=coll.nodeType, quadType=coll.quadType)

print("Approximation for k=1 :")
print(approx.getQDelta(k=1))
print("Approximation for k=2 :")
print(approx.getQDelta(k=2))

Approximation for k=1 :
[[0.08858796 0.         0.         0.        ]
 [0.         0.40946686 0.         0.        ]
 [0.         0.         0.78765946 0.        ]
 [0.         0.         0.         1.        ]]
Approximation for k=2 :
[[0.04429398 0.         0.         0.        ]
 [0.         0.20473343 0.         0.        ]
 [0.         0.         0.39382973 0.        ]
 [0.         0.         0.         0.5       ]]

However this corresponds to some very specific methods, in practice most of the time any \(k\) argument given to the getQDelta method will be ignored. Finally, note that all \(Q_\Delta\)-generator objects have a genCoeffs method, which per default provides the same result as the getQDelta method, but has some more functionalities (covered in a later tutorial).

🔔 As for the genQDeltaCoeffs function before, you can also use the qGen keyword argument for generic calls that will take only the name of the QDeltaGenerator as parameter.

Now, you can use those \(Q\)-coefficients and \(Q_\Delta\) approximations to build a Spectral Deferred Correction type time-stepper …