SISDC for ODEs

Let us first introduce the SISDC for ODEs for which we shall follow the original development in [19]. Assume we want to solve the ODE

u⁰(t) =f u(t), t

, u(a) =ua, t∈[a, b]. (4.7)

Before we continue, some definitions are appropriate. Lettibe a grid in the interval, a=t0 < t1 < ... < tM =b, and let u= (u0, ..., uM) be the values of u(t) at the

a

t0 t1 t_m−1 tm tM

b

nodes,ui=u(ti). Define theMth order Lagrange interpolant as

˜

u(t) =LM(u, t) = XM i=0

ui πi(t), (4.8)

where the fundamental polynomialsπi(t) are defined by πi(t) =

YM

j=0 j6=i

t−tj

ti−tj. (4.9)

The Lagrange interpolant (4.8) is a continuous function, hence we can differentiate and integrate it.

We define the integration operatorI^j+1_j as I^j+1_j u=

Z tj+1

tj

˜

u(τ) dτ, j= 0, ..., M−1. (4.10)

The integration operator integrate ˜u(t) betweentj andtj+1, i.e. the gray part in Figure 4.1. We evaluate this integral numerically using a quadrature rule on the

a

t0 t1 tj tj+1 tM

b

Figure 4.1:

form

I^j+1_j u= XM

i=0

αjiui, (4.11)

4.1.1. SISDC for ODEs 51 hence we may think ofI^j+1_j as a vector and the evaluation of the integral as an inner product ofI^j+1_j and the vectoru. Let us analyze how accurate this quadrature rule can be: We have M+ 1 nodes in between a and b, represented in the vector u.

The position of the M+ 1 nodes are set, but we have M+ 1 weights which we can choose freely. Hence we can determine theM+ 1 weights in order to integrate a polynomial of degree at most M exactly. Note that none of the nodes in the quadrature is actually inside the integration interval. The Lagrange interpolant ˜u in Equation (4.8) is a polynomial of degreeM, hence it can integrated exactly.

Notice that there areM integrations of this type, and thereby M quadrature rules, one for each interval. All the interval integrals can be determined by a matrix-vector product, theI^j+1_j being the (j+ 1)’th row in the matrix.

The grid{ti}is typically a Chebyshev grid or similar used in spectral methods, hence the Lagrange interpolant (4.8) and above integration has spectral accuracy.

Returning to the SISDC, the spectral deferred correction method is based on the Picard integral form of the ODE (4.7)

u(t) = Z t

a

f u(τ), τ

dτ+ua. (4.12)

Given an approximation of the solution v(t), let the error be δ(t), hence u(t) = v(t) +δ(t). Substitutingu(t) in the Integral form (4.12) gives

v(t) +δ(t) = Z t

a

f v(τ) +δ(τ), τ

dτ +ua. (4.13)

This is also a Picard form of an ODE, hence we could actually solve for the error δ(t), and use the result to update the approximation v. If we want p’th order accuracy, we would have to evaluate the integral top’th order accuracy, i.e. using ap’th order time integration scheme. But if we can evaluate this equation top’th order, we could just as well evaluate the original Picard form (4.12) top’th order, hence we would have p’th order accuracy already, and we would not be needing this equation.

Continue as follows: Insertingv(t) foru(t) in Equation (4.12) will not hold, but the remainder can be considered as a kind of error indicator. Hence a measure of the error is the residual equation

(t) = Z t

a

f v(τ), τ

dτ +ua−v(t). (4.14)

Subtracting Eq. (4.14) from Eq. (4.13) gives a correction equation δ(t) =

which is also on the same form as (4.12), a Picard integral form of an ODE. This equation tells how well the error measure (t) approximates the true error δ(t).

This correction equation we can solve in the same manner as we solved the original ODE, and the result can be used for updating the solution. However, we do not need ap’th order scheme for the integral in Equation (4.15).

The argument goes as follows: Leth=b−abe the interval length. Assumev(t) is orderr accurate, i.e. u(t) = v(t) +O(h^r). Then δ(t) = O(h^r). The magnitude of (t) can be found by considering the residual equation (4.14), assuming for notational simplicity thatf does not depend on time

(t) =

where we have Taylor expanded f at u(τ). We have here assumed f to be suf-ficiently smooth, such that f⁰(u) is bounded. Consider the integral in Equa-tion (4.15), call itδ(t)

Assume that we evaluate the residual equation (4.14) exactly, getting the exact (t). If we also evaluate δ(t) exactly, we would get the exact error δ(t) and the exact solutionu(t),

u(t) =v(t) +δ(t) =v(t) +δ(t) +(t). (4.18) If we evaluate δ(t) using a order stime integration scheme, then we would only get a O(h^s) correct result, i.e. ˜δ(t) = δ(1 +O(h^s)). Using this to update our solution to a new approximation, i.e. v^∗(t) =v(t) +δ(t), we would get

v^∗(t) =v(t) +δ(t)(1 +O(h^s)) +(t)

=u(t) +δ(t)O(h^s) =u(t) +O(h^r+s), (4.19)

4.1.1. SISDC for ODEs 53 hence the new solution v^∗(t) will be r orders more accurate than the original approximationv(t).

In practice, we cannot evaluate the residual equation (4.14) exactly. But if we just evaluate it to orderO(h^r+s), then that would be sufficient, not adding lower order errors tov^∗(t) in Equation (4.19). The residual equation integratesf(v(τ)), but since the approximate solutionv(t) is known, we do not need a time integration scheme as for the integral in the correction equation (4.15). Any standardO(h^r+s) quadrature rule may be used.

The articles found on spectral deferred correction [19, 49] does not include any description on how the method works, like the one presented here. They only state how to evaluate the correction equation. The presentation should give an idea of, why the method is formulated that way, and how it works.

We will now apply this setup to the i’th interval in the grid, hi = ti+1 − ti. Consider first Equation (4.12) for the i’th interval, using an explicit Euler approximation for the integral we get the approximation

vi+1=hif(vi, ti) +ui. (4.20)

Similarly, using an explicit Euler approximation for the correction equation (4.15) gives

δi+1=δi+hi(f(vi+δi, ti)−f(vi, ti)) +i+1−i. (4.21) Subtract the residual equation (4.14) at timeti from that at timeti+1, to get

i+1−i= Z ti+1

ti

f(v(τ), τ) dτ −vi+1+vi. (4.22) Combining with (4.21) gives a direct equation for the updated solution v^∗_i+1 = vi+1+δi+1

v^∗_i+1=v_i^∗+hi(f(v_i^∗, ti)−f(vi, ti)) + Z ti+1

ti

f(v(τ), τ) dτ. (4.23) To complete the correction procedure, we need to specify how to calculate the integral in (4.23). This is where we need an accurate evaluation, which we base on the quadrature rule (4.11). Letfi =f(vi, ti), and let ˜f be the Lagrange interpolant of thefi’s. We can now integrate using (4.10)

Z ti+1

which completes the scheme. A similar expression exists for the implicit Euler approximation

v^∗_i+1=v_i^∗+hi f(v^∗_i+1, ti+1)−f(vi+1, ti+1)

+Iⁱ⁺¹_i f. (4.25)

In [49] the method above is extended in a semi-implicit manner, where the right hand sidef is split into different parts treated explicitly and implicitly in a straight forward manner.

One can repeat the process with another correction step, using the new v^∗_i as vi, which will raise the order of the approximation by one each time. This can continue as long as the integral (4.24) is evaluated sufficiently accurately. Using M+ 1 points in the Lagrange interpolant provides an error of size O(h^M+1) for the integral, and an O(h^M+1) global accuracy for the solution. A complete error analysis is provided in [19].

In document A Level Set (Sider 62-66)