Transcription
The automatic solution of PDEsusing a global spectral method2014 CBMS-NSF Conference, 29th June 2014Alex TownsendPhD studentUniversity of Oxford(with Sheehan Olver)Supervised by Nick Trefethen, University of Oxford.Supported by EPSRC grant EP/P505666/1.
IntroductionChebfunPiecewise smooth(2010)Endpointsingularities(2010*)Blow upfunctions(2011)Reducing regularityChebfunLinear rdinary differential ainsLinear PDEsNonlinearPDEsTwo dimensionsAlex Townsend @ Oxford1/21
IntroductionChebop: Spectral collocation for ODEsIn 2008: Overload the MATLAB backslash command \ for operators [Driscoll,Bornemann, & Trefethen 2008].L chebop(@(x,u) diff(u,2)-x.*u,[-30 30]); % Airy equationL.lbc 1; L.rbc 0;% Set boundary conditionsu L \ 0; plot(u)% Solve and plot50 5 10 30Alex Townsend @ Oxford 20 1001020302/21
IntroductionSpectral collocation basicsGiven values on a grid, what are the values of the derivative on that same grid?:u’ u’u’u’u1 u2u5u9x1 x2x5x912 0 u1 u1 Dn . . , 0 unun59Dn diffmat(n).For example, u0 (x) cos(x)u(x) is represented asLn Dn diag (cos(x1 ), . . . , cos(xn )) Rn n .Alex Townsend @ Oxford3/21
IntroductionWhy do spectral methods get a bad press?1. Dense matrices.2. Ill-conditioned matrices.3. When has it converged? Tricky.See, for example: [Canuto et al. 07], , [Fornberg 98], [Trefethen 00].Alex Townsend @ Oxford4/21
IntroductionWhy do spectral methods get a bad press?1. Dense matrices.2. Ill-conditioned matrices.3. When has it converged? Tricky.See, for example: [Canuto et al. 07], , [Fornberg 98], [Trefethen 00].Alex Townsend @ Oxford4/21
IntroductionWhy do spectral methods get a bad press?1. Dense matrices.2. Ill-conditioned matrices.3. When has it converged? Tricky.Condition number1510101051001005001000See, for example: [Canuto et al. 07], , [Fornberg 98], [Trefethen 00].Alex Townsend @ Oxford4/21
IntroductionWhy do spectral methods get a bad press?1. Dense matrices.2. Ill-conditioned matrices.3. When has it converged? Tricky.Condition numberError in solution15100101010 510510 1010010050010002004006008001000See, for example: [Canuto et al. 07], , [Fornberg 98], [Trefethen 00].Alex Townsend @ Oxford4/21
A fast and well-conditioned spectral methodDifferentiation operatorWork with coefficients: Spectral methods do not have to result in dense, illconditioned matrices. (Just don’t discretize the differentiation operator faithfully.)The idea is to use simple relations between Chebyshev polynomials: 1 (Uk Uk 2 ) , k 2, 2 dTk 1 kUk 1 , k 1,Tk k 1, 2 U1 , 0,dxk 0, U0 ,k 0. 0 1 2 ,D 3 . . . 1 1 0 2 1 0 1 22 .S 11 0 2 2 . . . . . . . . . Olver & T., A fast and well-conditioned spectral method, SIAM Review, 2013.Alex Townsend @ Oxford5/21
A fast and well-conditioned spectral methodMultiplication operatorTj Tk 1T j k 2 2a0 a1 a2 a3 a 1 2a0 a1 a21 M[a] a2 a1 2a0 a12 a3 a2 a1 2a0 . . . {zToeplitz 0. . . . . . a1 . . . 1 a2. . . 2 a3 . . . .} 1Tj k20 0 . . . . a2 a3 a4 . . . .a3 a4 a5 . . .a4 a5 a6 . . . . . .{z}0Hankel rank-1Multiplication is not a dense operator in finite precision. It is m-banded: mXXa(x) ak Tk (x) ãk Tk (x) O( ),k 0Alex Townsend @ Oxfordk 06/21
A fast and well-conditioned spectral methodWhat about this new spectral method?1. Almost banded matrices.2. Well-conditioned matrices.3. When has it converged? Trivial.Other approaches: [Clenshaw 57], [Greengard 91], [Shen 03].Alex Townsend @ Oxford7/21
A fast and well-conditioned spectral methodWhat about this new spectral method?1. Almost banded matrices.2. Well-conditioned matrices.3. When has it converged? Trivial.Other approaches: [Clenshaw 57], [Greengard 91], [Shen 03].Alex Townsend @ Oxford7/21
A fast and well-conditioned spectral methodWhat about this new spectral method?1. Almost banded matrices.2. Well-conditioned matrices.3. When has it converged? Trivial.Condition number21011001005001000Other approaches: [Clenshaw 57], [Greengard 91], [Shen 03].Alex Townsend @ Oxford7/21
A fast and well-conditioned spectral methodWhat about this new spectral method?1. Almost banded matrices.2. Well-conditioned matrices.3. When has it converged? Trivial.Condition numberError in solution210010110 1010010050010002004006008001000Other approaches: [Clenshaw 57], [Greengard 91], [Shen 03].Alex Townsend @ Oxford7/21
A fast and well-conditioned spectral methodFirst exampleu0 (x) x 3 u(x) 100 sin(20,000x 2 ),u( 1) 0.The exact solution is!Z xt4x4100e 4 sin(20,000t 2 )dt .u(x) e 4 11.2N chebop(@(x,u) · · · );N.lbc 0; u N \ f;Adaptively selects thediscretisation size.Forms a chebfun object[Chebfun V4.2].kũ uk 1.5 10 15 .1u(x)0.8degree(u) 203910.60.4time 15.5s0.20 0.2 1Alex Townsend @ Oxford 0.50x0.518/21
A fast and well-conditioned spectral methodAnother example1u(x) 0,u( 1) 1.1 50,000x 2The exact solution with a 50,000 is !tan 1 ( ax) tan 1 ( a)u(x) exp . au0 (x) 5101.0051 10Absolute erroru(x)degree(u) 50930.99510 15100.990.985 1Alex Townsend @ Oxford 20 0.50x0.5110 1Old chebopPreconditioned collocationUS method 0.500.51x9/21
A fast and well-conditioned spectral methodA high-order exampleu(10) (x) cosh(x)u(8) (x) cos(x)u(2) (x) x 2 u(x) 0u( 1) 0, u0 ( 1) 1, u(k ) ( 1) 0, k 2, 3, 4.00.5100.4 50.310 un 1 un 20.2u(x)0.10 0.1 10Z10 1510(ũ(x) ũ( x))2 1 0.2 1.3 10 14 . 20 0.3! 12110 0.4 0.5 1 25 0.50xAlex Townsend @ Oxford0.511020406080100120140n10/21
Chebop and Chebop2Convenience for the userL chebop(@(x,u) diff(u,2)-x.*u,[-30 30]); % Airy equationL.lbc 1; L.rbc 0;% Set boundary conditionsu L \ 0;% u chebfunConverthandle intodiscretisationinstructionsImposebcs andsolveà x bConstructdiscretisationA Rn nConvergedto thesolution?yesConstructa chebfunno, increase nL chebop2(@(x,y,u) laplacian(u) (1000 y)*u);% Helmholtz with gravityL.lbc 1; L.rbc 1; L.ubc 1; L.dbc 1;% Set boundary conditionsu L \ 0;% u chebfun2Converthandle intodiscretisationinstructionsConstruct amatrix equationwith an ny nxsolution matrixImposebcs andsolvematrixequationConvergedto thesolution?yesConstructa chebfun2no, increase nx or ny or bothAlex Townsend @ Oxford11/21
Chebop and Chebop2Convenience for the userL chebop(@(x,u) diff(u,2)-x.*u,[-30 30]); % Airy equationL.lbc 1; L.rbc 0;% Set boundary conditionsu L \ 0;% u chebfunConverthandle intodiscretisationinstructionsImposebcs andsolveà x bConstructdiscretisationA Rn nConvergedto thesolution?yesConstructa chebfunno, increase nL chebop2(@(x,y,u) laplacian(u) (1000 y)*u);% Helmholtz with gravityL.lbc 1; L.rbc 1; L.ubc 1; L.dbc 1;% Set boundary conditionsu L \ 0;% u chebfun2Converthandle intodiscretisationinstructionsConstruct amatrix equationwith an ny nxsolution matrixImposebcs andsolvematrixequationConvergedto thesolution?yesConstructa chebfun2no, increase nx or ny or bothAlex Townsend @ Oxford11/21
Interpreting user-defined inputAutomatic differentiationImplemented byforward-mode operatoroverloadingInterpret anonymousfunction as a sequenceof elementaryoperationsCan also calculateFréchet derivativesKey people:Ásgeir Birkisson andToby DriscollAlex Townsend @ Oxforduxx uyy 50u yu .* uxxu*uyyyuuu12/21
Low rank approximationNumerical rankFor A Cm n , SVD gives best rank k wrt 2-norm [Eckart & Young 1936]A min(m,n)Xσj uj vj j 1kXσj uj vj ,σk 1 tol.j 1For Lipschitz smooth bivariate functions [Schmidt 1909, Smithies 1937]f (x, y) Xσj uj (y)vj (x) j 1kXσj uj (y)vj (x).j 1For compact linear operators acting on functions of two variables,!L Xj 1Alex Townsend @ Oxfordyσj Lj Lxj kXyσj Lj Lxj .j 113/21
Low rank approximationNumerical rankFor A Cm n , SVD gives best rank k wrt 2-norm [Eckart & Young 1936]A min(m,n)Xσj uj vj j 1kXσj uj vj ,σk 1 tol.j 1For Lipschitz smooth bivariate functions [Schmidt 1909, Smithies 1937]f (x, y) Xσj uj (y)vj (x) j 1kXσj uj (y)vj (x).j 1For compact linear operators acting on functions of two variables,!L Xj 1Alex Townsend @ Oxfordyσj Lj Lxj kXyσj Lj Lxj .j 113/21
Low rank approximationNumerical rankFor A Cm n , SVD gives best rank k wrt 2-norm [Eckart & Young 1936]A min(m,n)Xσj uj vj j 1kXσj uj vj ,σk 1 tol.j 1For Lipschitz smooth bivariate functions [Schmidt 1909, Smithies 1937]f (x, y) Xσj uj (y)vj (x) j 1kXσj uj (y)vj (x).j 1For compact linear operators acting on functions of two variables,!L Xj 1Alex Townsend @ Oxfordyσj Lj Lxj kXyσj Lj Lxj .j 113/21
Low rank approximationDo the low rank stuff before discretizationLow rank-then-discretize: Instead of low rank techniques after discretization,do them before.For example, Helmholtz is of rank 2 2 u K 2 u (uxx K2K2K2K2u) (uyy u) (D2 I) I I (D2 I).2222Let A be your favourite ODE discretization of D2 K22 I,then (typically)AXI IXA T .In general, if L is of rank k we havekXAj XBjT Fj 1Alex Townsend @ Oxford14/21
Low rank approximationDo the low rank stuff before discretizationLow rank-then-discretize: Instead of low rank techniques after discretization,do them before.For example, Helmholtz is of rank 2 2 u K 2 u (uxx K2K2K2K2I) I I (D2 I).u) (uyy u) (D2 2222Let A be your favourite ODE discretization of D2 K22 I,then (typically)AXI IXA T .In general, if L is of rank k we havekXAj XBjT Fj 1Alex Townsend @ Oxford14/21
Low rank approximationDo the low rank stuff before discretizationLow rank-then-discretize: Instead of low rank techniques after discretization,do them before.For example, Helmholtz is of rank 2 2 u K 2 u (uxx K2K2K2K2I) I I (D2 I).u) (uyy u) (D2 2222Let A be your favourite ODE discretization of D2 K22 I,then (typically)AXI IXA T .In general, if L is of rank k we havekXAj XBjT Fj 1Alex Townsend @ Oxford14/21
Low rank approximationComputing the rank of a partial differential operatorRecast differential operators as polynomials: Once you have polynomialscomputing the rank is easy.The rank ofL Ny XNxXaij (x, y)i 0 j 0 i j y i x jequals a TT-rank [Oseledets 2011] (between {x, s} and {y, t}) ofh(x, s, y, t) Ny XNxXi jaij (s, t)y x i 0 j 0Rank 1:ODEsTrivial PDEsAlex Townsend @ OxfordRank 2:Laplace, HelmholtzTransport, Heat, WaveBlack-ScholeskXcj (t, y)rj (s, x).j 1Rank 3:BiharmonicLots here.15/21
Low rank approximationComputing the rank of a partial differential operatorRecast differential operators as polynomials: Once you have polynomialscomputing the rank is easy.The rank ofL Ny XNxXaij (x, y)i 0 j 0 i j y i x jequals a TT-rank [Oseledets 2011] (between {x, s} and {y, t}) ofh(x, s, y, t) Ny XNxXi jaij (s, t)y x i 0 j 0Rank 1:ODEsTrivial PDEsAlex Townsend @ OxfordRank 2:Laplace, HelmholtzTransport, Heat, WaveBlack-ScholeskXcj (t, y)rj (s, x).j 1Rank 3:BiharmonicLots here.15/21
Low rank approximationComputing the rank of a partial differential operatorRecast differential operators as polynomials: Once you have polynomialscomputing the rank is easy.The rank ofL Ny XNxXaij (x, y)i 0 j 0 i j y i x jequals a TT-rank [Oseledets 2011] (between {x, s} and {y, t}) ofh(x, s, y, t) Ny XNxXi jaij (s, t)y x i 0 j 0Rank 1:ODEsTrivial PDEsAlex Townsend @ OxfordRank 2:Laplace, HelmholtzTransport, Heat, WaveBlack-ScholeskXcj (t, y)rj (s, x).j 1Rank 3:BiharmonicLots here.15/21
Low rank approximationConstruct a nx by ny generalised Sylvester matrix equationIf the PDE is Lu f , where L is of rank-k then we solve for X Cny nx in,kXσj Aj XBjT F,Aj Cny ny ,Bj Cnx nx .j 1X solution’s coefficients Aj Alex Townsend @ OxfordyAj , Bj 1D spectral discretization of Lj , Lxj 16/21
Low rank approximationMatrix equation solversRank 1: A1 XB1T F. Solve A1 Y F, then B1 X T Y T .Rank 2: A1 XB1T A2 XB2T F. Generalised Sylvester solver (RECSY)[Jonsson & Kågström, 2002].Rank k, k 3: Solve N N system using almost banded structure.2/2O(N3Execution time110)O(N2)10010N)(O 110blue rank 1green rank 2red rank 3 210 310110Alex Townsend @ Oxford2 N1031017/21
ExamplesHelmholtz equation 2 u 2ω2 u 0,u(x, 1) f (x, 1),u( 1, y) f ( 1, y),where f (x, y) cos(ωx) cos(ωy).Wave numbers vs. polynomial degreeω 5060015000.50ny400300200 0.5100 1 1Alex Townsend @ Oxford 0.50x0.5100100200300400500ω18/21
ExamplesVariable helmholtz equationN chebop2(@(x,y,u) laplacian(u) 10000(1/2 sin(x)ˆ2).*cos(y)ˆ2.*u);N.lbc 1; N.rbc 1; N.ubc 1; N.dbc 1;u N \ chebfun2(@(x,y) cos(x.*y));0Absolute maximum error10 510 1010 1510200N 1,050,625,Alex Townsend @ Oxforderror 1.47 10 13 ,400 600N8001000time 44.2s.19/21
ExamplesWave and Klein–Gordon equationN chebop2(@(u) diff(u,2,1) - diff(u,2,2) 5*u); % u tt - u xx 5uN.dbc @(x,u) [u-exp(-10*x) diff(u)]; N.lbc 0; N.rbc 0;u N \ 0;Alex Townsend @ Oxford20/21
ConclusionSpectral methods do not have to be ill-conditioned. (Don’t discretizedifferentiation faithfully.)Spectral methods are extremely convenient and flexible.As of 2014, global spectral methods are heavily restricted to a fewgeometries.Thank you for listeningAlex Townsend @ Oxford21/21
University of Oxford (with Sheehan Olver) Supervised by Nick Trefethen, University of Oxford. Supported by EPSRC grant EP/P505666/1. Introduction Chebfun Chebfun . Preconditioned collocation US method Alex Townsend @ Oxford 9/21. A fast and well-conditioned spectral method A high-order example u(10) .
