20 - SN Partial Differential Equations and Applications Webinars

Springer Nature

5.16K subscribers

Join Markus Melenk as he investigates high order numerical methods for fractional diffusions in polygons.

Don’t want to miss another webinar? Register here: https://go.sn.pub/lXgEax​

Have an article on fractional diffusions? Submit Your Article: https://go.sn.pub/GT525V​

Learn more about publishing with SN PDEA: https://go.sn.pub/UFulhi

00:05

See the city that you see in my background an,

00:08

which is, by the way,

00:09

the birth town of Albert Einstein.

00:12

I'm a member of the editorial board of Springer Nature

00:16

PDA.

00:17

As my great pleasure to chair,

00:19

today's women are an introduce our today's speaker,

00:23

professor Marcus Mailing. Markus started his scientific career at the

00:29

University of Maryland at College Park,

00:32

where he received his master and also his PhD,

00:35

University of Maryland for people in numerical mathematics.

00:40

This rings a Bell. In fact,

00:42

both. Marcus Master and PhD thesis have been supervised by

00:47

Eva Baboushka,

00:48

who will turn 95 if a well informed on March

00:52

22nd this year.

00:54

Doctor Melink then moved to Switzerland,

00:57

where he worked as a postdoc in the group over

01:00

Chris of Schwab and received his habitation at the ETH

01:04

Zurich in the year 2000.

01:06

After acting professorships in Regensburg,

01:09

Keelan Chemnitz All Germany, Marcus moved to the Max Planck

01:13

Institute in Leipsic and then as reader to the University

01:17

of Reading since 2005 he is professor of Computational Mathematics

01:22

at the Technical University in Vienna.

01:25

Professor Malik is a world known expert in numerical analysis.

01:30

His work bridges between the analysis and numerical approximations of

01:35

white classes of piddies and gave deep insights into their

01:40

close connection,

01:41

in particular with respect to regularity of solutions and its

01:46

relation to convergence,

01:48

speed of numerical schemes is work has drawn significant attention

01:53

with about 2800 citations.

01:55

In the math scinet. Today,

01:58

Professor Mailing is going to speak on high order numerical

02:02

methods for fractional diffusions in polygons.

02:05

Marcus, the floor is yours.

02:07

Thank you very much Carsten for your kind words.

02:11

I'd like to thank also the organizers for giving me

02:15

so the opportunity to present some of my recent work

02:20

here,

02:20

and I'd like to acknowledge my coworkers in this.

02:25

This research led Bonnie from Heriot Watt and Custer Schwab

02:29

from ETA Zurich,

02:30

and in the last part of the talk I will

02:33

talk on extension to time dependent problems and that is

02:37

joint work with Alexander Reeder from the University of Vienna.

02:42

OK, so here's the outline.

02:44

So after an introduction I will say a few things

02:47

about fairly elementary,

02:49

actually about approximation on geometric measures.

02:52

Then I will present two methods.

02:54

So the first one is a galerkin method for the

02:57

extension problem,

02:58

the so-called extension problem, and the second method,

03:02

localization two is a solution is based on discretizing a

03:06

solution formula for the for the problem,

03:09

and then the last part I will talk about extensions

03:13

to time dependent problems.

03:15

OK, So what is? Fractional differential operators,

03:20

so where they became popular in recent years,

03:23

and I think the reason is that sort of applications

03:26

dries then.

03:27

So if I want to describe anomalous diffusion if wants

03:32

to solve.

03:33

Problems in computational finance and material Sciences on image processing.

03:37

It turns out that fractional differential operators are actually a

03:41

powerful tool to model these phenomena.

03:46

So what is the fractional Laplacian?

03:48

So in full space there are several essentially equivalent definitions,

03:52

so the first one is maybe the best known one.

03:55

It's in terms of four year.

03:58

Transform so you take the Fourier transform of a function.

04:02

You apply it symbol, you multiply by its symbols side

04:06

to the two S,

04:07

and then you transform back so there would be a

04:11

natural definition of the fractional Laplacian.

04:15

Full space this is another one by singular integral as

04:19

a singular integral operator says hyper singular integral operator to

04:24

be precise.

04:25

Or you could view it as a directly to Norman

04:28

operator for an elliptic problem in 1/2 space.

04:32

So let's see Caffarelli Sylvester extension.

04:35

Or you could use the done for Taylor Calculus.

04:38

We will also do that they even more.

04:42

The definitions. In bounded domains and this is of course,

04:47

but we are interested in here for the.

04:50

For computational purposes. There are these definitions extent,

04:54

but not necessarily in equivalent way.

04:57

So there are two particular cases.

05:00

The first one is the integral fractional Laplacian.

05:03

When simply restricts any of these definitions to the bounded

05:08

domain,

05:08

and then there is a second case and that is

05:12

the case of interest here.

05:14

It is the so called spectral fractional Laplacian.

05:17

That means one starts from the from eigenfunction decompositions of

05:23

the operator under.

05:24

Of the integer operator, so the Laplacian and uses that

05:29

spectral decomposition to define the spectral one.

05:33

So I give you a more details about the spectral

05:37

fractional.

05:38

Operator here on this slide,

05:40

so we're looking at L to the SU equal to

05:44

F with directly conditions and the differential operator L is

05:49

minus.

05:49

Divergent matrix a gradual where a is a uniformly symmetric

05:54

positive definite operator.

05:56

The whole thing takes place in Omega,

05:59

which is a Polygon. And we assume that a is

06:02

smooth.

06:03

In fact, we were able to place a buyer the

06:06

identity matrix in a minute.

06:07

OK, so how how is this operator defined?

06:11

As I said, we use the spectral decomposition so the

06:16

5K and Lambda K are the eigen pairs of the

06:19

differential operator L and we normalize them.

06:23

And then the application of the differential operator L to

06:27

the S is given by this formula and it works

06:30

like that.

06:31

So you take a function U.

06:32

You decompose it in the eigenfunctions 5K times scalar.

06:37

Coefficients UK and then the action of L to the

06:42

SU is simply a sum over OK.

06:44

Value Lambda K raised to the power S coefficient UK

06:48

times the eigenfunction Phi K.

06:53

OK, so with this definition this problem on top here

06:57

is easily solved.

06:59

If you decompose F in the eigen functions.

07:04

Then the solution you simply the sum sum of RKU

07:07

K5K where the UK is given by Lambda K to

07:10

the minus S *

07:11

F. OK, so we have a solution formula already and

07:15

you may want to.

07:16

So what's the big deal about numerical methods?

07:19

Well in practice of course the eigenfunctions are not known

07:23

to us,

07:23

so this definition is not really numerically tractable.

07:29

This problem comes with.

07:32

Energy with spaces. HS or H Sigma and they are

07:36

simply given by sort of decay conditions on the coefficients.

07:40

When you expand the function in terms of the eigen

07:43

functions.

07:44

So this is here. The formal definition for relevant to

07:48

us.

07:48

We will only be 3 cases.

07:50

Sigma equal to 0 corresponds to the L2 space Sigma

07:54

equal to 1 corresponds to the space H10.

07:57

So the 0 boundary condition is encoded here.

08:00

In this I value problem in the case S equal

08:03

to Sigma equal to S corresponds to the energy space

08:07

of this elliptic problem.

08:09

So in fact the L2,

08:10

DSHS, glyptic.

08:13

And I said already, we will focus you on the

08:16

case.

08:17

Matrix A is Laplace is identity,

08:19

so we're looking at the fractional Laplacian.

08:23

OK, So what is the goal of today's talk?

08:25

So I want to design numerical methods based on differential

08:29

operators on local operators,

08:30

and I want to discuss high order methods,

08:33

and in particular I'm shooting for exponential convergence for Christmas

08:38

analytic data.

08:39

I said already, I will do that for two different

08:43

types of methods.

08:44

And I want to highlight here one key things of

08:47

the one unifying theme for these two different discretizations.

08:52

And that is the measures that I will use a

08:54

geometrically refined towards vertices and edges.

08:58

So in this 2D case here with a smooth boundary.

09:01

So you see you have measures that are fairly dramatically

09:06

refined towards the.

09:07

Towards the boundary. So they are long thin elements.

09:12

And here in the case of polygons you have the

09:14

same thing.

09:15

Here you have these small thin elements,

09:17

but near the boundary you actually have further geometric refinement

09:22

towards the corners in order to resolve corner singularities.

09:27

OK, so let me start with something very simple about

09:32

geometric meshes.

09:33

So that yeah, so the first question that I want

09:36

to discuss is why should I care about high order

09:39

methods.

09:39

So this is a very simple questions.

09:41

Let's have a smooth function here.

09:43

This black curve and I want to approximate it with

09:46

let's say a piecewise linear function.

09:48

So I do have piecewise linear interpolant,

09:50

or I want to approximate with a higher order polynomial.

09:54

So in this particular case,

09:55

are polynomial of degree 4 Enter chosen Chebyshev interpolation?

09:59

So already qualitatively, you see that you polynomial of degree

10:03

4 so much better.

10:03

Proximation, this piecewise linear one.

10:07

OK, so complete error analysis for this kind of interpolation

10:11

is of course available,

10:12

so for the piecewise linear approximation the error scales like

10:17

H squared,

10:17

where H is the mesh spacing of these interpolation nodes.

10:22

So that leads to algebraic convergence.

10:25

So the error in some norm here is end to

10:28

the minus two,

10:29

where N is the number of knots that we had

10:32

here.

10:33

So he had 12345678.

10:35

In the Chebyshev interpolation. So if you have smooth functions

10:39

then you get eat,

10:41

for example two to the minus P some something that

10:44

is even smaller times apiece plus first derivative of F.

10:47

So as P group grows.

10:49

So this goes to zero very quickly,

10:51

so you can actually expect exponential convergence in the polynomial

10:57

degree.

10:57

OK, so that tells you that.

11:00

I order approximation could be very accurate and then need

11:04

quantitatively you see that here.

11:06

So for example, if you want to have with a

11:09

piecewise linear is an error of 10 to the minus

11:12

four,

11:12

you need something like 200 knots.

11:15

So while you tweet with 10 to the minus four

11:18

here with Chebyshev interpolation,

11:20

you get by with less than 9 pretty normal degree

11:23

9,

11:23

so 10 points.

11:25

OK, so you clearly see the advantage of higher order

11:29

methods.

11:30

OK, so is there a catch?

11:31

Well, you see that here in the error estimate.

11:34

So you need of course for polynomial degree P you

11:37

need.

11:38

Of course higher order derivatives.

11:40

So what happens when that is not the case?

11:44

So for example, if we try to approximate the function

11:46

X to the Alpha,

11:47

where Alpha is between zero and one and the function

11:50

is of course not smooth at X equal to 0

11:52

because derivatives will blow up and then a Chebyshev interpolation,

11:56

it does not work so great.

11:58

So actually you see an algebraic convergence of 1 on

12:02

P.

12:02

So the reason is that the function is not smooth.

12:06

So what should we do?

12:09

Well, the answer is on the next slide.

12:11

If I still want to have good approximation by polynomials,

12:14

I should use piecewise polynomial approximation on appropriate measures,

12:19

and these are the so called geometric measures that will

12:22

appear frequently in this talk and This is why I

12:25

want to sort of spend some time on that.

12:28

So here we have the 1D case,

12:29

so we approximate the function X to the Alpha save,

12:32

and we know that the singularity is at X equal

12:36

to 0.

12:36

So let's so we need to refine the mesh there

12:40

in order to resolve that singularity.

12:42

So for example, on the next mesh we refine and

12:45

then by splitting into two and then we find again

12:49

by splitting again into two,

12:51

but only the left hand element,

12:53

and then we refine again,

12:55

but only the again the leftmost element and so forth.

12:59

So a geometric mesh. More generally,

13:02

speaking with layers and grading Factor,

13:05

Sigma is given by the points just left end .0

13:08

and then Sigma to the L Sigma to the L

13:11

Sigma to the R

13:12

1 and so forth up to Sigma to 0,

13:14

which is one OK. So these are the measures,

13:17

so they are strongly refined towards an appoint.

13:21

So how should we approximated by polynomials?

13:24

Well, let's say on the courses mesh we use polynomials

13:28

of degree one on the next match we will use

13:31

polynomial.

13:32

We two on each of these elements and then three

13:35

on each of these elements,

13:37

and so forth. So you see,

13:39

we simultaneously increase the approximation order,

13:43

the polynomial degree and refine the mesh near one point

13:48

and for this.

13:50

Case then, so for extra,

13:51

the Alpha is piecewise. Chebyshev interpolation actually converges exponentially in

13:57

the polynomial degree P and the number of nodes.

14:00

So the number of.

14:02

And the ones that we have for this interpolation problem

14:06

is P times Peter plus one.

14:07

So the number of elements times the polynomial degree know

14:11

dimensional polynomial P.

14:12

And here's the plot.

14:15

So we have the arrow versus the number of knots

14:18

in here.

14:18

This red curve tells you what shows you that you

14:21

really get this exponential convergence in the polynomial EP.

14:28

OK, so that extends to 2D,

14:29

so so This is why the geometric measure idea sort

14:33

of generalizes to do and higher dimensions.

14:35

So for example, in elliptic problems in 2D is a

14:39

typical solutions have corner singularity,

14:41

so the solution behaves like R to the elf,

14:44

attempt something smooth and R is the distance from a

14:48

point.

14:48

Let's say the vertex of the domain.

14:50

If you want to do the approximation there you do

14:54

a geometric refinement towards that trouble some point.

14:58

You see you so it's a self similar refinement towards

15:01

that point.

15:02

So the elements near the singularity,

15:04

again very very small.

15:07

And by selecting our proportional to the polynomial degree,

15:11

you get exponential convergence.

15:15

Edge Singularity says this is a sort of archetype of.

15:18

This is something similar, let's say in 2D you have

15:21

a function of X&Y and it's smooth in X,

15:24

but you have a nasty why dependence,

15:26

so it's white with the Alpha.

15:28

In that case, again, you do a geometric mesh refinement

15:32

towards the edge Y equal to 0 by so you

15:35

do a geometric mesh in One Direction and.

15:38

Node mesh refinement in the other direction,

15:40

so you get this on isotropic elements.

15:42

Again, you get exponential convergence.

15:46

OK, so these are was a long prelude on geometric

15:49

measures.

15:50

Let me start now with the the first numerical methods.

15:53

The Caffarel based on the Caffarelli Silvestre extension.

15:57

This is So what is the Caffarelli Silvestre extension?

16:01

So we enter the function U on this domain Omega,

16:04

let's embed this Omega into a bigger into a by

16:08

increasing the dimension by one so we extend it into

16:11

our three.

16:12

And Omega is now simply the bottom of an infinite

16:16

cylinder.

16:16

See which is Omega Cross zero Infinity?

16:20

OK, and the extra variable that we introduced now is

16:23

called Y and the rest of the talk.

16:25

OK, so we have this infinite cylinder and we will

16:28

seek a function capital U on this infinite cylinder so

16:32

there will be a relation between the function little U.

16:35

This function that we are interested in their function capital

16:39

U and that is that you is that race.

16:42

On the bottom of the cylinder of this function capital

16:45

you.

16:46

So what happens to capital you capital?

16:48

You solve an equation. I hold you equal to zero

16:52

in the cylinder,

16:53

and Allison elliptic operator. So it's minus the versions of

16:57

Y to the 1

16:58

2 S grad U. So it's an elliptic use of

17:01

an elliptic equation.

17:03

OK, so we need to satisfy here these boundary conditions.

17:07

This is accounted for by saying that requiring that U

17:11

is 0 on the mantle of the cylinder and we

17:13

somehow have to bring in the data F and that

17:16

is headed to normal derivative of U on the bottom

17:20

of the cylinder is up to a scaling factor.

17:22

DS the function F.

17:25

So, and this is the Caffarelli Sylvester extension.

17:29

So what do we know about the function U?

17:32

I've written that down here once more.

17:35

Think for a moment about the case as equal to

17:38

1/2,

17:38

then 1 2 S 0.

17:39

So you see that actually capital U solves is harmonic.

17:43

It solves Laplace equation.

17:46

And more generally, actually you expect that you is analytic

17:49

in the cylinder simply because the coefficient Y to the

17:52

1

17:53

2 S is analytic.

17:55

It's an elliptic problem in a corner domain because we

17:59

have here these.

18:00

These edges also here on the in the green edges,

18:03

and so you expect algebraic singularities at the edges and

18:07

vertices of the.

18:08

Cylinder.

18:10

And you also, since the coefficient Y to the minus

18:14

2 S,

18:14

degenerates. You expect an algebraic singularity at Y equal to

18:20

0.

18:20

OK, so a lot of these things are simply expectations,

18:23

So what can one show?

18:26

So, for example, the analyticity in the function in the

18:30

variable Y can be shown explicitly,

18:32

so here's the reference. So this is the theorem that

18:37

comes with it.

18:38

So.

18:39

Why is it? Why does it seem Taylor said it's

18:43

analytic?

18:43

So you see the function capital U.

18:46

If you take the elf derivative then it behaves like

18:49

L factorial and that happens for all natural numbers L.

18:53

So that is a hallmark of analyticity.

18:58

Another thing that we know is that actually decays exponentially

19:02

in Y.

19:02

That is what this factor tells us.

19:06

And there's an algebraic singularity at Y equal to 0,

19:09

and you see that here in this highlighted term,

19:12

because for each differentiation in.

19:15

For the Earth derivative you have why to the minus

19:18

L,

19:19

so each time you differentiate you pull out a factor

19:24

one on one on Y.

19:25

OK, One knows a little bit more about the solution

19:28

behavior.

19:29

So for example, the trace on Omega,

19:31

so the function little yo behaves like the distance to

19:34

the boundary raised to a power that is not an

19:37

integer.

19:37

So here's a picture of what that looks like.

19:40

So you see the typical algebraic singular edge singularity of

19:44

the solution,

19:44

so that is what's known.

19:48

How should we discretize?

19:51

How should we solve for you?

19:53

Well, we follow here this reference and do a classical

19:57

finite element method on the cylinder.

19:59

So the bilinear form that comes with it is is

20:02

given here.

20:03

The simply obtained by multiplying by test function and integrating

20:08

and integrating by parts.

20:10

There's an energy space that is simply given here by

20:14

those functions that have finite energy.

20:17

So the bilinear form AC is finite and we need

20:20

to introduce we need to incorporate the boundary conditions.

20:24

So the vanishing on the mantle.

20:27

So this is this function,

20:28

H1. Circle, OK, the variational formulation is then completely classical,

20:33

so it's given here by this expression.

20:37

And to get lurking methods is also fairly simple.

20:41

You simply replace the infinite dimensional space H1 by a

20:47

closed subspace VN.

20:49

And this is a nice elliptic setting,

20:51

so actually one has twice the optimality in this energy

20:55

North.

20:55

So the finite element solution UN exists is unique and

20:59

it's error measured in the energy norm is the minimum

21:02

overall approx possible approximations from the finite element space VN

21:07

in the energy nor. And remember that ultimately we're interested

21:11

in that race,

21:12

so we need to relate the energy norm.

21:15

Two norm on the Omega,

21:16

but we have a trace estimate,

21:18

so the HS. Norm. Is actually bounded by the energy

21:22

norm,

21:22

so convergence in the energy non tells us something about

21:25

conversions in the Energy North.

21:27

OK, so sorry, this quasar optimality result tells us that

21:31

we need to worry about.

21:32

We only need to worry about this.

21:34

The selection of the space VN.

21:37

How should we choose it?

21:39

And we will choose it as a tensor product so

21:43

it is a tensor product of a space of piecewise

21:46

polynomials of degree P on a mesh T which lives

21:49

here in Omega.

21:52

10s arise with a space of piecewise polynomials of degree

21:57

P on a mesh GL that lives in the on

22:01

the 1D.

22:01

Interval 0 code 8Y.

22:06

OK, so let's first focus on the mesh GL to

22:11

you.

22:11

We know that the solution decays exponentially in Y,

22:14

so that's good news because that means we can simply

22:18

cut off at some point.

22:19

Cody, why, and probably not make a big error.

22:23

The next thing we know is that the solution has

22:26

a singularity at Y,

22:28

and it's an algebraic singularity.

22:30

So from the prelude we know that a geometric mesh

22:33

refinements was Y equal to 0.

22:35

Is a good choice, so this is the mesh.

22:38

The mesh GL Geo that we will use for the

22:41

rest of the talk.

22:42

So the interval zero curly Y is discretized with a

22:45

mesh with L layers,

22:47

so L + 1 elements.

22:48

So the smallest element is Sigma to the other times

22:52

Y and so forth.

22:53

Sigma to the Y and then why?

22:55

OK, so this is the space that we will use.

22:59

So the question is, how should we choose to mesh

23:03

in Omega so that E.

23:07

So the morally speaking we know already what's going on

23:10

because we know that we should capture these algebraic singular

23:14

to use these edge singularity is near the boundary of

23:17

the Omega.

23:19

So now we expect geometric measures to work well.

23:24

Technical trouble is that we don't have a full regularity

23:28

theory for the solution of the extended problem,

23:31

so This is why we proved an exponential convergence result

23:35

in a different way.

23:37

So let's so this is the finite element space that

23:41

we want to use.

23:42

To motivate the the correct choice of space.

23:46

SPT so of the triangulation team.

23:49

Let's look at this semi discretization case first.

23:53

So let's replace the SPT just by the infinite dimensional

23:59

subspace H10.

24:01

OK, so and then this is again a closed subspace.

24:05

So EU semi will be the galerkin approximation to the.

24:10

Function you.

24:13

OK, then one can show this is essentially the analyticity

24:17

of the function.

24:18

You and the exponential decay that actually the semi discrete

24:22

approximation error is bounded by these three terms.

24:26

So the first one is exponentially in the polynomial degree

24:29

P,

24:29

so that is disappear.

24:32

The second term is exponential.

24:35

In the number L of layers of the geometric mesh

24:38

here.

24:39

So essentially each of the minus B prime L is

24:42

related to the size of the smallest element here by

24:46

Y equal to 0.

24:47

And the third term here each of the minus B

24:50

double prime curly.

24:51

Why is related to the cut off that we make

24:54

an error by simply cutting off?

24:55

But that error is small as you see.

25:00

Balancing these three terms tells us already what how we

25:03

should choose Curly.

25:04

Why Alan P? There should be proportional to each other.

25:09

So that will fix for the rest of the talk.

25:12

The relation between Curly Weick LNP.

25:17

OK, so yes again, the geometric mesh that people use.

25:21

Semi discreet space.

25:25

Um, I want to get my hands on the this

25:29

semi discreet solution.

25:32

So the.

25:34

On this use me. I want to get my hands

25:37

on it.

25:38

For that it is clever to actually find a good

25:42

basis of this space SP GL.

25:45

OK, so for that we use we solve an eigenvalue

25:48

problem you'll see on the next slide why this is

25:51

a clever choice.

25:52

So we saw the Viis and immune eigenvalues mu eyes.

25:56

So this eigenvalue problem. So it's a 1D eigenvalue problem

26:00

on the interval zero Y.

26:01

OK, here you recognize the.

26:03

Of course the weight that comes with the extension problem.

26:08

NB We normalize things here.

26:11

OK, So what I'm plotting here in the left curve

26:15

is the the eigenvalues mu I of this eigenvalue problem.

26:19

For the case of a geometric mesh with 12 layers

26:22

and polynomial degree 12.

26:24

So you have roughly 100 eigenmodes.

26:28

And then plot and plotting them.

26:30

Here the value over the number.

26:32

So and you see that most of them,

26:34

so they ranges from 10 to the minus 22 essentially

26:38

one.

26:38

So you see that most of the eigenvalues are very

26:42

very small.

26:43

OK, one can quantify that.

26:45

Here's an quantitative estimate. So the mu eyes are between

26:49

L squared and Sigma 2,

26:51

the two L. So this is essentially small T square

26:55

off the smallest element in the mesh GL.

27:00

OK, so we have these eigenvalue these eigenfunctions VI.

27:04

So now we make the answers that is semi discrete.

27:08

Solution is functionally unknown functions UI which live in H10

27:13

times these.

27:14

These piecewise polynomials VI of Y.

27:17

OK, no one can compute.

27:19

It's a fairly elementary calculation that the UI satisfy this

27:24

differential equation or this variational formulation.

27:28

So mu I times the integral over graduate Grad V

27:32

plus the integral of a UI times V is.

27:35

V times.

27:37

F * V The integral of F *

27:39

V. OK, so this is a variational problem.

27:43

You see it's an elliptic problems for the for the

27:46

US.

27:46

Essentially it is minus UI Laplace.

27:48

You plus you equals F.

27:51

Now you notice that the mule has most of them

27:54

are small,

27:54

so this is actually a singularly perturbed problem.

27:57

OK, this convenient to abbreviate this bilinear form here.

28:02

And the next thing we of course not interested in

28:05

this semi discrete approximation.

28:07

Ultimately we are interested in the fully discrete galerkin approximation,

28:11

and then it's a calculation that a fully discrete approximation

28:15

so the approximation from.

28:17

From the space SP of T tenderized with SP of

28:21

G as this form,

28:23

so the UN is the sum over Piu I times

28:27

VI.

28:27

Well, the P eyes are simply the galerkin approximations of

28:32

the functions UI in the space SP of T.

28:35

So they satisfy this best approximation property.

28:41

So that is in a way we have a solution

28:44

formula for the fully discrete problem,

28:46

but in terms of the semi discrete approximation,

28:49

because so the semi discrete approximation,

28:52

then we apply a galerkin approximation to that and then

28:56

there we have a fully discrete approximation.

28:59

So it's numerically also interesting that we have a nice

29:03

decoupling technique but for the for all purposes here.

29:07

We see that actually the.

29:11

This formula tells us how we should choose the mesh.

29:14

Namely, we should choose it in such a way that

29:18

the.

29:19

Pi of US can be can approximate well the solutions

29:24

you are.

29:25

And remember that the sulat you eyes are solutions of

29:29

elliptic prompts of singular perturbation type.

29:32

So ultimately we have to design the mesh such that

29:36

we can solve.

29:37

Then we can approximate a bunch of singularly perturbed problems

29:41

well,

29:42

and the range of this singular perturbation problem of parameter

29:46

is fairly extreme,

29:47

so that it can be 10 to the minus 15

29:51

to size 1.

29:52

OK, so we've reduced things to understanding solution of singular

29:56

perturbation problems.

29:58

Um, I've done that in my life,

30:00

so let's look at first the problem of a single

30:03

scale.

30:03

So we look at minus epsilon squared plus you,

30:06

plus you equal to F.

30:08

So here's the 1D case.

30:10

So the typical solution behavior for smooth right hand side

30:13

is like that.

30:14

So you have a smooth function solution and then in

30:18

a small O of epsilon neighborhood you have a sharp

30:22

layer.

30:22

OK, so how do you resolve that numerically?

30:25

Well, you probably need a small element of size epsilon

30:29

near the boundary layer,

30:30

so at the left end point and then the right

30:33

endpoint in order to resolve that layer.

30:36

If you approximate with piecewise polynomials of degree PA,

30:40

good choice is actually to have this width P times

30:45

epsilon.

30:46

If you do that 2D,

30:47

then the correct mesh looks like that.

30:49

So you have, let's say a course measure elements and

30:52

then here you have a long thin element where the

30:55

width is again P times epsilon.

30:57

The reason why you need long thin elements is that

31:00

these boundary layers essentially all the trouble.

31:03

Some effects are 1D, so they very sharply in the

31:06

direction normal to the boundary,

31:08

but they're very smoothly in the direction tangential to the

31:12

boundary.

31:12

So This is why long thin elements are good choice.

31:16

Four corner domains things a little bit more complicated,

31:20

so you need again. For example,

31:22

here. These long thin elements and near the vertices you

31:26

need to resolve the corner singularities and this is what

31:29

happens here on and off.

31:31

P times epsilon neighborhood there you need a geometric mesh

31:35

refinement as I've shown you in the prelude,

31:38

so again itself similar refinement towards it corner.

31:41

So what you see here is a blow up off

31:44

of this region.

31:45

If you do that, you get exponential convergence in the

31:50

polynomial degree on these kind of measures uniformly in the

31:56

in the.

31:56

Singular perturbation parameter epsilon.

32:00

OK, so this is a single case upscale case.

32:02

What happens if epsilon is in a range of?

32:06

In a range between epsilon mean and once a well.

32:09

Of course you unit sense sense.

32:11

You need to use the union of all these measures

32:14

and the way to conveniently realize that is you do

32:17

geometric mesh refinement towards the boundary.

32:20

So for example in the ones in case you do

32:23

a geometric mesh with the smallest element being Sigma to

32:26

the L,

32:27

where the Alice chosen such that Sigma to the L

32:30

resolves the smallest length scale,

32:32

so is proportional to epsilon min,

32:34

and then it's at least plausible that sort of.

32:37

Each of. Each. For each epsilon that is between epsilon

32:42

minion epsilon one,

32:44

you have a you have a mesh that looks like

32:46

that sort of.

32:47

This is finer mesh than the one that you would

32:51

need.

32:51

Same thing in 2D. You would use a geometric mesh

32:55

refinement towards the boundary with the smallest element be having

32:59

with Sigma to the L.

33:00

And I'll tell you something about what happens on the

33:04

in corner domains on the next slide.

33:06

OK, then one can show that the.

33:09

HP Finite element method converges again exponentially in the polynomial

33:14

degree P in the energy norm uniformly in.

33:17

In the.

33:19

In the singular perturbation parameter epsilon.

33:24

OK, so how should we do?

33:26

Use a chosen measures for corner domains.

33:31

The simplest case is of course,

33:32

a square.

33:35

Then Mom would use a tensor product mesh.

33:38

In the general case.

33:40

So for example here for this L shaped domain,

33:43

I would propose to use to construct these measures we

33:46

have mesh patches,

33:47

so you start with a.

33:48

Course triangulation, like here and then.

33:51

Each of these elements is subdivided further.

33:55

So each of these patches.

33:57

So for example, let's start with this one here.

33:59

So you want a geometric refinement towards here,

34:02

this boundary edge. So you do that here by geometric

34:06

mesh refinement towards the edge.

34:08

So what happens here at the at the corner?

34:11

So you have two edges,

34:12

so you need to do geometric mesh refinement towards both

34:16

edges,

34:16

and so you use in a way it ends up

34:19

product of geometric measures.

34:22

This corner here. So there is a singularity here.

34:25

So you need a geometric mesh refinement towards the singular

34:28

point,

34:29

and again we've seen that self similar refinement towards the

34:32

singularity already.

34:34

And then there's the last kind of patches you need.

34:36

A transition between the boundary layer Patch and the corner

34:40

late Connor Patch,

34:41

and this is what I call the mixed Patch,

34:43

so it is simply is meant as a transition between

34:47

this one and this one.

34:48

And the other elements sort of.

34:50

Their inside are not refined anymore.

34:55

OK, so let's go back to our.

35:00

Caffarelli Sylvester extension problem. So we're going back now to

35:05

the fractional diffusion problem.

35:07

So if we use. Let's say in one the the

35:11

mesh T that is geometrically refined towards the.

35:15

And points on 2D. Let's say a geometric mesh that

35:19

is geometrically refined towards heat boundary then.

35:23

And this is the finite element space that are used

35:27

now.

35:27

So piecewise polynomials on these geometric meshes 10s arise with

35:33

the.

35:33

SP of GL.

35:36

So this is the discretization dimension in the Y variable.

35:39

Then one can actually show this is here from this

35:43

reference.

35:44

If the geometry is smooth 1D or analytic,

35:46

then one can actually show exponential convergence in the polynomial

35:50

degree.

35:50

If the polynomial degree is chosen proportional to the number

35:55

of layers and proportional to the cut off parameter Y.

35:59

OK, so this is exponential convergence in P.

36:01

If you check what the problem sizes,

36:04

so the VN then you find out it is P

36:06

to the spatial dimension plus 3.

36:09

Turn 2D for example. This is Peter the five.

36:14

In corner domains so the polygons.

36:16

So now strictly too deep you need to do the

36:19

geometric refinement towards edges and vertices.

36:22

So for example here if this is the boundary and

36:25

then this is the message that we will use.

36:28

So you do a geometric mesh refinement towards the edges

36:32

and here near the.

36:34

Vertex of the domain. You do a geometric mesh refinement

36:39

towards the vertex.

36:40

Then you get again exponential convergence in the polynomial degree

36:44

P and the problem size now is Peach is A6

36:47

rather than P to the five.

36:51

OK, so we calculated this also numerically.

36:54

So what I'm plotting here for the,

36:56

let's say, the square domain,

36:58

so you see here the mesh that we use.

37:00

The error in the energy norm versus the polynomial creepy

37:03

for different values of S.

37:05

And you clearly see exponential convergence in the polynomial degree.

37:10

I should say that the numerical results are obtained with

37:15

NGS software developed by your insurance group in Vienna.

37:20

We did also calculations from L shaped domain and here

37:23

for a slit domain.

37:24

And qualitatively you see the same thing.

37:26

Exponential convergence in the energy norm.

37:31

OK, I want to tell you a few things about

37:34

the second method,

37:36

the so-called sync quadrature method.

37:39

So this is now a new method.

37:41

Starting point is your idea of bonito and pashak to

37:45

start from a representation formula for the solution.

37:49

So if the function little you that we seek has

37:53

actually a form out to the minus F *

37:56

S * F. So we have a solution operator SL

37:59

to the minus S has this form,

38:01

so it is an integral over Lambda.

38:05

Lambda identity plus differential operator which is here the fractional

38:09

upper foot which is the Laplace operator to the minus

38:12

one.

38:13

OK, so we can rewrite this integral.

38:15

So how should we discretize that?

38:17

Well, we start with discretizing the integral.

38:21

So the simplest thing to do is.

38:24

To do a discretization with the trapezoidal rule so this

38:29

infinite integral from minus Infinity to Infinity.

38:32

This result on a uniform mesh with M + 1

38:35

points.

38:36

The quadrature points YJR J *

38:38

M to the minus one half and J runs from

38:41

minus M2M.

38:42

So in other words, is a uniform mesh with spacing

38:46

1 /

38:46

sqrt M OK and then the approximation to this this

38:50

solution operator is Q to the minus S which is

38:53

simply a finite sum.

38:55

Of solution operators for elliptic problems.

38:58

So it's identity plus epsilon squared,

39:00

LA plus or minus plus.

39:02

Races apart minus and the epsilon is given by this

39:06

formula here.

39:07

OK, of course when discretized,

39:09

this semi discretization incurs an error,

39:12

but the error is actually exponentially small in the number

39:15

of quadrature points,

39:17

so it's each of the minus square root of M.

39:19

So in other words.

39:21

That is under control.

39:24

So here is again the solution formula.

39:26

Of course, for a fully discrete method one needs to

39:30

discretize these solution operators.

39:32

So in other words, we need to solve.

39:36

Elliptic problems, and we do that by a high order

39:40

finite element methods sobia.

39:42

Piecewise polynomial approximation.

39:46

By putting on something creepy on a mesh tee.

39:49

So how should we choose the mesh T?

39:52

Well, so we're discretizing an elliptic problem identity plus epsilon

39:58

squared times Laplace.

39:59

Now these epsilons here. Here's the formula.

40:02

Most of them are small,

40:04

or at least the troublesome ones are small,

40:06

so we need to be able to solve singularly pitter

40:09

problems.

40:10

Now we have told you already how to do that.

40:14

We can do that with.

40:16

Geometrically, find measures where the measures are.

40:21

Optically refined towards the boundary.

40:24

And this is what what one can do.

40:27

And then one obtains a scheme.

40:30

That converges exponentially in P,

40:33

so in the appropriate norm,

40:35

if we do the two things,

40:37

the first thing is that for the discretization of the

40:42

geometry T.

40:45

So often Omega we use the geometric mesh with L

40:48

layers of refinement where the number of layers of refinement

40:52

is proportional to P and proportional to square root of

40:56

M, because E to the minus sqrt M /

40:59

2 is the smallest length scale that is induced by

41:04

the by the quadrature.

41:06

And.

41:08

Of course we need not to remember that.

41:13

If we have an error here,

41:14

each of the minus square root of M,

41:16

so we need to balance errors.

41:18

So square root M should be proportional to P.

41:21

If you do the bean counting,

41:23

then actually the problem size for smooth domains is also

41:27

P2 to five and for polygons is also P2 to

41:30

six,

41:30

so it's the same problem size as previously.

41:35

We did some numerics with for this.

41:38

So here's some for the certain parameter choices.

41:42

We computed the error. And we did not compute the

41:46

HSR robot.

41:47

We computed something similar to the.

41:51

Energy norm error that we use for the first method

41:54

so it's the integral of F times exact solution minus

41:58

numerical approximation.

41:59

So we we plotted the error here versus the polynomial

42:03

degree and you observe again exponential convergence here for this

42:08

square domain.

42:09

Should we wonder is this better than the first method?

42:13

So here I'm plotting for the same polynomial degrees?

42:16

Is the convergence of the first method selection method and

42:20

it seems like the alert the Balakrishnan formula here.

42:23

So the second method is slightly better,

42:26

but one should be careful here with comparing parameters,

42:29

so if when.

42:32

For this particular choice of parameters,

42:34

if one plus the error versus the number of linear

42:37

systems that one actually has to solve.

42:39

So for example you for this Balakrishnan formula,

42:42

we had a bunch for each quarter point we had

42:44

to solve an elliptic problem,

42:46

so these are fully decoupled problems and similar decoupling trick

42:50

can be done here for the galerkin method.

42:52

So a fair measure is really the number of linear

42:55

systems that one has to solve for and then actually

42:58

at least for this parameter choice,

43:00

the Gallery method comes out slightly better.

43:04

We did the same calculation here for the.

43:09

L shaped domain name for the slip domain.

43:12

We compare again here. The Balakrishnan formula with the galerkin

43:17

method,

43:18

and also for the arrow versus number of linear systems

43:22

that need to be solved.

43:24

So the same qualitatively the same behavior As for the

43:28

square domain hands.

43:29

OK, the last few minutes I want to switch topic

43:34

and talk about time dependent fractional problems because these.

43:39

Techniques that are developed here for the for the elliptic

43:43

case extend to the to the time dependent problem.

43:47

So the first cases here actually the fractional heat equation.

43:51

So you have U dot plus the fraction Laplacian applied

43:55

to you equal to F.

43:57

We're it's a spectral fractional Laplacian here,

44:00

and so we apply Billingsley conditions and of course we

44:04

need initial condition,

44:06

so at T equals to zero we impose conditions 0.

44:10

So if we know how to discretize the fractional Laplacian.

44:14

So how do we deal with the the time dependence?

44:18

Well, essentially one can use any kind of.

44:24

Spatial sorry time discretization technique,

44:28

so we in particular use actually a DG method in

44:31

time,

44:32

but you could also do a whole quarter method if

44:35

you wished.

44:36

OK, So what is the new feature here?

44:40

The new feature here is that at time T equal

44:44

to 01 actually has a singularity that is typical of

44:48

these parabolic problems.

44:50

So we need to do again,

44:52

Jim.

44:53

Refinement towards cheaper to zero and what we propose here.

44:57

Following here, certain Schwab is actually to do a geometric

45:01

mesh refinement towards time T equal to 0 because the

45:05

solution has an algebraic singularity,

45:07

there a fractional time derivative with camera being not an

45:11

integer like in the Balakrishnan Formula One actually has a

45:15

solution formula and this is the solution.

45:18

Formula is maybe a little bit daunting,

45:20

but go through the step by step.

45:23

So the first thing is that.

45:24

Just like in the publication of Formula One has a

45:27

contour integral,

45:28

so there's a need to build.

45:30

It essentially goes from minus Infinity plus Infinity,

45:33

and we will use again sing quadrature for that.

45:36

So just like in the elliptic case,

45:38

that integral is discretized with us in quadrature.

45:41

Then we have here these object E gamma,

45:43

wanandi, gamma, gamma, those armytek Leffler function so that entire

45:47

functions explicitly known.

45:48

So in a way they are under control.

45:51

The next object here is the Z plus.

45:54

Laplace raises up on this one,

45:56

so these are of course elliptic solution solution operators for

46:01

elliptic.

46:02

Problems and the choice of the contour here is such

46:05

that these are these lead to singularly perturbed problems.

46:09

So when we discretize this elliptic solution operator,

46:13

we use that with the techniques that we've seen before

46:16

and isotropic refinement towards the boundary is there.

46:20

The key thing, and then here,

46:22

for the inhomogeneity, there's another integral in.

46:26

In the time derivative in the time,

46:28

so from zero to T over Tau.

46:30

And yes again and algebraic singularity.

46:33

That can again be dealt with with an HP quadrature

46:37

in the in the towel variable.

46:39

So in other words we have the techniques that we've

46:43

developed into earlier in this talk can actually be combined.

46:48

To get an exponential convergence scheme and this is what

46:51

can actually prove.

46:53

So if the geometry is analytic,

46:55

the initial data is analytic.

46:58

And F hasn't is analytic and we do the appropriate

47:02

mesh refinement.

47:03

So T is refined towards the boundary geometrically.

47:07

And to think portraiture has P squared quadrature points and

47:12

the discretization is done with an HP quadrature,

47:16

then one get one gets exponential convergence in in P.

47:21

And this is a numerical example.

47:26

Where one sees the exponential convergence.

47:30

OK, so let me wrap up.

47:32

I've shown you two discretization techniques.

47:35

The first one was using the increasing the spatial dimension

47:40

by one and using a standard galerkin method on this

47:45

infinite cylinder.

47:46

So when you use the geometric mesh that is refined

47:50

towards Y and on Omega one use a nicer measure

47:53

setter and also tropically refined towards the boundary and in

47:58

polygons one needs to do additional refinement towards the vertices

48:02

of the Polygon and then one gets exponential convergence for

48:07

the fully discrete problem the convergence is root exponential,

48:11

so enter the minus each of them enter the minus

48:15

one 5th for smooth domains.

48:17

And each of the end to the minus 164 polygons

48:20

and the second discretization was we discretize first solution formula

48:25

with a sync monitor and then the resulting decouples singularly

48:30

perturbed problems are solved with high order finite element method

48:35

we use for those the same geometrically refined measures As

48:39

for the previous method and we get the same exponential

48:43

convergence results.

48:45

Of course I should say here that parameter B is

48:49

not the same one,

48:50

but sort of qualitatively one gets the same convergence result,

48:54

and for these use these elliptic solvers can be used

48:58

to.

49:00

You design solvers for parabolic problems,

49:02

so particularly for the fractional heat equation for time fraction

49:06

of the problem,

49:07

and with that I conclude my talk.

49:08

I thank you very much for your attention,

49:10

and here a few more references about the things that

49:13

I talked about today.

49:14

Thank you very much.

49:17

Thank you Marcus for this very very nice talk exactly

49:22

in in the flavor of what I said before bridging

49:26

between analysis and numerical schemes for difficult problems and the

49:31

the presentation is open for questions.

49:42

Anton, you simply. Are you unmuted?

49:46

No, not yet.

49:49

I think you can unmute yourself.

49:51

Is that correct Stanley?

49:55

He says that he doesn't have a mic,

49:56

so I guess we'll have to.

49:59

Just read the question in the in the chat.

50:02

Can you see it focus annoy you?

50:05

Sorry, given that BK and HP fem errors appear pretty

50:09

similar,

50:10

which is more performant in the sense of Wall Clock

50:14

time to solution.

50:18

I think they were fairly similar in Wall Clock time.

50:22

So, because ultimately these are.

50:29

Well, no, they're not dissimilar because one has in the

50:32

for the HP method,

50:34

one has fewer.

50:36

Fewer problems to solve, but of course I mean these

50:39

are the number of linear systems that one has to

50:42

solve.

50:43

I mean, they're completely decoupled,

50:45

so there's this completely can be done in parallel,

50:48

so I it is not as critical as it may

50:50

sound he goes OK.

50:51

I mean, here's the differences between,

50:54

let's say 20 and 40,

50:55

so it's a factor too,

50:56

so.

51:00

So it's not that critical and in the end these

51:03

problems are fairly small.

51:05

Each of these problems is fairly small because we use

51:08

I don't know polynomial degree,

51:10

something like 7 to 10 and the measures are.

51:13

I mean, you see, it's not that many elements.

51:18

And there is a.

51:21

Chat from Squire says please give Boundary equation for L

51:25

shaped colored solution.

51:29

Boundary sorry, can you run this for me again?

51:32

Please give Boundary equation for L shaped colored solution.

51:36

I don't know what this means,

51:38

only conditions or well the boundary conditions are 00.

51:41

That's the same problem we solve.

51:44

It's always dealership conditions.

51:49

So homogeneous directly boundary conditions,

51:52

right? Yeah further questions.

52:00

At the very beginning, when you when you introduced the

52:05

the fractional Laplacian,

52:07

you gave the different possibilities of.

52:11

Of defining it and one was via the eigen decomposition

52:15

of L.

52:16

If you would have a different basis for characterizing let's

52:21

say standard,

52:22

several of spaces like like wavelets or whatever,

52:26

could you then also define the the S power BI

52:30

that replacing the eigendecomposition the eigenfunctions?

52:38

I have not must admit I have not thought about

52:40

that.

52:47

But I will. I wouldn't be too surprised if you

52:51

could do that,

52:52

but I need to think about that because in connection

52:56

with that and the next question would be,

52:59

is there anything known about optimality when you compare the

53:04

error with respect to the best approximation given in degrees

53:08

of freedom?

53:15

I mean you are asking for end with result,

53:17

right? Something like that. Is there something known for those

53:21

kind of problems?

53:22

Yes, there is something now that goes actually goes back

53:25

to an very old paper by myself.

53:31

Well, essentially I these are elliptic problems,

53:35

so you and you have a nice nice size of

53:39

home morphism between.

53:41

If you look at the I mean you have a

53:43

nice isoform morphism between the the energy.

53:46

So the solution and the and the right hand side

53:50

F.

53:50

And so you as soon as you know how well

53:53

you can approximate the right hand side as soon as

53:56

you know it have as soon as you have an

53:59

end with result for the data,

54:01

you can buy this isomorphism,

54:03

create an end with result for the for the solution,

54:06

but of course then the approximation space is kind of

54:11

funny because it will be the.

54:13

Pull back from the from a good space for the

54:17

for the data.

54:18

Yeah, so for example, result is that if the data

54:21

are analytic,

54:22

the end with result gives you exponential convergence with a

54:27

slightly better exponent.

54:29

So we got something like end to the.

54:33

What was it? Enter the.

54:35

1 /, 6 for polygons and you you should for

54:38

polygons you should actually get to enter the minus 1

54:41

/

54:41

2. So there's some room,

54:44

but of course, then these spaces are not pull up

54:47

piecewise polynomial,

54:49

so it's not quite clear what to do OK?

54:53

Good further questions.

55:00

One point I did not quite get.

55:02

I mean, you made it clear that there is a

55:05

difference between the geometry of the domain and and in

55:09

elliptic problems we know that these reentrant corners are are

55:13

an issue. Are there any kind of restrictions to the

55:16

to the shape of the domain?

55:18

Because you said polygonal domain and but it could then

55:23

be more or less arbitrary.

55:25

Yes, well technically I mean.

55:28

Curvilinear polygons are allowed to,

55:30

so that is something that I did not emphasize,

55:33

but yes, that is possible to the slit domain is

55:37

technically speaking not allowed,

55:39

but I'm pretty sure the the theory could be generalized

55:43

to that.

55:43

I mean, you see in a way we went about

55:46

proving this.

55:47

We went about it in a roundabout way because we

55:50

realized we use singular results for singular perturbation problems,

55:54

and this is simply, this is the tool that we

55:58

had so.

55:59

That is also limits here the the result that we

56:02

can prove the initial conditions and two D because that's

56:06

what we use in principle.

56:08

Of course, what we what should hold is that.

56:12

The extension problem should one should be able to write

56:18

down.

56:18

Should characterize my should be able to characterize the solution

56:23

in terms of well,

56:24

these published Google Spaces, so these countably normed spaces,

56:28

when should be able to do that.

56:30

But I mean, this is hard task.

56:32

I mean full proof is not yet given for the

56:35

classical Laplacian.

56:36

I mean we're all waiting,

56:38

and here this isn't. This is,

56:40

we have weights, right? I mean it's not even right.

56:44

Nice right? Elliptic operator? So This is why we took

56:47

the detour.

56:48

Using this decoupling technique.

56:53

So you would need to to consider weight weighted best

56:56

of spaces or whatever.

56:58

So yes, yes guess fun yeah?

57:01

It's feasible it is feasible.

57:04

OK, good further questions. So this does not seem to

57:08

be the case,

57:09

so that let me thank Marcus once again for this

57:12

very nice talk you see also in the in the

57:15

chat:

57:15

was saying fantastic presentations or just that I I forward

57:19

that piece of information to Marco's.

57:21

And then I would hand over the Mike back to

57:24

to Stanley.

57:25

I guess right that about covers it.

57:27

Thanks everybody for coming. Very very interesting talk and please

57:31

join us next week.

57:32

We have a talk by actually my colleague at Florida

57:35

Tech Kenisha Pereira.

57:37

I just put the. Title of the talk in next

57:40

week.

57:40

Talk into the chat. So yes,

57:42

thanks again, Marcus and thanks everyone for coming.

57:45

Thank you, bye. Night.