Designing Time Filters Using a General Linear Method Approach Sigal - - PowerPoint PPT Presentation

Designing Time Filters Using a General Linear Method Approach

Sigal Gottlieb (UMassD) and Zachary J. Grant (ORNL) Joint work with: Victor P. DeCaria & William J. Layton (UPitt) International Conference on Advances in Applied Mathematics at Tel Aviv University December 18-20, 2018

Sigal Gottlieb Time Filtering Abarbanel Memorial at TAU 1 / 26

Overview

1

In memory of Saul (Shalom) Abarbanel

2

Background: Time Filters

3

Time Filters as General Linear Methods

4

Methods found using our GLM approach and optimization code

5

Conclusions

Sigal Gottlieb Time Filtering Abarbanel Memorial at TAU 2 / 26

Memories of Shalom

Uncle Shalom to me and my brothers Many Friday night dinners over 40+ years Travel and NSF panels Advice

Sigal Gottlieb Time Filtering Abarbanel Memorial at TAU 3 / 26

Motivation

Consider Layton’s filtered implicit Euler scheme: y(1) = un + ∆tF(y(1)) un+1 = y(1) − 1 3

• un−1 − 2un + y(1)

. This is an implicit Euler step (which is first order, but has nice stability properties) followed by a linear combination of some steps, which raises the order. The motivation for this approach is that we have legacy codes which have an implicit Euler scheme. Adding the ”filtering” step allows us to raise the

• rder without touching the delicate workings of the implicit legacy code.

Sigal Gottlieb Time Filtering Abarbanel Memorial at TAU 4 / 26

Generalization of Layton’s filtered Euler

We can generalize Layton’s filtered Euler method by adding a ”pre-processing” as well as a ”post-processing” step: ˜ un = dun−1 + (1 − d)un y(1) = ˜ un + ∆tF(y(1)) un+1 = 1 3 − 2d

• −un−1 + 2(1 − d)un + 2y(1)

. This extra freedom allows an entire family of second order methods that have implicit Euler as the driver.

Sigal Gottlieb Time Filtering Abarbanel Memorial at TAU 5 / 26

Generalization of Layton’s filtered Euler

With d = 0, this becomes Layton’s implicit Euler filter show above: y(1) = un + ∆tF(y(1)) un+1 = y(1) − 1 3

• un−1 − 2un + y(1)

. With d = 1, this becomes: y(1) = un−1 + ∆tF(y(1)) un+1 = 2y(1) − un−1

Sigal Gottlieb Time Filtering Abarbanel Memorial at TAU 6 / 26

Numerical results

To compare the effect of these filters, consider the problem dy dy = −10y2.

∆t Implicit Euler IE + P2

Layton Filter

IE + P2

LF Filter

5.0e-03 4.2642e-04 – 4.6646e-05 – 1.0503e-05 – 2.5e-03 2.1692e-04 0.97 1.1608e-05 2.00 2.7556e-06 1.93 1.25e-03 1.0944e-04 0.98 2.9001e-06 2.00 7.0661e-07 1.96 8.33e-04 7.3182e-05 0.99 1.2890e-06 1.99 3.1679e-07 1.97 6.25e-04 5.5000e-05 0.99 7.2557e-07 1.99 1.8138e-07 1.93

Sigal Gottlieb Time Filtering Abarbanel Memorial at TAU 7 / 26

Numerical Results

−3.4 −3.2 −3 −2.8 −2.6 −2.4 −2.2 −7 −6.5 −6 −5.5 −5 −4.5 −4 −3.5 −3

Implicit Euler Implicit Euler + Layton Filter Implicit Euler + LF Filter

Implicit Euler compared to Layton’s filter (d = 0) and the leapfrog filter (d = 1)

Sigal Gottlieb Time Filtering Abarbanel Memorial at TAU 8 / 26

General Linear Methods

Both the original method and the family of pre- and post-filtered methods can be written as a GLM. yn

i = k

• ℓ=1

diℓun−k+ℓ + ∆t

k−1

• ℓ=1

ˆ aiℓF(un−k+ℓ) + ∆t

s

• j=1

aijF(yn

j )

1 ≤ i ≤ s un+1 =

k

• ℓ=1

θℓun−k+ℓ + ∆t

k−1

• ℓ=1

ˆ bℓF(un−k+ℓ) + ∆t

s

• j=1

bjF(yn

j ),

the un−k+j are the previous steps, while the yn

j are intermediate stages

used to compute the next solution value un+1.

Sigal Gottlieb Time Filtering Abarbanel Memorial at TAU 9 / 26

Layton’s Method as a General Linear Method

Our first example is Layton’s filtered implicit Euler scheme : y(1) = un + ∆tF(y(1)) un+1 = y(1) − 1 3

• un−1 − 2un + y(1)

. We can write this as a one-stage two-step GLM: y(1) = un + ∆tF(y(1)) un+1 = −1 3un−1 + 4 3un + 2 3∆tF(y(1)), which has the Butcher GLM tableau: D =

• 1
• , Θ =
• − 1

3 4 3

• , A =
• 1
• , b =

2

3

• Sigal Gottlieb

Time Filtering Abarbanel Memorial at TAU 10 / 26

Generalized Layton’s Method as a General Linear Method

The family of second order methods with implicit Euler as the driver: ˜ un = dun−1 + (1 − d)un y(1) = ˜ un + ∆tF(y(1)) un+1 = 1 3 − 2d

• −un−1 + 2(1 − d)un + 2y(1)

can be written as a GLM y(1) = dun−1 + (1 − d)un + ∆tF(y(1)) un+1 = 2d − 1 3 − 2d un−1 + 4(1 − d) 3 − 2d un + 2 3 − 2d F

• y(1)

, with tableau D =

• d

1 − d

• , Θ =
• 2d−1

3−2d 4(1−d) 3−2d

• , A =
• 1
• , b =
• 2

3−2d

• Sigal Gottlieb

Time Filtering Abarbanel Memorial at TAU 11 / 26

Generalized Layton’s Method as a General Linear Method

The free parameter d gives us flexibility in the choice of filter. Two extreme cases to consider would be d = 0, which results in the Layton’s implicit Euler filter show above, and d = 1, which results in the scheme: y(1) = un−1 + ∆tF(y(1)) un+1 = 2y(1) − un−1 which can be written in classic GLM form y(1) = un−1 + ∆tF(y(1)) un+1 = un−1 + 2∆tF(y(1)), with the Butcher tableau: D =

• 1
• , Θ =
• 1
• , A =
• 1
• , b =
• 2
• .

Sigal Gottlieb Time Filtering Abarbanel Memorial at TAU 12 / 26

Producing Time Filters with a GLM approach

The idea behind time filtering is that certain core methods are used, and are pre- and post-processed. Consider the core method with s − 1 stages and k steps: yn

1 = un

yn

i = k

• ℓ=1

˜ diℓun−k+l + ∆t

k−1

• ℓ=1

ˆ aiℓF(un−k+ℓ) + ∆t

s

• j=1

aijF(yn

j ),

2 ≤ i ≤ s un+1 = yn

s = k

• ℓ=1

˜ θℓun−k+ℓ + ∆t

k−1

• ℓ=1

˜ ˆ bℓF(un−k+ℓ) + ∆t

s

• j=1

˜ bjF(yn

j )

Note that the final row coefficients are the same as the prior row coefficients

˜ θl = ˜ dsℓ, ˜ ˆ bℓ = ˆ asℓ, ˜ bj = asj.

Sigal Gottlieb Time Filtering Abarbanel Memorial at TAU 13 / 26

Pre-filters using a GLM approach

To pre-process the method, we modify un. To do this, we change the initial stage so that un is replaced by a linear combination of the old steps yn

1 = un −

→ yn

1 = k

• ℓ=1

d1ℓun−k+ℓ and then use yn

1 instead of un in all the middle stages of the core method:

yn

i = ˜

dikyn

1 + k−1

• ℓ=1

˜ diℓun−k+ℓ + ∆t

k−1

• ℓ=1

ˆ aiℓF(un−k+ℓ) + ∆t

s

• j=1

aijF(yn

j )

=

k

• ℓ=1

diℓun−k+l + ∆t

k−1

• ℓ=1

ˆ ailF(un−k+ℓ) + ∆t

s

• j=1

aijF(yn

j ),

where diℓ = ˜ di1d1ℓ + ˜ diℓ for ℓ = 1, k − 1, and dik = ˜ di1d1k.

Sigal Gottlieb Time Filtering Abarbanel Memorial at TAU 14 / 26

Post-filters using a GLM approach

To post-process the method, we simply modify the final line to un+1 =

k

• ℓ=1

θℓun−k+ℓ + ∆t

k−1

• ℓ=1

ˆ bℓF(un−k+ℓ) + ∆t

s

• j=1

bjF(yn

j )

where the coefficients θℓ, ˆ bℓ, bj can be chosen freely.1 We can choose these coefficients so that the final stage can be written

• nly as combinations of old steps un−k+ℓ and stages yn

j .

The coefficients θ, ˆ b and b are referred to as the post-processing coefficients.

1subject only to the order conditions Sigal Gottlieb Time Filtering Abarbanel Memorial at TAU 15 / 26

Time filtering using a GLM approach

The overall method can now be written as yn

1 = k

• ℓ=1

d1ℓun−k+ℓ yn

i = k

• ℓ=1

diℓun−k+ℓ + ∆t

k−1

• ℓ=1

ˆ ailF(un−k+ℓ) + ∆t

s

• j=1

aijF(yn

j )

2 ≤ i ≤ s un+1 =

k

• ℓ=1

θℓun−k+ℓ + ∆t

k−1

• l=1

ˆ bℓF(un−k+ℓ) + ∆t

s

• j=1

bjF(yn

j ).

where diℓ = ˜ di1d1ℓ + ˜ diℓ for ℓ = 1, k − 1, and dik = ˜ di1d1k, and d1ℓ, θℓ, ˆ bℓ, bj can be chosen for all ℓ and j.

Sigal Gottlieb Time Filtering Abarbanel Memorial at TAU 16 / 26

What’s the point?

Obviously, this is a fun game. But what does it give us? Approach the design of the pre- and post-processor as an

• ptimization problem:

Given the coefficients of the core method, we treat d1ℓ, θℓ, ˆ bℓ, bj as free parameters The equality constraints are given by the order conditions we wish the new overall GLM to satisfy Additional requirements for low storage etc. can be imposed on the structure of the coefficients We define our objective function to be minimized (error constants) or maximized (linear stability region)

We coded this optimization problem in MATLAB and quickly found methods that have desirable stability, accuracy, and efficiency properties. New GLMs based on trapezoid rule, BDF methods, fully implicit Runge–Kutta

Sigal Gottlieb Time Filtering Abarbanel Memorial at TAU 17 / 26

Filtered Trapezoid Rule

Use the implicit trapezoid rule as the core method

y (1) = un y (2) = un + 1 2∆tF(y (1)) + 1 2∆tF(y (2)) un+1 = y (2)

The filtered method takes the form

y (1) =

k

• ℓ=1

d1ℓun−k+ℓ y (2) = y (1) + 1 2∆tF(y (1)) + 1 2∆tF(y (2)) un+1 =

k

• ℓ=1

θℓun−k+j + 2 ˆ b1(y2 − y1).

Sigal Gottlieb Time Filtering Abarbanel Memorial at TAU 18 / 26

Filtered Trapezoid Rule

In the optimization code we can require the filtered scheme to satisfy third

• r fourth order accuracy conditions, to produce the P3 or P4 filters:

∆t Trapezoid Rule + P3 Filter + P4 Filter 5.00e-03 4.41e-06 – 5.02e-07 – 1.24e-07 – 2.50e-03 1.25e-06 1.81 7.47e-08 2.74 9.25e-09 3.74 1.25e-03 3.36e-07 1.89 1.02e-08 2.86 6.35e-10 3.86 8.33e-04 1.53e-07 1.93 3.14e-09 2.91 1.29e-10 3.91 6.25e-04 8.73e-08 1.95 1.34e-09 2.94 4.17e-11 3.94

Sigal Gottlieb Time Filtering Abarbanel Memorial at TAU 19 / 26

Filtered Trapezoid Rule

−3.6 −3.4 −3.2 −3 −2.8 −2.6 −2.4 −2.2 −11 −10 −9 −8 −7 −6 −5

Implicit Trapezoid Rule Trap Rule +P3 Filter Trap Rule +P4 Filter

Sigal Gottlieb Time Filtering Abarbanel Memorial at TAU 20 / 26

Filtered BDF schemes

Start with the backward Euler, which is the BDF1 scheme: y1 = un + ∆tF(y1) Now we can use a BDF2-type method to postprocess un+1 = −1 3un−1 + 4 3un + 2 3∆tF(y1) this becomes the second order scheme:

y1 = un + ∆tF(y1) un+1 = −1 3un−1 + 4 3un + 2 3∆tF(y1) = −1 3un−1 + 2 3un + 2 3y1 = y1 − 1 3(un−1 − 2un + y1) ≈ y1 − 1 3∆t2un+1

tt

Sigal Gottlieb Time Filtering Abarbanel Memorial at TAU 21 / 26

Predictor Corrector BDF-2 to BDF-3

Start with the BDF2 scheme as the core: y1 = −1 3un−1 + 4 3un + 2 3∆tF(y1) Now use a BDF3 idea to post-process2 un+1 = 2 11un−2 − 9 11un−1 + 18 11un + 6 11∆tF(y1) this becomes the third order scheme:

y1 = −1 3un−1 + 4 3un + 2 3∆tF(y1) un+1 = 2 11un−2 − 9 11un−1 + 18 11un + 6 11∆tF(y1) = 2 11un−2 − 6 11un−1 + 6 11un + 9 11y1 ≈ y1 − 2 11∆t3un+1

ttt

2We can think of this as a predictor-corrector approach Sigal Gottlieb Time Filtering Abarbanel Memorial at TAU 22 / 26

Predictor Corrector BDF-3 to BDF-4

Start with the BDF3 core method y(1) = 2 11un−2 − 9 11un−1 + 18 11un + 6 11∆tF(y1) and use the BDF4-inspired post-processor to produce

y (1) = 2 11un−2 − 9 11un−1 + 18 11un + 6 11∆tF(y1) un+1 = − 3 25un−3 + 16 25un−2 − 36 25un−1 + 48 25un + 12 25∆tF(y1) = − 3 25un−3 + 12 25un−2 − 18 25un−1 + 12 25un + 22 25y1 = y1 − 3 25(un−3 − 4un−2 + 6un−1 − 4un + y1) ≈ y1 − 3 25∆t4un+1

tttt

Sigal Gottlieb Time Filtering Abarbanel Memorial at TAU 23 / 26

Filtered fully implicit Runge–Kutta (2,2)

Consider the L-stable fully implicit Lobatto IIIC scheme: y (1) = un + 1 2∆tF(y (1)) − 1 2∆tF(y (2)) un+1 = un + 1 2∆tF(y (1)) + 1 2∆tF(y (2)) The 2 Step time-filtered scheme can be written as : ˆ u = 0.373461706729200un−1 + 0.626538293270800un y (1) = ˆ u + 1 2∆tF(y (1)) − 1 2∆tF(y (2)) y (2) = ˆ u + 1 2∆tF(y (1)) + 1 2∆tF(y (2)) un+1 = a1un−1 + a2un + a3y (1) + a4y (2)

where a1 = −0.075425887737539 a2 = 0.551112405533260 a3 = −0.596071637983322 a4 = 1.120385120187601.

This scheme is A stable but not L stable.

Sigal Gottlieb Time Filtering Abarbanel Memorial at TAU 24 / 26

Conclusions

Time filtering is an approach that can be helpful to improve the accuracy of methods within legacy codes (and maybe for other reasons as well) The time-filtering approach can be applied to multistep, multi-stage (Runge–Kutta), or multistep multi-stage methods We re-wrote the time-filtering procedure as a GLM We constructed an optimization code that can quickly find filters with a variety of desirable properties We found methods that are now being tested by Layton and DeCaria

• n interesting problems (NSE)

Sigal Gottlieb Time Filtering Abarbanel Memorial at TAU 25 / 26

The NSF-funded ICERM at Brown University

ICERM holds workshops, semester programs, small collaborative groups, and summer programs. We are actively seeking your proposals. Ask me about ICERM!

Sigal Gottlieb Time Filtering Abarbanel Memorial at TAU 26 / 26