-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy path2020-05-19.html
executable file
·243 lines (189 loc) · 15.9 KB
/
2020-05-19.html
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
<!DOCTYPE html>
<html lang="en">
<head>
<meta charset="UTF-8">
<meta http-equiv="X-UA-Compatible" content="IE=Edge">
<meta name="description" content="">
<meta name="keywords" content="">
<meta name="author" content="">
<meta name="viewport" content="width=device-width, initial-scale=1, maximum-scale=1">
<title>Exponential integrators for stiff PDEs</title>
<!--
Template 2085 Neuron
http://www.tooplate.com/view/2085-neuron
-->
<link rel="stylesheet" href="css/bootstrap.min.css">
<link rel="stylesheet" href="css/all.css">
<link rel="stylesheet" href="css/magnific-popup.css">
<!-- Main css -->
<link rel="stylesheet" href="css/style.css">
<link href="https://fonts.googleapis.com/css?family=Lora|Merriweather:300,400" rel="stylesheet">
<!-- Global site tag (gtag.js) - Google Analytics -->
<script async src="https://www.googletagmanager.com/gtag/js?id=UA-108246943-1"></script>
<script>
window.dataLayer = window.dataLayer || [];
function gtag(){dataLayer.push(arguments);}
gtag('js', new Date());
gtag('config', 'UA-108246943-1');
</script>
<!-- Google AdSense -->
<script async src="https://pagead2.googlesyndication.com/pagead/js/adsbygoogle.js?client=ca-pub-8299683413193339" crossorigin="anonymous"></script>
<!--<script async src="https://fundingchoicesmessages.google.com/i/pub-8299683413193339?ers=1" nonce="AuY4vKpDElA9N7YwmFfMjg"></script>
<script nonce="AuY4vKpDElA9N7YwmFfMjg">(function() {function signalGooglefcPresent() {if (!window.frames['googlefcPresent']) {if (document.body) {const iframe = document.createElement('iframe'); iframe.style = 'width: 0; height: 0; border: none; z-index: -1000; left: -1000px; top: -1000px;'; iframe.style.display = 'none'; iframe.name = 'googlefcPresent'; document.body.appendChild(iframe);} else {setTimeout(signalGooglefcPresent, 0);}}}signalGooglefcPresent();})();</script>-->
</head>
<body>
<!-- PRE LOADER -->
<div class="preloader">
<div class="sk-spinner sk-spinner-wordpress">
<span class="sk-inner-circle"></span>
</div>
</div>
<!-- Navigation section -->
<div class="navbar navbar-default navbar-static-top" role="navigation">
<div class="container">
<div class="navbar-header">
<button class="navbar-toggle" data-toggle="collapse" data-target=".navbar-collapse">
<span class="icon icon-bar"></span>
<span class="icon icon-bar"></span>
<span class="icon icon-bar"></span>
</button>
<a href="index.html" class="navbar-brand">H. Montanelli</a>
</div>
<div class="collapse navbar-collapse">
<ul class="nav navbar-nav navbar-right">
<li><a href="about.html">About</a></li>
<li><a href="blog.html">Blog</a></li>
<li><a href="publications.html">Publications</a></li>
<li><a href="research.html">Research</a></li>
<li><a href="students.html">Students</a></li>
<li><a href="talks.html">Talks</a></li>
<li><a href="teaching.html">Teaching</a></li>
<li><a href="https://www.paypal.com/donate/?hosted_button_id=UCLCSJLFL433E"><b>Donate</b></a></li>
</ul>
</div>
</div>
</div>
<!-- Home Section -->
<section id="home" class="main-about parallax-section">
<div class="overlay"></div>
<div class="container">
<div class="row">
<div class="col-md-12 col-sm-12">
<h1>Exponential integrators for stiff PDEs</h1>
</div>
</div>
</div>
</section>
<!-- Blog Single Post Section -->
<section id="about">
<div class="container">
<div class="row">
<div class="col-md-offset-1 col-md-10 col-sm-12">
<p><i>May 19, 2019 — Support my next blog post, <a href="https://www.paypal.com/donate/?hosted_button_id=UCLCSJLFL433E">buy me a coffee</a> ☕.</i></p>
<center>
<img src="/blog/iconic.png" class="img-responsive">
</center>
<p><a href="https://arxiv.org/pdf/1604.08900.pdf">One of the papers</a> I wrote during my Ph.D. at Oxford with <a href="https://www.numerical.rl.ac.uk/people/nbootland/">Niall Bootland</a> was just accepted in <a href="https://www.journals.elsevier.com/mathematics-and-computers-in-simulation">Mathematics and Computers in Simulation</a>. In this post, I review the main ideas of the paper.</p>
<h2>Periodic stiff PDEs</h2>
<p>We are interested in computing smooth solutions in a periodic domain in 1D/2D/3D of stiff PDEs of the form
$$
u_t = \mathcal{L}u + \mathcal{N}(u), \quad u(0,\boldsymbol{x})=u_0(\boldsymbol{x}),
$$
where \(u(t,\boldsymbol{x})\) is a function of time \(t\) and space \(\boldsymbol{x}\), \(\mathcal{L}\) is a linear differential operator with constant coefficients, and \(\mathcal{N}\) is a nonlinear operator of lower order with constant coefficients.</p>
<p>In applications, PDEs of this kind typically arise when two or more different physical processes are combined, and many PDEs of interest in
science and engineering take this form. For example, the <a href="https://en.wikipedia.org/wiki/Korteweg%E2%80%93de_Vries_equation">Korteweg–de Vries equation</a> \(u_t = -u_{xxx} - uu_x\), the starting point of the study of nonlinear waves and solitons, couples third-order linear dispersion with first-order convection, and the <a href="http://en.wikipedia.org/wiki/Allen–Cahn_equation">Allen–Cahn equation</a> \(u_t = \epsilon u_{xx} + u - u^3\) couples second-order linear
diffusion with a nondifferentiated cubic reaction term.</p>
<p>Often a system of equations rather than a single scalar equation is involved, for example in <a href="https://en.wikipedia.org/wiki/Reaction%E2%80%93diffusion_system">reaction-diffusion systems</a> such as the Gray–Scott and Schnakenberg equations, which involve two components coupled together. (The importance of coupling of nonequal diffusion constants in science was made famous by <a href="https://en.wikipedia.org/wiki/Alan_Turing">Alan Turing</a> in 1952 in the most highly-cited of all his papers.) Fourth-order terms also arise, for example in the <a href="https://en.wikipedia.org/wiki/Cahn%E2%80%93Hilliard_equation">Cahn–Hilliard equation</a>, whose solutions describe structures of alloys, and in the
<a href="https://en.wikipedia.org/wiki/Kuramoto%E2%80%93Sivashinsky_equation">Kuramoto–Sivashinsky equation</a>, related to combustion problems among others, whose solutions are chaotic. The figure at the top of this post shows six examples of solutions of such PDEs.</p>
<h2>Discretization in Fourier space</h2>
<p>Our paper describes and compares specialized methods that take advantage of two special features of the family of PDEs above. The first one is the periodic boundary conditions.
This allows us to discretize the spatial component with a <a href="https://en.wikipedia.org/wiki/Spectral_method">Fourier spectral method</a> on \(N\) points; the PDE becomes a system of \(N\) ODEs,
$$
\widehat{u}' = \mathbf{L}\widehat{u} + \mathbf{N}(\widehat{u}), \quad \widehat{u}(0)=\widehat{u}_0,
$$
where \(\widehat{u}(t)\) is the vector of \(N\) Fourier coefficients of the <a href="https://en.wikipedia.org/wiki/Trigonometric_interpolation">trigonometric interpolant</a> of \(u(t,X)\) at time \(t\), and \(\mathbf{L}\) (a \(N\times N\) matrix) and \(\mathbf{N}\) are the discretized versions of \(\mathcal{L}\) and \(\mathcal{N}\) in Fourier space.</p>
<p>For example, in 1D on \([0, 2\pi]\) with \(\mathcal{L}u=u_{xx}\) and an even number \(N\) of equispaced grid points \(\{x_j=2\pi j/N\}_{j=0}^{N-1}\), we look for a solution \(u(t,x)\) of the form
$$
u(t,x) \approx \sum_{k=-N/2}^{N/2}{\hspace{-0.3cm}}'{\;\,} \widehat{u}_k(t) e^{ikx}
$$
with <a href="https://en.wikipedia.org/wiki/Fourier_series">Fourier coefficients</a>
$$
\widehat{u}_k(t) = \frac{1}{N}\sum_{j=0}^{N-1}u(t,x_j)e^{-ikx_j},\;\; -\frac{N}{2}\leq k\leq \frac{N}{2}-1,
$$
and with \(\widehat{u}_{N/2}(t) = \widehat{u}_{-N/2}(t)\).
(The prime on the summation sign signifies that the terms \(k=\pm N/2\) are halved.) Since <a href="https://en.wikipedia.org/wiki/Fast_Fourier_transform">FFT</a> codes only store \(N\) coefficients, the vector \(\widehat{u}(t)\) is defined as
$$
\widehat{u}(t) = \Big(\frac{\widehat{u}_{-N/2}}{2}+\frac{\widehat{u}_{N/2}}{2}, \widehat{u}_{-N/2+1}(t),\ldots,\widehat{u}_{N/2-1}(t)\Big)^T.
$$
For this PDE, \(\mathbf{L}=\mathbf{D}_N^{(2)}\) is the (diagonal) <i>second-order Fourier differentiation matrix</i> with entries \(-k^2\), \(-N/2\leq k\leq N/2-1\). Note that stiffness is related to \(\mathbf{L}\) having large <a href="https://en.wikipedia.org/wiki/Eigenvalues_and_eigenvectors">eigenvalues</a> since the stability of spectral methods for time-dependent PDEs requires that the eigenvalues, scaled by the time-step, lie in the stability region of the time-stepping formula.</p>
<h2>A simple exponential integrator</h2>
<p>The second special feature of these PDEs is that they are is semilinear, i.e., the higher-order terms of the equation are linear. <a href="https://en.wikipedia.org/wiki/Exponential_integrator">Exponential integrators</a> are a class of numerical methods for systems of ODEs that are aimed at taking advantage of this. The linear part \(\mathbf{L}\), responsible for the stiffness, is integrated exactly using the <a href="https://en.wikipedia.org/wiki/Matrix_exponential">matrix exponential</a> while a numerical scheme is applied to \(\mathbf{N}\).</p>
<p>One of the simplest exponential integrators, commonly known as the Exponential Time Differencing (ETD) Euler method, is given by
$$
\widehat{u}^{n+1} = e^{h\mathbf{L}}\widehat{u}^n + h\varphi_1(h\mathbf{L})\mathbf{N}(\widehat{u}^n),
$$
where \(h=t_{n+1}-t_n\) is the time-step and
$$
\varphi_1(z) = \frac{e^z-1}{z}.
$$
Note that, in practice, the nonlinear evaluations \(\mathbf{N}(\widehat{u}^n)\) is carried out in value space, that is, \(\mathbf{N}(\widehat{u}^n)\) means \(\mathbf{F}(\mathbf{N}(\mathbf{F}^{-1}\widehat{u}^n))\), with discrete Fourier transform \(\mathbf{F}\). Since the matrix \(\mathbf{L}\) diagonal, computing the matrix exponential and matrix-vector products cost only \(\mathcal{O}(N)\) operations. The dominant cost per time-step is therefore the cost of an FFT, i.e., \(\mathcal{O}(N\log N)\) operations.</p>
<h2>Results of our numerical experiments and Chebfun</h2>
<p>We compare in our paper 30 exponential integrators of fourth and higher order on 11 model problems in 1D, 2D, and 3D. Our conclusion is that it is hard to do much better than one of the simplest of these formulas, the ETDRK4 scheme of Cox and Matthews.</p>
<p>Our numerical experiments were performed using <a href="https://www.mathworks.com/products/matlab.html">MATLAB</a> and have been embedded within <a href="https://www.chebfun.org/">Chebfun</a>. More specifically, the <code>spin</code>, <code>spin2</code> and <code>spin3</code> codes implement a Fourier spectral method and exponential integrators to solve PDEs in 1D, 2D, and 3D periodic domains. The simplest way to see <code>spin</code> in action is to type simply <code>spin('ks')</code> (for the Kuramoto–Sivashinsky equation) or <code>spin2('gl')</code> (for the 2D Ginzburg–Landau equation) to invoke an example computation.
It is also possible to define your own PDE using the <code>spinop</code> class. Check out the <a href="https://www.chebfun.org/docs/guide/guide19.html">documentation</a> for details.</p>
<hr>
<h4>Blog posts about spectral methods</h4>
<p>2020 <a href="2020-05-19.html">Exponential integrators for stiff PDEs</a></p>
<p>2018 <a href="2018-12-05.html">Computer-assisted proofs for PDEs</a></p>
<p>2018 <a href="2018-02-27.html">Spherical caps in cell polarization</a></p>
<p>2018 <a href="2018-01-25.html">Solving nonlocal equations on the sphere</a></p>
<p>2018 <a href="2018-01-04.html">Gibbs phenomenon and Cesàro mean</a></p>
<p>2017 <a href="2017-10-26.html">Solving PDEs on the sphere</a></p>
<p>2017 <a href="2017-10-09.html">When planets dance</a></p>
</div>
</div>
</div>
</section>
<!-- Footer Section -->
<footer>
<div class="container">
<div class="row">
<div class="col-md-4 col-md-offset-1 col-sm-6">
<h3>Contact</h3>
<p style="color:white"><i class="fa fa-send-o"></i>[email protected]</p><br>
<a href="https://www.inria.fr/en"><img src="/images/inria.png" class="img-responsive"></a>
</div>
<div class="col-md-4 col-md-offset-1 col-sm-6">
<h3 style="opacity:0;">Contact</h3>
<p style="color:white">Copyright © 2024 Hadrien Montanelli</p><br>
<a href="https://www.ip-paris.fr/en"><img src="/images/polytechnique.png" class="img-responsive"></a>
</div>
<div class="clearfix col-md-12 col-sm-12">
<hr style="color:white;">
</div>
<div class="col-md-12 col-sm-12">
<ul class="social-icon">
<li><a href="https://github.com/Hadrien-Montanelli" class="fa-brands fa-github" style="font-size:30px;color:white"></a></li>
<li><a href="https://scholar.google.com/citations?user=Bjmkfe8AAAAJ&hl=en&oi=sra/" class="fa-brands fa-google" style="font-size:30px;color:white"></a></li>
<li><a href="https://www.linkedin.com/in/hadrien-montanelli/" class="fa-brands fa-linkedin" style="font-size:30px;color:white"></a></li>
<li><a href="https://orcid.org/0009-0005-1742-9828" class="fa-brands fa-orcid" style="font-size:30px;color:white"></a></li>
<li><a href="https://twitter.com/drmontanelli" class="fa-brands fa-x-twitter" style="font-size:30px;color:white;"></a></li>
</ul>
</div>
</div>
</div>
</footer>
<!-- Back top -->
<a href="#back-top" class="go-top"><i class="fa fa-angle-up"></i></a>
<!-- SCRIPTS -->
<script src="js/jquery.js"></script>
<script src="js/bootstrap.min.js"></script>
<script src="js/particles.min.js"></script>
<script src="js/app.js"></script>
<script src="js/jquery.parallax.js"></script>
<script src="js/smoothscroll.js"></script>
<script src="js/custom.js"></script>
<script src="https://polyfill.io/v3/polyfill.min.js?features=es6"></script>
<script id="MathJax-script" async src="https://cdn.jsdelivr.net/npm/mathjax@3/es5/tex-mml-chtml.js"></script>
</body>
</html>