An asymptotic preserving method for Levy Fokker Planck equation with fractional diffusion limit
Wang, Li
We develop a numerical method for the Levy-Fokker-Planckequation with the fractional diffusive scaling. There are two main challenges. One comes from a two-fold nonlocality, that is, the need to apply the frac-tional Laplacian operator to a power law decay distribution. The other comes from long-time/small mean-free-path scaling, which calls for a uniform stablesolver. To resolve the first difficulty, we use a change of variable to convert theunbounded domain into a bounded one and then apply Chebyshev polyno-mial based pseudo-spectral method. To resolve the second issue, we propose an asymptotic preserving scheme based on a novel micro-macro decomposition that uses the structure of the test function in proving the fractionaldiffusion limit analytically.
 So, today I will talk about asymptotic preserving schemes for Levy-Folker-Plunk equation, and this is a joint work with Ujoshi, who is currently a PhD student at Minnesota. So here's the outline. First I will give a brief introduction on the Levy-Folker-Plunk equation, and then I will introduce our numerical methods and show you some numerical results and draw the conclusion in the end. So this is the Levy-Folker-Plunk equation. As you can see that compared with the classical Folker-Plunk equation, the difference is that it replaces the diffusion with this fractional diffusion. And there are several ways of defining the fractional diffusion, and here I list two of them. So the first one is this integral formulation, or we call it risk formulation. So this Pv is the principal value, and this c is a constant depending on the dimension and this power s. Or you can define it through the Fourier transform. So the Fourier index of this operator is the modulus of k to the power of 2s. So clearly from the second definition, you can see that when s equals 1, it reduces to the classical Folker-Plunk operator, or the classical diffusion operator. And also from the first one, you can see the difference between this one and the classical diffusion is that this one is a non-local operator and the classical one is a local one. So the motivation for considering this equation is for its application in plasma physics. But I believe that the application is far beyond physics. So as long as your microscopic dynamics is governed by a Levy flight instead of Brownian motion, then this kind of operator will appear when you model that system. So first let's look at several properties of this Levy-Folker-Plunk operator. So the first one is the conservation of mass, so that's common, similar, the same as the classical one. But the equilibrium is different. So here the equilibrium instead of being a Gaussian type distribution, it is actually this power law decay. So it's a slow decaying function. And also we can consider the entropy dissipation of this operator. So basically we can consider a very general one. So if you define this entropy function like this, so phi you can consider as a smooth convex function, then we can plug in this f over m so that it can be viewed as the relative entropy, then we have the exponential decay of the relative entropy. So this is given by this Gentile and Invert. So if we want to extract the macroscopic dynamics out of this kinetic equation, we should consider a long time and small mean free path scaling. And it turns out that the correct scaling for this system is the following. So you can see that instead of rescale t to epsilon square times t, we have to consider this fractional power. And so the reason for this fractional power is to, is the following. So if you consider the classical case where we have this classical diffusion, then immediately we see that when epsilon goes to zero, f will converge to this local equilibrium. And if you integrate the whole equation in v, then you get the conservation of mass and with this flux J. And to determine this flux, you just multiply v to the equation and integrate again. And this is what you will get. And sending epsilon to zero, then you'll see that this J over epsilon to the leading order is this term. So plug it back into the conservation equation, you end up with this diffusion equation with the diffusion coefficient or diffusion matrix of this form. So from here, you can see that if we, so here this one is the equilibrium for the classical case. But if you want to plug in the equilibrium for the fractional Laplacian case, then this diffusion matrix becomes infinite. So that's why you cannot consider the classical diffusive scaling. So with this anonymous scaling, we have the following results given by Cessbara-Mellet and Trivissa. So it says that as epsilon goes to zero, f will converge to this local equilibrium, which is rho times m, in the following week topology. And this rho satisfies the fractional Laplacian equation. So you can clearly see that the fractional Laplacian or fractional diffusion in v in the original kinetic equation transfers to the fractional Laplacian in rho in the limit. And the way to see how we can get this limit is actually quite simple. So again, first, we know that when we send up to zero, f will converge to the local equilibrium. So here is the same as before, except that the local equilibrium here is this power law decay function. And then to derive the equation for this density rho, they have to construct this special test function. So basically, instead of allowing these x and v to be two independent variables, they require that these x and v are related in the following way. So once you plug in this special form of test function, so you multiply this test function to the original equation and you integrate in x, v, and t, immediately you see that because of the special choice of this test function, these two terms cancels. And this fractional diffusion in v becomes the fractional diffusion in x. So you can simplify this equation in the following way. And then sending epsilon to zero, it will converge to the equation for rho. And you see that this equation is nothing but the weak formulation of the fractional diffusion for the density rho. All right, so this is how we derive the fractional diffusion limit. Now let's talk about numerics. So our goal is to compute this Levy-Foucault-Plunk equation and hopefully that when epsilon goes to zero, it will become a solver for this fractional diffusion equation. And there are two main challenges. So one is the non-locality that comes from computing this Levy-Foucault-Plunk operator itself. So here because of the interplay of these two components, we know that this f will converge to the local equilibrium, which is a slow decaying function. So that means in principle we need to compute the fractional laplacian on a slow decaying function. And actually computing this non-local operator is an active area of research and there has been several quite good results in the past few years. So usually the method to compute this fractional laplacian or non-local operator in general can be categorized into these two classes. So the first one is this finite difference and finite volume method. So the idea is that you have to truncate your domain and either the function is fast decaying so that you can ignore the tail or if the function is slow decaying and the tail information is important, what you can do is that you add the tail information back, but that highly relies on the analytical behavior of the tail. So which is not feasible for our case because in our case that unless f reaches the equilibrium or unless it is at the initial condition, we don't know how the tail behaves. Also people also consider adding some artificial boundary condition when you do the truncation, but that would generate instability. So then the spectral method has been developed and it is more applicable. So the idea of the spectral method is it uses the non-local basis or in this case is the orthogonal polynomials and with the hope that if you apply the fractional laplacian to the basis function then it is easy or at least not hard to compute. So for instance Mao and Shen has developed this Hermit polynomial basis based spectral method and the idea there is that since the Hermit polynomial is invariant under Fourier transform so you can see that if you want to compute the fractional laplacian of a Hermit polynomial you can just transfer to the Fourier space and apply this Fourier symbol and then transfer back. But it turns out that if we want to use this method to compute a slow decaying function we need a lot of modes. So later on people realize that the Chebyshev polynomial is more performs well, performs better when you want to approximate a slow decaying function. So there are these two works that are based on this Chebyshev polynomial based spectral method and so these two work just treats the fractional laplacian applying on Chebyshev polynomial in a different way. And also you can refer to these two review papers. So here we will adopt approach in this work but I want to mention that our choice is not unique. Actually there are still room to improve. So the second thing is the stiffness introduced by this epsilon and for that we aim to develop the asymptotic preserving scheme. So here I listed a few results that for the linear Boltzmann equation in which case we still have the fractional diffusion but I want to emphasize that here the major difficulty is we do not have a strong convergence or we do not have a Hilbert expansion based convergence. So we know that Hilbert expansion actually gives us a lot of guidance when we design numerical methods because it gives us an idea of what's the magnitude of each term and how should we group the terms. But as you can see that when we derive in the derivation of fractional diffusion limit here we use this tricky choice of test function which gives us limited information on how we can construct the numerical methods. Alright so let me first show you how we compute this fractional laplacian and as I mentioned this is adopted from this work. So it goes as the following steps. So first the first step is we do a change of variable so this V is in the whole domain and if we do the change of variable not only a change of variable in the following way then c is in between minus 1 and 1 and this Lv you can consider that's a scaling parameter. So basically the larger the Lv you can see that the points are more spread out so usually if you want to compute a slow decay function then you want to choose this Lv to be large. But of course if you choose Lv to be large then there will be less points in the body so you have to strike a balance there. And then once we are in the finite domain from minus 1 to 1 then we can consider the Chebyshev polynomials and then if we do further change of variable if we let this arc cos i can see to be q then the Chebyshev polynomial just reduces to this trigonometric function. And now if we just work with the new variable q and keep in mind that the relation between q and original variable V is the following then we have to rewrite the fractional laplacian in this new variable q which looks a little bit complicated but the bottom line is it is an analytical expression. Now if we want to apply this fractional laplacian to a function so remember now the function if we write it into the new variable q then q is in the domain from 0 to pi so we can do an even extension and extend the function to the domain from 0 to 2 pi and in this way if we were to expand this f tilde q in the basis which is the cosine kq and because this is an even function so that expansion turns out to be quite easy it's just the discrete Fourier transform. And then the next thing is that we have to apply the fractional laplacian to this f of q and then it reduces to compute the fractional laplacian of the basis function and in this case the basis function is just this Fourier basis. But the thing is this fractional laplacian is no longer the original one so you cannot simply say that this one equals the modulus of k to the power of 2s times the Fourier basis instead we have to resort to the previous new formulation of this fractional laplacian. So actually we have so they have derived this analytical formula for this fractional laplacian applying on the Fourier basis but here I'm not going to list it in these slides it will be tedious but it can be pre-computed. Alright so let me show you the performance of this method by applying it to two kinds of f so the left hand side is this exponential decay function and the right hand side is this power law decay function. So for these two choices we have the analytical formula where we apply the fractional laplacian and this is the error between the numerical result and the analytical formula. So basically you can tell that the slow decay function is much harder to compute because with the same number of modes on the accuracy for this one is much lower than this one so you can see that if you want to reach the same accuracy then you will need a lot more modes for the slow decay function which is not surprising. Okay so if now we go back to the Levy-Folker Planck operator we have another term but this term is easy so we just rewrite it in terms of the new variable Q. Okay so now that we know how to compute the Levy-Folker Planck operator we have to deal with the stiffness and our idea is still we want to use this micro-micro decomposition but the key idea is that instead of saying that the macro part is rho times m we allow the macro part to depend on v as well but not in an arbitrary way. So although we have this dependence on v so the dependence is very mild so we require this x and v actually has the following relation. So this is inspired by this weak formulation. So once we have this special form of decomposition and we plug it into the original equation and this is what we get and again due to the special form of this eta we see that these two red term cancelled and the blue term changes to the fractional aplasia in x and this term goes to zero because m is the equilibrium and this i is just you can just consider as the residual where you try to apply sort of product rule where you apply this fractional aplasia. Okay so after this simplification the equation becomes the following so this equation is completely equivalent to the original one but then we do a splitting. So we split it into two equations one is for eta for the macro part and one is for the micro part. So the rationale behind this splitting is that we just keep in mind that we want to what we want is the in the end is the fractional aplasia and you can see that this eta the eta term in the original equation they are of the same order and they almost gives us what the diffusion equation should look like and then the rest of the thing is put into the micro part. So you can see that evolution is that eta evolves by itself but the row evolves according to eta and the initial condition is decomposed in the following way so eta is of course related to row but in the in the special form and and then the rest of the thing is given to g. So there is a slight thing we should be careful is that the whole derivation here relies on the fact that this eta is of this special form so we but we know that the but the evolution that eta obeys is the following equation so we need to make sure that eta that the the form of eta is preserved along the dynamics which is not hard to show. Alright so so we check a little bit about this splitting system so we have to make sure that this splitting is not too wild so we check the energy stability of this splitting and it turns out that the energy for the f so that which is they reconstructed or the total or just the original function f and also the energy for the macro and the micro part individually are uniformly bounded in time. Okay now the next thing is we need to solve this splitting system and if you look at the equation for g you will you will notice that it is as hard to solve as the original equation for f so our idea is we do an operator splitting so basically we split the the convection part from the collision part so in the first stage of the splitting we just focus on the collision part and then in the second stage we focus on the convection part but we also do a auxiliary term so this is like subtract add and subtract and so we add we subtract this gamma g in the first splitting and add it back in the second splitting. So for this splitting method compared to the original one the the advantage is quite obvious so first of all because of this auxiliary term it alleviates the eoconditioning so because we know that this L-S has this non-empty non-space so if you were to invert the system when epsilon is very small then the system to be inverted will be very stiff so actually it can be shown from this condition number so you can see that when epsilon is very small the condition number just grow dramatically but with this additional term the condition number just grow mildly so that's one advantage and the second advantage is to reduce the computational cost so here instead of just because so if you consider the original equation and if you treat it them implicitly then you have to deal with x and y simultaneously so that means you have to invert a very large system but once we do this the operator splitting then in the first stage we just need to do to invert an equation only for v so so x we can treat it locally and the same for the second one we just treat v locally so the computational cost is largely reduced and the third one is for the asymptotic property concern which I will explain in the next slides so basically we have the following result that the numerical solution obtained by this reconstruction will satisfies the following equation when epsilon goes to zero and to see why this is true first let's from this reconstructed solution we can see that so the first one this eta because the eta satisfies the fractional diffusion so it can be replaced in the following way and then with a pretty cheap argument we can swap the average with this fractional equation so basically we can move this average inside and then the average or the integral of f is rho and so you can see that this is pretty much what we want with two additional terms depending on g so if you can show that g or the average of g is a small term then we are good so we resort to the equation for g so here the first from the first stage we just need to control the magnitude of g so as long as it is bounded and then from the second stage it is clearly that it will be a small term so it will be order epsilon to the power of 2s so from this derivation you can also see the advantage of doing this operator splitting is that it gives us a little freedom so for the first stage we don't require so we don't need g to be small just immediately so for the first stage it just needs to be bounded but then the second stage gives the power of driving g to a small magnitude so that's why in the end it converges to the equation for the fractional flash okay so let me show you some numerical results so the first one is the spatially homogeneous case and here we check the long time behavior so here when s equals 0.5 we have the analytical form of the equilibrium so here the left hand side is we compare the numerical solution with this analytical formula and here we can only plot the body part so you can see that the body part they match very well and so it is with the tail part and then the last figure that's for the exponential decay of the relative entropy and we also check two other cases where s equals 0.6 and 0.8 and in these two cases we do not have the analytical formula to compare with so we only see whether we obtain the correct tail behavior and as you can see that the tail is proportional to this power and the next is for the spatially inhomogeneous case and you can see that we have the uniform convergence for a different choices of epsilon and these two figures are for different s and this is to check the energy stability and again this is for different choices of epsilon so for the last one you can see that there is a sharp transition this is due to the presence of initial layer and also we compare the solution with the solution to the original kinetic equation so we use a IMAX scheme with a very fine mesh and when epsilon equals 1 this is the comparison for our solution with the solution obtained by this IMAX scheme and the last one is to show the comparison with the diffusion limit and this one is to show the asymptotic preserving property basically if you decrease your epsilon then the difference between your f and the local equilibrium should be of other epsilon. Alright so to conclude that we have developed this asymptotic preserving scheme and the key idea is the spectral method with no local basis that is to treat the V-focal Planck operator and to design this asymptotic preserving scheme we need to use this micro micro decomposition and also to solve the split system we need to use the operator splitting so with that I'm going to end my talk and thanks for your attention. Thanks for your nice work and talk so you mentioned the so Cess Brung and Menon and Trevisa and several years ago proved this fractional diffusion limit and if you fact about the two years or three years ago Cess Brung and Menon generalized their work to the half space with non-local boundary condition so my question is is there any possible and to I mean numerically capture this non-local boundary behavior with this APS scheme? That's my question. We haven't tried that yet but yes I believe that we I don't know if this method can directly apply but I believe that we can develop sort of this APS scheme for that kind for the equation with boundary yes. Yeah okay thank you so because this boundary condition is very interesting because the boundary condition is kind of a non-local Neumann boundary condition so yeah. Yeah yeah yeah so I know that when you have fractional aplasia and you have the boundary domain then yes the boundary condition is highly non-trivial. Okay yeah thanks. How is your method sensitive to the operator you choose? Which operator? Yeah which operator? Maybe you'll have a message on chat you cannot say it? If it's not a log... No I see. If it is not well so I have to say that the way we do this micromicro decomposition that is specifically designed for the passage from Levy-Folker Planck to the fractional diffusion. So if it is not this operator then the mechanism for driving to the fractional diffusion might be different than the design of this micromicro decomposition will need to change so in that case I would say that yes it is specific to the Levy-Folker Planck operator. He says that you can unite the working teams to create the work. Yeah, significant difference there. Can you do if you are the type of people of the ATC in the world like aurance buildings
 So, today I will talk about asymptotic preserving schemes for Levy-Folker-Plunk equation, and  this is a joint work with Ujoshi, who is currently a PhD student at Minnesota.  So here's the outline. First I will give a brief introduction on the Levy-Folker-Plunk  equation, and then I will introduce our numerical methods and show you some numerical results  and draw the conclusion in the end. So this is the Levy-Folker-Plunk equation. As you  can see that compared with the classical Folker-Plunk equation, the difference is that it replaces  the diffusion with this fractional diffusion. And there are several ways of defining the  fractional diffusion, and here I list two of them. So the first one is this integral  formulation, or we call it risk formulation. So this Pv is the principal value, and this  c is a constant depending on the dimension and this power s. Or you can define it through  the Fourier transform. So the Fourier index of this operator is the modulus of k to the  power of 2s. So clearly from the second definition, you can see that when s equals 1, it reduces  to the classical Folker-Plunk operator, or the classical diffusion operator. And also  from the first one, you can see the difference between this one and the classical diffusion  is that this one is a non-local operator and the classical one is a local one. So the motivation  for considering this equation is for its application in plasma physics. But I believe that the  application is far beyond physics. So as long as your microscopic dynamics is governed by  a Levy flight instead of Brownian motion, then this kind of operator will appear when  you model that system. So first let's look at several properties of this Levy-Folker-Plunk  operator. So the first one is the conservation of mass, so that's common, similar, the same  as the classical one. But the equilibrium is different. So here the equilibrium instead  of being a Gaussian type distribution, it is actually this power law decay. So it's  a slow decaying function. And also we can consider the entropy dissipation of this operator.  So basically we can consider a very general one. So if you define this entropy function  like this, so phi you can consider as a smooth convex function, then we can plug in this  f over m so that it can be viewed as the relative entropy, then we have the exponential decay  of the relative entropy. So this is given by this Gentile and Invert.  So if we want to extract the macroscopic dynamics out of this kinetic equation, we should consider  a long time and small mean free path scaling. And it turns out that the correct scaling  for this system is the following. So you can see that instead of rescale t to epsilon square  times t, we have to consider this fractional power. And so the reason for this fractional  power is to, is the following. So if you consider the classical case where we have this classical  diffusion, then immediately we see that when epsilon goes to zero, f will converge to  this local equilibrium. And if you integrate the whole equation in v, then you get the  conservation of mass and with this flux J. And to determine this flux, you just multiply  v to the equation and integrate again. And this is what you will get. And sending epsilon  to zero, then you'll see that this J over epsilon to the leading order is this term.  So plug it back into the conservation equation, you end up with this diffusion equation with  the diffusion coefficient or diffusion matrix of this form. So from here, you can see that  if we, so here this one is the equilibrium for the classical case. But if you want to  plug in the equilibrium for the fractional Laplacian case, then this diffusion matrix  becomes infinite. So that's why you cannot consider the classical diffusive scaling.  So with this anonymous scaling, we have the following results given by Cessbara-Mellet  and Trivissa. So it says that as epsilon goes to zero, f will converge to this local equilibrium,  which is rho times m, in the following week topology. And this rho satisfies the fractional  Laplacian equation. So you can clearly see that the fractional Laplacian or fractional  diffusion in v in the original kinetic equation transfers to the fractional Laplacian in rho  in the limit. And the way to see how we can get this limit is actually quite simple. So  again, first, we know that when we send up to zero, f will converge to the local equilibrium.  So here is the same as before, except that the local equilibrium here is this power law  decay function. And then to derive the equation for this density rho, they have to construct  this special test function. So basically, instead of allowing these x and v to be two  independent variables, they require that these x and v are related in the following way. So  once you plug in this special form of test function, so you multiply this test function  to the original equation and you integrate in x, v, and t, immediately you see that because  of the special choice of this test function, these two terms cancels. And this fractional  diffusion in v becomes the fractional diffusion in x. So you can simplify this equation in  the following way. And then sending epsilon to zero, it will converge to the equation  for rho. And you see that this equation is nothing but the weak formulation of the fractional  diffusion for the density rho. All right, so this is how we derive the fractional diffusion  limit. Now let's talk about numerics. So our goal is to compute this Levy-Foucault-Plunk  equation and hopefully that when epsilon goes to zero, it will become a solver for this  fractional diffusion equation. And there are two main challenges. So one is the non-locality  that comes from computing this Levy-Foucault-Plunk operator itself. So here because of the interplay  of these two components, we know that this f will converge to the local equilibrium, which  is a slow decaying function. So that means in principle we need to compute the fractional  laplacian on a slow decaying function. And actually computing this non-local operator  is an active area of research and there has been several quite good results in the past  few years. So usually the method to compute this fractional laplacian or non-local operator  in general can be categorized into these two classes. So the first one is this finite difference  and finite volume method. So the idea is that you have to truncate your domain and either  the function is fast decaying so that you can ignore the tail or if the function is  slow decaying and the tail information is important, what you can do is that you add  the tail information back, but that highly relies on the analytical behavior of the tail.  So which is not feasible for our case because in our case that unless f reaches the equilibrium  or unless it is at the initial condition, we don't know how the tail behaves. Also people  also consider adding some artificial boundary condition when you do the truncation, but  that would generate instability. So then the spectral method has been developed and it  is more applicable. So the idea of the spectral method is it uses the non-local basis or in  this case is the orthogonal polynomials and with the hope that if you apply the fractional  laplacian to the basis function then it is easy or at least not hard to compute. So for  instance Mao and Shen has developed this Hermit polynomial basis based spectral method  and the idea there is that since the Hermit polynomial is invariant under Fourier transform  so you can see that if you want to compute the fractional laplacian of a Hermit polynomial  you can just transfer to the Fourier space and apply this Fourier symbol and then transfer  back. But it turns out that if we want to use this method to compute a slow decaying  function we need a lot of modes. So later on people realize that the Chebyshev polynomial  is more performs well, performs better when you want to approximate a slow decaying function.  So there are these two works that are based on this Chebyshev polynomial based spectral  method and so these two work just treats the fractional laplacian applying on Chebyshev  polynomial in a different way. And also you can refer to these two review papers. So here  we will adopt approach in this work but I want to mention that our choice is not unique.  Actually there are still room to improve. So the second thing is the stiffness introduced  by this epsilon and for that we aim to develop the asymptotic preserving scheme. So here I  listed a few results that for the linear Boltzmann equation in which case we still have the fractional  diffusion but I want to emphasize that here the major difficulty is we do not have a strong  convergence or we do not have a Hilbert expansion based convergence. So we know that Hilbert  expansion actually gives us a lot of guidance when we design numerical methods because it  gives us an idea of what's the magnitude of each term and how should we group the terms.  But as you can see that when we derive in the derivation of fractional diffusion limit  here we use this tricky choice of test function which gives us limited information on how we  can construct the numerical methods. Alright so let me first show you how we compute this  fractional laplacian and as I mentioned this is adopted from this work. So it goes as the  following steps. So first the first step is we do a change of variable so this V is in  the whole domain and if we do the change of variable not only a change of variable in  the following way then c is in between minus 1 and 1 and this Lv you can consider that's  a scaling parameter. So basically the larger the Lv you can see that the points are more  spread out so usually if you want to compute a slow decay function then you want to choose  this Lv to be large. But of course if you choose Lv to be large then there will be less points  in the body so you have to strike a balance there. And then once we are in the finite  domain from minus 1 to 1 then we can consider the Chebyshev polynomials and then if we do  further change of variable if we let this arc cos i can see to be q then the Chebyshev  polynomial just reduces to this trigonometric function. And now if we just work with the  new variable q and keep in mind that the relation between q and original variable V is the following  then we have to rewrite the fractional laplacian in this new variable q which looks a little  bit complicated but the bottom line is it is an analytical expression.  Now if we want to apply this fractional laplacian to a function so remember now the function  if we write it into the new variable q then q is in the domain from 0 to pi so we can do  an even extension and extend the function to the domain from 0 to 2 pi and in this way  if we were to expand this f tilde q in the basis which is the cosine kq and because this  is an even function so that expansion turns out to be quite easy it's just the discrete  Fourier transform. And then the next thing is that we have to apply the fractional laplacian  to this f of q and then it reduces to compute the fractional laplacian of the basis function  and in this case the basis function is just this Fourier basis. But the thing is this  fractional laplacian is no longer the original one so you cannot simply say that this one  equals the modulus of k to the power of 2s times the Fourier basis instead we have to  resort to the previous new formulation of this fractional laplacian. So actually we have  so they have derived this analytical formula for this fractional laplacian applying on  the Fourier basis but here I'm not going to list it in these slides it will be tedious  but it can be pre-computed. Alright so let me show you the performance of this method  by applying it to two kinds of f so the left hand side is this exponential decay function  and the right hand side is this power law decay function. So for these two choices we have  the analytical formula where we apply the fractional laplacian and this is the error  between the numerical result and the analytical formula. So basically you can tell that the  slow decay function is much harder to compute because with the same number of modes on the  accuracy for this one is much lower than this one so you can see that if you want to reach  the same accuracy then you will need a lot more modes for the slow decay function which  is not surprising. Okay so if now we go back to the Levy-Folker  Planck operator we have another term but this term is easy so we just rewrite it in terms  of the new variable Q. Okay so now that we know how to compute the Levy-Folker Planck  operator we have to deal with the stiffness and our idea is still we want to use this  micro-micro decomposition but the key idea is that instead of saying that the macro part  is rho times m we allow the macro part to depend on v as well but not in an arbitrary  way. So although we have this dependence on v so the dependence is very mild so we require  this x and v actually has the following relation. So this is inspired by this weak formulation.  So once we have this special form of decomposition and we plug it into the original equation  and this is what we get and again due to the special form of this eta we see that these  two red term cancelled and the blue term changes to the fractional aplasia in x and this term  goes to zero because m is the equilibrium and this i is just you can just consider as  the residual where you try to apply sort of product rule where you apply this fractional  aplasia. Okay so after this simplification the equation becomes the following so this  equation is completely equivalent to the original one but then we do a splitting. So we split  it into two equations one is for eta for the macro part and one is for the micro part.  So the rationale behind this splitting is that we just keep in mind that we want to what  we want is the in the end is the fractional aplasia and you can see that this eta the  eta term in the original equation they are of the same order and they almost gives us  what the diffusion equation should look like and then the rest of the thing is put into  the micro part. So you can see that evolution is that eta evolves by itself but the row  evolves according to eta and the initial condition is decomposed in the following way so eta  is of course related to row but in the in the special form and and then the rest of the  thing is given to g. So there is a slight thing we should be careful is that the whole  derivation here relies on the fact that this eta is of this special form so we but we know  that the but the evolution that eta obeys is the following equation so we need to make  sure that eta that the the form of eta is preserved along the dynamics which is not  hard to show. Alright so so we check a little bit about this splitting system so we have  to make sure that this splitting is not too wild so we check the energy stability of this  splitting and it turns out that the energy for the f so that which is they reconstructed  or the total or just the original function f and also the energy for the macro and the  micro part individually are uniformly bounded in time. Okay now the next thing is we need  to solve this splitting system and if you look at the equation for g you will you will  notice that it is as hard to solve as the original equation for f so our idea is we do  an operator splitting so basically we split the the convection part from the collision  part so in the first stage of the splitting we just focus on the collision part and then  in the second stage we focus on the convection part but we also do a auxiliary term so this  is like subtract add and subtract and so we add we subtract this gamma g in the first  splitting and add it back in the second splitting. So for this splitting method compared to the  original one the the advantage is quite obvious so first of all because of this auxiliary  term it alleviates the eoconditioning so because we know that this L-S has this non-empty  non-space so if you were to invert the system when epsilon is very small then the system  to be inverted will be very stiff so actually it can be shown from this condition number  so you can see that when epsilon is very small the condition number just grow dramatically  but with this additional term the condition number just grow mildly so that's one advantage  and the second advantage is to reduce the computational cost so here instead of just  because so if you consider the original equation and if you treat it them implicitly then you  have to deal with x and y simultaneously so that means you have to invert a very large  system but once we do this the operator splitting then in the first stage we just need to do  to invert an equation only for v so so x we can treat it locally and the same for the  second one we just treat v locally so the computational cost is largely reduced and  the third one is for the asymptotic property concern which I will explain in the next slides  so basically we have the following result that the numerical solution obtained by this  reconstruction will satisfies the following equation when epsilon goes to zero and to  see why this is true first let's from this reconstructed solution we can see that so  the first one this eta because the eta satisfies the fractional diffusion so it can be replaced  in the following way and then with a pretty cheap argument we can swap the average with  this fractional equation so basically we can move this average inside and then the average  or the integral of f is rho and so you can see that this is pretty much what we want  with two additional terms depending on g so if you can show that g or the average of g  is a small term then we are good so we resort to the equation for g so here the first from  the first stage we just need to control the magnitude of g so as long as it is bounded  and then from the second stage it is clearly that it will be a small term so it will be  order epsilon to the power of 2s so from this derivation you can also see the advantage  of doing this operator splitting is that it gives us a little freedom so for the first  stage we don't require so we don't need g to be small just immediately so for the first  stage it just needs to be bounded but then the second stage gives the power of driving  g to a small magnitude so that's why in the end it converges to the equation for the fractional  flash okay so let me show you some numerical results so the first one is the spatially  homogeneous case and here we check the long time behavior so here when s equals 0.5 we  have the analytical form of the equilibrium so here the left hand side is we compare the  numerical solution with this analytical formula and here we can only plot the body part so  you can see that the body part they match very well and so it is with the tail part  and then the last figure that's for the exponential decay of the relative entropy and we also  check two other cases where s equals 0.6 and 0.8 and in these two cases we do not have  the analytical formula to compare with so we only see whether we obtain the correct tail  behavior and as you can see that the tail is proportional to this power and the next  is for the spatially inhomogeneous case and you can see that we have the uniform convergence  for a different choices of epsilon and these two figures are for different s and this is  to check the energy stability and again this is for different choices of epsilon so for  the last one you can see that there is a sharp transition this is due to the presence of  initial layer and also we compare the solution with the solution to the original kinetic  equation so we use a IMAX scheme with a very fine mesh and when epsilon equals 1 this is  the comparison for our solution with the solution obtained by this IMAX scheme and the last  one is to show the comparison with the diffusion limit and this one is to show the asymptotic  preserving property basically if you decrease your epsilon then the difference between your  f and the local equilibrium should be of other epsilon.  Alright so to conclude that we have developed this asymptotic preserving scheme and the key  idea is the spectral method with no local basis that is to treat the V-focal Planck  operator and to design this asymptotic preserving scheme we need to use this micro micro decomposition  and also to solve the split system we need to use the operator splitting so with that  I'm going to end my talk and thanks for your attention.  Thanks for your nice work and talk so you mentioned the so Cess Brung and Menon and  Trevisa and several years ago proved this fractional diffusion limit and if you fact  about the two years or three years ago Cess Brung and Menon generalized their work to  the half space with non-local boundary condition so my question is is there any possible and  to I mean numerically capture this non-local boundary behavior with this APS scheme?  That's my question.  We haven't tried that yet but yes I believe that we I don't know if this method can directly  apply but I believe that we can develop sort of this APS scheme for that kind for the equation  with boundary yes.  Yeah okay thank you so because this boundary condition is very interesting because the  boundary condition is kind of a non-local Neumann boundary condition so yeah.  Yeah yeah yeah so I know that when you have fractional aplasia and you have the boundary  domain then yes the boundary condition is highly non-trivial.  Okay yeah thanks.  How is your method sensitive to the operator you choose?  Which operator?  Yeah which operator?  Maybe you'll have a message on chat you cannot say it?  If it's not a log...  No I see.  If it is not well so I have to say that the way we do this micromicro decomposition that  is specifically designed for the passage from Levy-Folker Planck to the fractional diffusion.  So if it is not this operator then the mechanism for driving to the fractional diffusion might  be different than the design of this micromicro decomposition will need to change so in that  case I would say that yes it is specific to the Levy-Folker Planck operator.  He says that you can unite the working teams to create the work.  Yeah, significant difference there.  Can you do if you are the type of people of the ATC in the world like aurance buildings
2021-03-23_Wang_0_47.png 2021-03-23_Wang_48_239.png 2021-03-23_Wang_240_323.png 2021-03-23_Wang_324_659.png 2021-03-23_Wang_660_1067.png 2021-03-23_Wang_1068_3671.png 2021-03-23_Wang_3672_5279.png 2021-03-23_Wang_5280_6083.png 2021-03-23_Wang_6084_8255.png 2021-03-23_Wang_8256_9407.png 2021-03-23_Wang_9408_10019.png 2021-03-23_Wang_10020_10799.png 2021-03-23_Wang_10800_11891.png 2021-03-23_Wang_11892_11987.png 2021-03-23_Wang_11988_13319.png 2021-03-23_Wang_13320_18167.png 2021-03-23_Wang_18168_19955.png 2021-03-23_Wang_19956_21347.png 2021-03-23_Wang_21348_21599.png 2021-03-23_Wang_21600_22487.png 2021-03-23_Wang_22488_23507.png 2021-03-23_Wang_23508_25043.png 2021-03-23_Wang_25044_26879.png 2021-03-23_Wang_26880_28319.png 2021-03-23_Wang_28320_28559.png 2021-03-23_Wang_28560_29315.png 2021-03-23_Wang_29316_29675.png 2021-03-23_Wang_29676_31307.png 2021-03-23_Wang_31308_32039.png 2021-03-23_Wang_32040_32927.png 2021-03-23_Wang_32928_35579.png 2021-03-23_Wang_35580_36755.png 2021-03-23_Wang_36756_37211.png 2021-03-23_Wang_37212_38483.png 2021-03-23_Wang_38484_39983.png 2021-03-23_Wang_39984_40031.png 2021-03-23_Wang_40032_41027.png 2021-03-23_Wang_41028_41603.png 2021-03-23_Wang_41604_42035.png 2021-03-23_Wang_42036_42479.png 2021-03-23_Wang_42480_43067.png 2021-03-23_Wang_43068_43523.png 2021-03-23_Wang_43524_43571.png 2021-03-23_Wang_43572_44279.png 2021-03-23_Wang_44280_48347.png 2021-03-23_Wang_48348_48431.png 2021-03-23_Wang_48432_48479.png 2021-03-23_Wang_48480_48527.png 2021-03-23_Wang_48528_48593.png
So, today I will talk about asymptotic preserving schemes for Levy-Folker-Plunk equation, and this is a joint work with Ujoshi, who is currently a PhD student at Minnesota. So here's the outline. First I will give a brief introduction on the Levy-Folker-Plunk equation, and then I will introduce our numerical methods and show you some numerical results and draw the conclusion in the end. So this is the Levy-Folker-Plunk equation. As you can see that compared with the classical Folker-Plunk equation, the difference is that it replaces the diffusion with this fractional diffusion. And there are several ways of defining the fractional diffusion, and here I list two of them. So the first one is this integral formulation, or we call it risk formulation. So this Pv is the principal value, and this c is a constant depending on the dimension and this power s. Or you can define it through the Fourier transform. So the Fourier index of this operator is the modulus of k to the power of 2s. So clearly from the second definition, you can see that when s equals 1, it reduces to the classical Folker-Plunk operator, or the classical diffusion operator. And also from the first one, you can see the difference between this one and the classical diffusion is that this one is a non-local operator and the classical one is a local one. So the motivation for considering this equation is for its application in plasma physics. But I believe that the application is far beyond physics. So as long as your microscopic dynamics is governed by a Levy flight instead of Brownian motion, then this kind of operator will appear when you model that system. So first let's look at several properties of this Levy-Folker-Plunk operator. So the first one is the conservation of mass, so that's common, similar, the same as the classical one. But the equilibrium is different. So here the equilibrium instead of being a Gaussian type distribution, it is actually this power law decay. So it's a slow decaying function. And also we can consider the entropy dissipation of this operator. So basically we can consider a very general one. So if you define this entropy function like this, so phi you can consider as a smooth convex function, then we can plug in this f over m so that it can be viewed as the relative entropy, then we have the exponential decay of the relative entropy. So this is given by this Gentile and Invert. So if we want to extract the macroscopic dynamics out of this kinetic equation, we should consider a long time and small mean free path scaling. And it turns out that the correct scaling for this system is the following. So you can see that instead of rescale t to epsilon square times t, we have to consider this fractional power. And so the reason for this fractional power is to, is the following. So if you consider the classical case where we have this classical diffusion, then immediately we see that when epsilon goes to zero, f will converge to this local equilibrium. And if you integrate the whole equation in v, then you get the conservation of mass and with this flux J. And to determine this flux, you just multiply v to the equation and integrate again. And this is what you will get. And sending epsilon to zero, then you'll see that this J over epsilon to the leading order is this term. So plug it back into the conservation equation, you end up with this diffusion equation with the diffusion coefficient or diffusion matrix of this form. So from here, you can see that if we, so here this one is the equilibrium for the classical case. But if you want to plug in the equilibrium for the fractional Laplacian case, then this diffusion matrix becomes infinite. So that's why you cannot consider the classical diffusive scaling. So with this anonymous scaling, we have the following results given by Cessbara-Mellet and Trivissa. So it says that as epsilon goes to zero, f will converge to this local equilibrium, which is rho times m, in the following week topology. And this rho satisfies the fractional Laplacian equation. So you can clearly see that the fractional Laplacian or fractional diffusion in v in the original kinetic equation transfers to the fractional Laplacian in rho in the limit. And the way to see how we can get this limit is actually quite simple. So again, first, we know that when we send up to zero, f will converge to the local equilibrium. So here is the same as before, except that the local equilibrium here is this power law decay function. And then to derive the equation for this density rho, they have to construct this special test function. So basically, instead of allowing these x and v to be two independent variables, they require that these x and v are related in the following way. So once you plug in this special form of test function, so you multiply this test function to the original equation and you integrate in x, v, and t, immediately you see that because of the special choice of this test function, these two terms cancels. And this fractional diffusion in v becomes the fractional diffusion in x. So you can simplify this equation in the following way. And then sending epsilon to zero, it will converge to the equation for rho. And you see that this equation is nothing but the weak formulation of the fractional diffusion for the density rho. All right, so this is how we derive the fractional diffusion limit. Now let's talk about numerics. So our goal is to compute this Levy-Foucault-Plunk equation and hopefully that when epsilon goes to zero, it will become a solver for this fractional diffusion equation. And there are two main challenges. So one is the non-locality that comes from computing this Levy-Foucault-Plunk operator itself. So here because of the interplay of these two components, we know that this f will converge to the local equilibrium, which is a slow decaying function. So that means in principle we need to compute the fractional laplacian on a slow decaying function. And actually computing this non-local operator is an active area of research and there has been several quite good results in the past few years. So usually the method to compute this fractional laplacian or non-local operator in general can be categorized into these two classes. So the first one is this finite difference and finite volume method. So the idea is that you have to truncate your domain and either the function is fast decaying so that you can ignore the tail or if the function is slow decaying and the tail information is important, what you can do is that you add the tail information back, but that highly relies on the analytical behavior of the tail. So which is not feasible for our case because in our case that unless f reaches the equilibrium or unless it is at the initial condition, we don't know how the tail behaves. Also people also consider adding some artificial boundary condition when you do the truncation, but that would generate instability. So then the spectral method has been developed and it is more applicable. So the idea of the spectral method is it uses the non-local basis or in this case is the orthogonal polynomials and with the hope that if you apply the fractional laplacian to the basis function then it is easy or at least not hard to compute. So for instance Mao and Shen has developed this Hermit polynomial basis based spectral method and the idea there is that since the Hermit polynomial is invariant under Fourier transform so you can see that if you want to compute the fractional laplacian of a Hermit polynomial you can just transfer to the Fourier space and apply this Fourier symbol and then transfer back. But it turns out that if we want to use this method to compute a slow decaying function we need a lot of modes. So later on people realize that the Chebyshev polynomial is more performs well, performs better when you want to approximate a slow decaying function. So there are these two works that are based on this Chebyshev polynomial based spectral method and so these two work just treats the fractional laplacian applying on Chebyshev polynomial in a different way. And also you can refer to these two review papers. So here we will adopt approach in this work but I want to mention that our choice is not unique. Actually there are still room to improve. So the second thing is the stiffness introduced by this epsilon and for that we aim to develop the asymptotic preserving scheme. So here I listed a few results that for the linear Boltzmann equation in which case we still have the fractional diffusion but I want to emphasize that here the major difficulty is we do not have a strong convergence or we do not have a Hilbert expansion based convergence. So we know that Hilbert expansion actually gives us a lot of guidance when we design numerical methods because it gives us an idea of what's the magnitude of each term and how should we group the terms. But as you can see that when we derive in the derivation of fractional diffusion limit here we use this tricky choice of test function which gives us limited information on how we can construct the numerical methods. Alright so let me first show you how we compute this fractional laplacian and as I mentioned this is adopted from this work. So it goes as the following steps. So first the first step is we do a change of variable so this V is in the whole domain and if we do the change of variable not only a change of variable in the following way then c is in between minus 1 and 1 and this Lv you can consider that's a scaling parameter. So basically the larger the Lv you can see that the points are more spread out so usually if you want to compute a slow decay function then you want to choose this Lv to be large. But of course if you choose Lv to be large then there will be less points in the body so you have to strike a balance there. And then once we are in the finite domain from minus 1 to 1 then we can consider the Chebyshev polynomials and then if we do further change of variable if we let this arc cos i can see to be q then the Chebyshev polynomial just reduces to this trigonometric function. And now if we just work with the new variable q and keep in mind that the relation between q and original variable V is the following then we have to rewrite the fractional laplacian in this new variable q which looks a little bit complicated but the bottom line is it is an analytical expression. Now if we want to apply this fractional laplacian to a function so remember now the function if we write it into the new variable q then q is in the domain from 0 to pi so we can do an even extension and extend the function to the domain from 0 to 2 pi and in this way if we were to expand this f tilde q in the basis which is the cosine kq and because this is an even function so that expansion turns out to be quite easy it's just the discrete Fourier transform. And then the next thing is that we have to apply the fractional laplacian to this f of q and then it reduces to compute the fractional laplacian of the basis function and in this case the basis function is just this Fourier basis. But the thing is this fractional laplacian is no longer the original one so you cannot simply say that this one equals the modulus of k to the power of 2s times the Fourier basis instead we have to resort to the previous new formulation of this fractional laplacian. So actually we have so they have derived this analytical formula for this fractional laplacian applying on the Fourier basis but here I'm not going to list it in these slides it will be tedious but it can be pre-computed. Alright so let me show you the performance of this method by applying it to two kinds of f so the left hand side is this exponential decay function and the right hand side is this power law decay function. So for these two choices we have the analytical formula where we apply the fractional laplacian and this is the error between the numerical result and the analytical formula. So basically you can tell that the slow decay function is much harder to compute because with the same number of modes on the accuracy for this one is much lower than this one so you can see that if you want to reach the same accuracy then you will need a lot more modes for the slow decay function which is not surprising. Okay so if now we go back to the Levy-Folker Planck operator we have another term but this term is easy so we just rewrite it in terms of the new variable Q. Okay so now that we know how to compute the Levy-Folker Planck operator we have to deal with the stiffness and our idea is still we want to use this micro-micro decomposition but the key idea is that instead of saying that the macro part is rho times m we allow the macro part to depend on v as well but not in an arbitrary way. So although we have this dependence on v so the dependence is very mild so we require this x and v actually has the following relation. So this is inspired by this weak formulation. So once we have this special form of decomposition and we plug it into the original equation and this is what we get and again due to the special form of this eta we see that these two red term cancelled and the blue term changes to the fractional aplasia in x and this term goes to zero because m is the equilibrium and this i is just you can just consider as the residual where you try to apply sort of product rule where you apply this fractional aplasia. Okay so after this simplification the equation becomes the following so this equation is completely equivalent to the original one but then we do a splitting. So we split it into two equations one is for eta for the macro part and one is for the micro part. So the rationale behind this splitting is that we just keep in mind that we want to what we want is the in the end is the fractional aplasia and you can see that this eta the eta term in the original equation they are of the same order and they almost gives us what the diffusion equation should look like and then the rest of the thing is put into the micro part. So you can see that evolution is that eta evolves by itself but the row evolves according to eta and the initial condition is decomposed in the following way so eta is of course related to row but in the in the special form and and then the rest of the thing is given to g. So there is a slight thing we should be careful is that the whole derivation here relies on the fact that this eta is of this special form so we but we know that the but the evolution that eta obeys is the following equation so we need to make sure that eta that the the form of eta is preserved along the dynamics which is not hard to show. Alright so so we check a little bit about this splitting system so we have to make sure that this splitting is not too wild so we check the energy stability of this splitting and it turns out that the energy for the f so that which is they reconstructed or the total or just the original function f and also the energy for the macro and the micro part individually are uniformly bounded in time. Okay now the next thing is we need to solve this splitting system and if you look at the equation for g you will you will notice that it is as hard to solve as the original equation for f so our idea is we do an operator splitting so basically we split the the convection part from the collision part so in the first stage of the splitting we just focus on the collision part and then in the second stage we focus on the convection part but we also do a auxiliary term so this is like subtract add and subtract and so we add we subtract this gamma g in the first splitting and add it back in the second splitting. So for this splitting method compared to the original one the the advantage is quite obvious so first of all because of this auxiliary term it alleviates the eoconditioning so because we know that this L-S has this non-empty non-space so if you were to invert the system when epsilon is very small then the system to be inverted will be very stiff so actually it can be shown from this condition number so you can see that when epsilon is very small the condition number just grow dramatically but with this additional term the condition number just grow mildly so that's one advantage and the second advantage is to reduce the computational cost so here instead of just because so if you consider the original equation and if you treat it them implicitly then you have to deal with x and y simultaneously so that means you have to invert a very large system but once we do this the operator splitting then in the first stage we just need to do to invert an equation only for v so so x we can treat it locally and the same for the second one we just treat v locally so the computational cost is largely reduced and the third one is for the asymptotic property concern which I will explain in the next slides so basically we have the following result that the numerical solution obtained by this reconstruction will satisfies the following equation when epsilon goes to zero and to see why this is true first let's from this reconstructed solution we can see that so the first one this eta because the eta satisfies the fractional diffusion so it can be replaced in the following way and then with a pretty cheap argument we can swap the average with this fractional equation so basically we can move this average inside and then the average or the integral of f is rho and so you can see that this is pretty much what we want with two additional terms depending on g so if you can show that g or the average of g is a small term then we are good so we resort to the equation for g so here the first from the first stage we just need to control the magnitude of g so as long as it is bounded and then from the second stage it is clearly that it will be a small term so it will be order epsilon to the power of 2s so from this derivation you can also see the advantage of doing this operator splitting is that it gives us a little freedom so for the first stage we don't require so we don't need g to be small just immediately so for the first stage it just needs to be bounded but then the second stage gives the power of driving g to a small magnitude so that's why in the end it converges to the equation for the fractional flash okay so let me show you some numerical results so the first one is the spatially homogeneous case and here we check the long time behavior so here when s equals 0.5 we have the analytical form of the equilibrium so here the left hand side is we compare the numerical solution with this analytical formula and here we can only plot the body part so you can see that the body part they match very well and so it is with the tail part and then the last figure that's for the exponential decay of the relative entropy and we also check two other cases where s equals 0.6 and 0.8 and in these two cases we do not have the analytical formula to compare with so we only see whether we obtain the correct tail behavior and as you can see that the tail is proportional to this power and the next is for the spatially inhomogeneous case and you can see that we have the uniform convergence for a different choices of epsilon and these two figures are for different s and this is to check the energy stability and again this is for different choices of epsilon so for the last one you can see that there is a sharp transition this is due to the presence of initial layer and also we compare the solution with the solution to the original kinetic equation so we use a IMAX scheme with a very fine mesh and when epsilon equals 1 this is the comparison for our solution with the solution obtained by this IMAX scheme and the last one is to show the comparison with the diffusion limit and this one is to show the asymptotic preserving property basically if you decrease your epsilon then the difference between your f and the local equilibrium should be of other epsilon. Alright so to conclude that we have developed this asymptotic preserving scheme and the key idea is the spectral method with no local basis that is to treat the V-focal Planck operator and to design this asymptotic preserving scheme we need to use this micro micro decomposition and also to solve the split system we need to use the operator splitting so with that I'm going to end my talk and thanks for your attention. Thanks for your nice work and talk so you mentioned the so Cess Brung and Menon and Trevisa and several years ago proved this fractional diffusion limit and if you fact about the two years or three years ago Cess Brung and Menon generalized their work to the half space with non-local boundary condition so my question is is there any possible and to I mean numerically capture this non-local boundary behavior with this APS scheme? That's my question. We haven't tried that yet but yes I believe that we I don't know if this method can directly apply but I believe that we can develop sort of this APS scheme for that kind for the equation with boundary yes. Yeah okay thank you so because this boundary condition is very interesting because the boundary condition is kind of a non-local Neumann boundary condition so yeah. Yeah yeah yeah so I know that when you have fractional aplasia and you have the boundary domain then yes the boundary condition is highly non-trivial. Okay yeah thanks. How is your method sensitive to the operator you choose? Which operator? Yeah which operator? Maybe you'll have a message on chat you cannot say it? If it's not a log... No I see. If it is not well so I have to say that the way we do this micromicro decomposition that is specifically designed for the passage from Levy-Folker Planck to the fractional diffusion. So if it is not this operator then the mechanism for driving to the fractional diffusion might be different than the design of this micromicro decomposition will need to change so in that case I would say that yes it is specific to the Levy-Folker Planck operator. He says that you can unite the working teams to create the work. Yeah, significant difference there. Can you do if you are the type of people of the ATC in the world like aurance buildings
equilibrium for the fractional Laplacian case, diffusive scaling. So with this anonymous And there are two main challenges. So one is the non-locality that comes from computing this Chebyshev polynomial based spectral method and so use this tricky choice of test function which we do the change of variable not only a change of next thing is that we have to apply the side is this power law decay function. So for then in the second stage we focus on the operator splitting then in the first stage we that we have the uniform convergence for a have developed this asymptotic preserving scheme we need to use the operator splitting so with decomposition that is specifically designed for