lecture2 dg

%23%20%2F%2F%2F%20script%0A%23%20requires-python%20%3D%20%22%3E%3D3.12%22%0A%23%20dependencies%20%3D%20%5B%0A%23%20%20%20%20%20%22marimo%22%2C%0A%23%20%20%20%20%20%22matplotlib%22%2C%0A%23%20%20%20%20%20%22mograder%22%2C%0A%23%20%20%20%20%20%22numpy%22%2C%0A%23%20%20%20%20%20%22scipy%22%2C%0A%23%20%5D%0A%23%20%2F%2F%2F%0A%0Aimport%20marimo%0A%0A__generated_with%20%3D%20%220.25.0%22%0Aapp%20%3D%20marimo.App(width%3D%22medium%22)%0A%0A%0A%40app.cell%0Adef%20_()%3A%0A%20%20%20%20import%20marimo%20as%20mo%0A%20%20%20%20import%20matplotlib.pyplot%20as%20plt%0A%20%20%20%20import%20numpy%20as%20np%0A%20%20%20%20from%20mograder.runtime%20import%20check%2C%20hint%0A%20%20%20%20from%20numpy.polynomial%20import%20legendre%0A%20%20%20%20from%20scipy.integrate%20import%20solve_ivp%0A%0A%20%20%20%20plt.rcParams.update(%0A%20%20%20%20%20%20%20%20%7B%0A%20%20%20%20%20%20%20%20%20%20%20%20%22figure.dpi%22%3A%20110%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%22axes.grid%22%3A%20True%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%22grid.alpha%22%3A%200.3%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%22font.size%22%3A%2011%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%22legend.fontsize%22%3A%209%2C%0A%20%20%20%20%20%20%20%20%7D%0A%20%20%20%20)%0A%0A%20%20%20%20def%20incomplete(*objs)%3A%0A%20%20%20%20%20%20%20%20%22%22%22True%20if%20any%20object%20is%20an%20unfinished-exercise%20placeholder%20(None%20or%20...).%22%22%22%0A%20%20%20%20%20%20%20%20for%20obj%20in%20objs%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20if%20obj%20is%20None%20or%20obj%20is%20Ellipsis%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20return%20True%0A%20%20%20%20%20%20%20%20%20%20%20%20if%20isinstance(obj%2C%20(tuple%2C%20list))%20and%20incomplete(*obj)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20return%20True%0A%20%20%20%20%20%20%20%20return%20False%0A%0A%20%20%20%20%23%20Period%20of%20the%20traveling%20wave%20sin(x%20-%20t)%20on%20the%20periodic%20domain%20%5B0%2C%202%CF%80)%0A%20%20%20%20T%20%3D%202%20*%20np.pi%0A%20%20%20%20return%20T%2C%20check%2C%20hint%2C%20incomplete%2C%20legendre%2C%20mo%2C%20np%2C%20plt%2C%20solve_ivp%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%20Lecture%202%20%E2%80%94%20Solving%20PDEs%3A%20(pseudo)spectral%20and%20nodal%20discontinuous%20Galerkin%20methods%0A%0A%20%20%20%20Last%20time%2C%204th-order%20finite%20differences%20gave%20errors%20%24%5Cpropto%20N%5E%7B-4%7D%24.%20To%20gain%20another%20factor%20of%2010%20in%20accuracy%20we%0A%20%20%20%20needed%20%2410%5E%7B1%2F4%7D%5Capprox%201.8%24%20times%20more%20points%20(and%2C%20because%20of%20the%20CFL%20condition%2C%20%24%5Capprox%203.2%24%20times%20more%20work%20in%201d).%0A%20%20%20%20Today%20we%20build%20a%20method%20whose%20error%20falls%20**exponentially**%20with%20the%20number%20of%20points%20for%20smooth%20solutions%3A%0A%20%20%20%20the%20**nodal%20discontinuous%20Galerkin%20(DG)%20method**.%20We%20evolve%20the%20same%20scalar-wave%20system%20with%20the%20same%20time%20integrator%2C%0A%20%20%20%20so%20the%20two%20methods%20can%20be%20compared%20directly.%0A%0A%20%20%20%20%7C%20Part%20%7C%20Topic%20%7C%20%E2%89%88%20time%20%7C%0A%20%20%20%20%7C---%7C---%7C---%7C%0A%20%20%20%20%7C%201%20%7C%20The%20idea%2C%20and%20a%20derivation%20of%20nodal%20DG%20in%201d%20(derivation%20steps%20D1%E2%80%93D3)%20%7C%2020%20min%20%7C%0A%20%20%20%20%7C%202%20%7C%20Polynomial%20interpolation%3A%20Legendre%20polynomials%20and%20Legendre%E2%80%93Gauss%E2%80%93Lobatto%20points%20(Exercise%204)%20%7C%2010%20min%20%7C%0A%20%20%20%20%7C%203%20%7C%20Differentiating%20the%20interpolant%3B%20why%20nodal%20methods%20make%20nonlinearities%20easy%20(Exercise%205)%20%7C%2010%20min%20%7C%0A%20%20%20%20%7C%204%20%7C%20Boundary%20corrections%20(numerical%20fluxes)%3A%20required%20properties%2C%20upwinding%20and%20dissipation%20%7C%2015%20min%20%7C%0A%20%20%20%20%7C%205%20%7C%20A%20multi-element%20DG%20solver%20for%20the%20scalar%20wave%20(Exercise%206)%20%7C%2015%20min%20%7C%0A%20%20%20%20%7C%206%20%7C%20Convergence%3A%20%24h%24-%20and%20%24p%24-refinement%20(Exercise%207)%20%7C%2010%20min%20%7C%0A%0A%20%20%20%20As%20before%3A%20complete%20the%20cells%20marked%20%60%23%20YOUR%20CODE%20HERE%60%2C%20and%20watch%20the%20coloured%20check%20boxes.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20---%0A%20%20%20%20%23%23%20Part%201%20%E2%80%94%20From%20finite%20differences%20to%20nodal%20DG%0A%0A%20%20%20%20**Finite%20differences%20approximate%20the%20*derivative*.**%20A%20stencil%20combines%20a%20few%20neighbouring%20point%20values%0A%20%20%20%20to%20estimate%20%24%5Cpartial_x%20u%24.%20The%20error%20is%20set%20by%20the%20truncation%20of%20a%20local%20Taylor%20series%2C%20so%20it%20falls%0A%20%20%20%20algebraically%2C%20%24%5Cpropto%5CDelta%20x%5Ep%24.%0A%0A%20%20%20%20**Spectral%20methods%20approximate%20the%20*solution*.**%20Instead%2C%20represent%20the%20solution%20itself%20by%20a%0A%20%20%20%20high-degree%20polynomial%2C%0A%0A%20%20%20%20%24%24%0A%20%20%20%20u(x)%20%5Capprox%20u_h(x)%20%3D%20%5Csum_%7Bj%3D0%7D%5E%7BN%7D%20u_j%5C%2C%5Cell_j(x)%2C%0A%20%20%20%20%24%24%0A%0A%20%20%20%20and%20then%20differentiate%20the%20approximation%20**exactly**%3A%20%24%5Cpartial_x%20u_h%20%3D%20%5Csum_j%20u_j%5C%2C%5Cell_j'(x)%24.%0A%20%20%20%20There%20is%20no%20separate%20%22derivative%20approximation%22.%20Once%20we%20have%20chosen%20how%20to%20approximate%20the%20solution%2C%0A%20%20%20%20its%20derivative%20follows%20exactly.%20For%20smooth%20(analytic)%20solutions%20the%20approximation%20error%2C%20and%20with%20it%20the%20derivative%20error%2C%0A%20%20%20%20decreases%20**exponentially**%20with%20%24N%24.%0A%0A%20%20%20%20**Discontinuous%20Galerkin**%20combines%20this%20with%20the%20geometric%20flexibility%20of%20finite%20volumes%3A%20split%20the%20domain%20into%0A%20%20%20%20%24K%24%20elements%2C%20use%20a%20polynomial%20of%20degree%20%24N%24%20in%20each%20one%2C%20and%20allow%20the%20solution%20to%20be%20**discontinuous**%20across%20element%0A%20%20%20%20boundaries.%20Neighbouring%20elements%20communicate%20only%20through%20a%20*numerical%20flux*%20at%20their%20shared%20faces%20(Part%204).%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.callout(%0A%20%20%20%20%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20**Three%20knobs%20control%20a%20DG%20discretization%3A**%0A%0A%20%20%20%201.%20the%20**size%20of%20each%20element**%20%24%5CDelta%20x%24%20(%24h%24-refinement)%2C%0A%20%20%20%202.%20the%20**number%20of%20points%20per%20element**%20%24N%2B1%24%2C%20i.e.%20the%20polynomial%20degree%20%24N%24%20(%24p%24-refinement)%2C%0A%20%20%20%203.%20the%20**location%20of%20the%20points**%20inside%20each%20element.%0A%0A%20%20%20%20Knobs%202%20and%203%20**together**%20are%20what%20give%20exponential%20convergence%3A%20increasing%20%24N%24%20only%20helps%20if%20the%20points%20are%20placed%0A%20%20%20%20well%20(Part%202%20shows%20that%20equally%20spaced%20points%20can%20make%20things%20*worse*).%20Knob%201%20on%20its%20own%20gives%20algebraic%0A%20%20%20%20convergence%20%24%5Cpropto%20%5CDelta%20x%5E%7BN%2B1%7D%24%2C%20just%20like%20a%20high-order%20finite-difference%20method.%0A%20%20%20%20%22%22%22)%2C%0A%20%20%20%20%20%20%20%20kind%3D%22info%22%2C%0A%20%20%20%20)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%23%20Derivation%20of%20nodal%20DG%20in%201d%0A%0A%20%20%20%20We%20consider%20a%20system%20in%20**flux-balance%20form**%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%5Cpartial_t%20u%20%2B%20%5Cpartial_x%20F(u)%20%3D%20S(u).%0A%20%20%20%20%24%24%0A%0A%20%20%20%20Our%20scalar-wave%20system%20is%20of%20this%20form%20with%0A%0A%20%20%20%20%24%24%0A%20%20%20%20u%20%3D%20%5Cbegin%7Bpmatrix%7D%5CPsi%5C%5C%20%5CPi%5C%5C%20%5CPhi%5Cend%7Bpmatrix%7D%2C%5Cqquad%0A%20%20%20%20F(u)%20%3D%20%5Cbegin%7Bpmatrix%7D0%5C%5C%20%5CPhi%5C%5C%20%5CPi-%5Cgamma_2%5CPsi%5Cend%7Bpmatrix%7D%2C%5Cqquad%0A%20%20%20%20S(u)%20%3D%20%5Cbegin%7Bpmatrix%7D-%5CPi%5C%5C%200%5C%5C%20-%5Cgamma_2%5CPhi%5Cend%7Bpmatrix%7D.%0A%20%20%20%20%24%24%0A%0A%20%20%20%20(Check%3A%20the%20%24%5CPhi%24%20row%20reads%20%24%5Cpartial_t%5CPhi%20%2B%20%5Cpartial_x(%5CPi%20-%20%5Cgamma_2%5CPsi)%20%3D%20-%5Cgamma_2%5CPhi%24.)%0A%0A%20%20%20%20**Step%201%20%E2%80%94%20Elements%20and%20the%20reference%20element.**%20Split%20%24%5B0%2C2%5Cpi)%24%20into%20%24K%24%20elements%0A%20%20%20%20%24%5COmega_k%3D%5Bx_k%2C%20x_k%2B%5CDelta%20x%5D%24%20and%20map%20each%20one%20to%20the%20reference%20element%20%24%5Cxi%5Cin%5B-1%2C1%5D%24%3A%0A%0A%20%20%20%20%24%24%0A%20%20%20%20x(%5Cxi)%20%3D%20x_k%20%2B%20%5Cfrac%7B1%2B%5Cxi%7D%7B2%7D%5C%2C%5CDelta%20x%2C%20%5Cqquad%20J%20%5Cequiv%20%5Cfrac%7Bdx%7D%7Bd%5Cxi%7D%20%3D%20%5Cfrac%7B%5CDelta%20x%7D%7B2%7D.%0A%20%20%20%20%24%24%0A%0A%20%20%20%20**Step%202%20%E2%80%94%20Nodal%20expansion.**%20Choose%20%24N%2B1%24%20nodes%20%24%5Cxi_0%3C%5Cdots%3C%5Cxi_N%24%20and%20the%20Lagrange%20polynomials%0A%20%20%20%20%24%5Cell_j(%5Cxi)%20%3D%20%5Cprod_%7Bm%5Cneq%20j%7D%5Cfrac%7B%5Cxi-%5Cxi_m%7D%7B%5Cxi_j-%5Cxi_m%7D%24%2C%20which%20satisfy%20%24%5Cell_j(%5Cxi_i)%20%3D%20%5Cdelta_%7Bij%7D%24.%20In%20each%20element%0A%0A%20%20%20%20%24%24%0A%20%20%20%20u_h(t%2C%5Cxi)%20%3D%20%5Csum_%7Bj%3D0%7D%5EN%20u_j(t)%5C%2C%5Cell_j(%5Cxi)%2C%5Cqquad%0A%20%20%20%20F_h%20%3D%20%5Csum_%7Bj%3D0%7D%5EN%20F%5Cbig(u_j%5Cbig)%5C%2C%5Cell_j(%5Cxi)%2C%5Cqquad%0A%20%20%20%20S_h%20%3D%20%5Csum_%7Bj%3D0%7D%5EN%20S%5Cbig(u_j%5Cbig)%5C%2C%5Cell_j(%5Cxi).%0A%20%20%20%20%24%24%0A%0A%20%20%20%20The%20unknowns%20%24u_j%24%20are%20simply%20the%20**values%20at%20the%20nodes**.%20Note%20that%20%24F%24%20and%20%24S%24%20are%20evaluated%20*pointwise*%20at%20the%0A%20%20%20%20nodes%2C%20however%20nonlinear%20they%20are.%20This%20is%20what%20makes%20nodal%20methods%20so%20convenient%20for%20nonlinear%20problems.%0A%0A%20%20%20%20**Step%203%20%E2%80%94%20Galerkin%20condition.**%20We%20cannot%20satisfy%20the%20PDE%20everywhere%20with%20a%20polynomial%2C%20so%20we%20require%20the%20residual%0A%20%20%20%20to%20be%20orthogonal%20to%20every%20basis%20function%3A%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%5Cint_%7B%5COmega_k%7D%5Cell_i(x)%5Cleft%5B%5Cpartial_t%20u_h%20%2B%20%5Cpartial_x%20F_h%20-%20S_h%5Cright%5Ddx%20%3D%200%2C%5Cqquad%20i%20%3D%200%2C%5Cdots%2CN.%0A%20%20%20%20%24%24%0A%0A%20%20%20%20**D1%20%E2%80%94%20your%20turn%3A%20integrate%20the%20flux%20term%20by%20parts**%20to%20move%20the%20derivative%20onto%20%24%5Cell_i%24%20(the%20*weak%20form*).%0A%20%20%20%20You%20will%20get%20a%20boundary%20term%20%24%5Cbig%5B%5Cell_i%5C%2Cn%20F%5Cbig%5D%24%20evaluated%20at%20the%20two%20faces%2C%20with%20outward%20normal%20%24n%3D%5Cmp1%24%20at%20the%0A%20%20%20%20left%2Fright%20face.%0A%0A%20%20%20%20**Step%204%20%E2%80%94%20Numerical%20flux.**%20At%20a%20face%20the%20solution%20is%20double-valued%3A%20%24u%5E%7B%5Crm%20int%7D%24%20from%20our%20element%20and%20%24u%5E%7B%5Crm%20ext%7D%24%0A%20%20%20%20from%20the%20neighbour.%20We%20replace%20%24nF%24%20in%20the%20boundary%20term%20by%20a%20**numerical%20flux**%20%24G(u%5E%7B%5Crm%20int%7D%2C%20u%5E%7B%5Crm%20ext%7D)%24.%0A%20%20%20%20This%20is%20the%20*only*%20place%20where%20elements%20talk%20to%20each%20other.%0A%0A%20%20%20%20**D2%20%E2%80%94%20your%20turn%3A%20integrate%20by%20parts%20once%20more**%20(undoing%20D1%20for%20the%20volume%20integral)%20to%20obtain%20the%20**strong%20form**.%0A%20%20%20%20The%20boundary%20term%20should%20now%20contain%20%24G%20-%20nF%5E%7B%5Crm%20int%7D%24.%0A%0A%20%20%20%20**Step%205%20%E2%80%94%20Quadrature.**%20Evaluate%20all%20integrals%20with%20the%20%24N%2B1%24-point%20**Legendre%E2%80%93Gauss%E2%80%93Lobatto%20(LGL)**%20quadrature%20rule%20on%0A%20%20%20%20the%20same%20nodes%2C%20%24%5Cint_%7B-1%7D%5E%7B1%7D%20f%5C%2Cd%5Cxi%5Capprox%5Csum_m%20w_m%20f(%5Cxi_m)%24%2C%20which%20is%20exact%20for%20polynomials%20of%20degree%20%24%5Cle%202N-1%24.%0A%0A%20%20%20%20**D3%20%E2%80%94%20your%20turn%3A%20show%20that%20the%20mass%20matrix**%20%24M_%7Bij%7D%20%3D%20%5Cint_%7B%5COmega_k%7D%5Cell_i%5C%2C%5Cell_j%5C%2Cdx%24%20**becomes%20diagonal**%2C%0A%20%20%20%20and%20find%20its%20entries.%20(Is%20the%20quadrature%20exact%20for%20%24%5Cell_i%5Cell_j%24%3F)%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(mo)%3A%0A%20%20%20%20response_text%20%3D%20%22**D1.**%20*(weak%20form)*%20...%20**D2.**%20*(strong%20form)*%20...%20**D3.**%20*(mass%20matrix)*%20...%20Write%20your%20derivation%20here%20by%20editing%20this%20cell%3B%20LaTeX%20such%20as%20%24M_%7Bij%7D%24%20works.%22%0A%20%20%20%20%23%23%23%20BEGIN%20SOLUTION%0A%20%20%20%20response_text%20%3D%20r%22%22%22%0A%20%20%20%20**D1%20(weak%20form).**%20%24%5Cint%20%5Cell_i%5C%2C%5Cpartial_x%20F_h%5C%2Cdx%20%3D%20%5Cbig%5B%5Cell_i%20F_h%5Cbig%5D_%7Bx_k%7D%5E%7Bx_k%2B%5CDelta%20x%7D%20-%20%5Cint%20F_h%5C%2C%5Cpartial_x%5Cell_i%5C%2Cdx%24%2C%20so%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%5Cint_%7B%5COmega_k%7D%5Cell_i%5C%2C(%5Cpartial_t%20u_h%20-%20S_h)%5C%2Cdx%20-%20%5Cint_%7B%5COmega_k%7DF_h%5C%2C%5Cpartial_x%5Cell_i%5C%2Cdx%20%2B%0A%20%20%20%20%5Csum_%7B%5Crm%20faces%7D%5Cell_i%5C%2C%20n%5C%2CF%20%3D%200%2C%0A%20%20%20%20%24%24%0A%0A%20%20%20%20where%20the%20face%20sum%20is%20%24%5Cell_i(x_k%2B%5CDelta%20x)F(x_k%2B%5CDelta%20x)%20-%20%5Cell_i(x_k)F(x_k)%24%2C%20i.e.%20%24n%3D%2B1%24%20on%20the%20right%20face%20and%20%24n%3D-1%24%20on%20the%20left.%0A%0A%20%20%20%20**D2%20(strong%20form).**%20Replace%20%24nF%5Cto%20G%24%20in%20the%20face%20term%2C%20then%20integrate%20%24-%5Cint%20F_h%5Cpartial_x%5Cell_i%24%20by%20parts%20again%3A%0A%20%20%20%20%24-%5Cint%20F_h%5Cpartial_x%5Cell_i%20%3D%20%5Cint%5Cell_i%5Cpartial_x%20F_h%20-%20%5Csum_%7B%5Crm%20faces%7D%5Cell_i%5C%2CnF%5E%7B%5Crm%20int%7D%24.%20Hence%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%5Cint_%7B%5COmega_k%7D%5Cell_i%5C%2C(%5Cpartial_t%20u_h%20%2B%20%5Cpartial_x%20F_h%20-%20S_h)%5C%2Cdx%20%2B%20%5Csum_%7B%5Crm%20faces%7D%5Cell_i%5C%2C%5Cbig(G%20-%20nF%5E%7B%5Crm%20int%7D%5Cbig)%20%3D%200.%0A%20%20%20%20%24%24%0A%0A%20%20%20%20**D3%20(mass%20matrix).**%20%24M_%7Bij%7D%20%3D%20J%5Cint_%7B-1%7D%5E1%5Cell_i%5Cell_j%5C%2Cd%5Cxi%20%5Capprox%20J%5Csum_m%20w_m%5C%2C%5Cell_i(%5Cxi_m)%5Cell_j(%5Cxi_m)%0A%20%20%20%20%3D%20J%5Csum_m%20w_m%5C%2C%5Cdelta_%7Bim%7D%5Cdelta_%7Bjm%7D%20%3D%20J%5C%2Cw_i%5C%2C%5Cdelta_%7Bij%7D%24.%0A%20%20%20%20The%20product%20%24%5Cell_i%5Cell_j%24%20has%20degree%20%242N%24%2C%20one%20more%20than%20LGL%20quadrature%20integrates%20exactly%2C%20so%20this%20is%20a%20(very%20accurate)%0A%20%20%20%20approximation%20called%20*mass%20lumping*.%20Its%20payoff%20is%20a%20diagonal%2C%20trivially%20invertible%20mass%20matrix.%0A%20%20%20%20%22%22%22%0A%20%20%20%20%23%23%23%20END%20SOLUTION%0A%20%20%20%20mo.md(response_text)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%23%20Result%3A%20the%20semi-discrete%20nodal%20DG%20scheme%0A%0A%20%20%20%20With%20LGL%20nodes%20the%20endpoints%20%24%5Cxi_0%3D-1%24%20and%20%24%5Cxi_N%3D%2B1%24%20are%20nodes%2C%20so%20%24%5Cell_i(-1)%20%3D%20%5Cdelta_%7Bi0%7D%24%20and%20%24%5Cell_i(%2B1)%3D%5Cdelta_%7BiN%7D%24.%0A%20%20%20%20Applying%20the%20quadrature%20to%20the%20strong%20form%20and%20dividing%20by%20the%20diagonal%20mass%20matrix%20%24M_%7Bii%7D%3DJ%20w_i%24%3A%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%5Cboxed%7B%0A%20%20%20%20%5Cfrac%7Bdu_i%7D%7Bdt%7D%20%3D%20-%5Cfrac%7B1%7D%7BJ%7D%5Csum_%7Bj%3D0%7D%5EN%20D_%7Bij%7D%5C%2CF(u_j)%20%2B%20S(u_i)%20-%0A%20%20%20%20%5Cfrac%7B%5Cdelta_%7Bi0%7D%7D%7BJ%5C%2Cw_0%7D%5Cbig(G%20-%20nF%5E%7B%5Crm%20int%7D%5Cbig)%5CBig%7C_%7B%5Crm%20left%7D%20-%0A%20%20%20%20%5Cfrac%7B%5Cdelta_%7BiN%7D%7D%7BJ%5C%2Cw_N%7D%5Cbig(G%20-%20nF%5E%7B%5Crm%20int%7D%5Cbig)%5CBig%7C_%7B%5Crm%20right%7D%0A%20%20%20%20%7D%0A%20%20%20%20%24%24%0A%0A%20%20%20%20where%20%24D_%7Bij%7D%20%3D%20%5Cell_j'(%5Cxi_i)%24%20is%20the%20**differentiation%20matrix**.%20The%20boundary%20correction%20%24G-nF%5E%7B%5Crm%20int%7D%24%20is%0A%20%20%20%20%22**lifted**%22%20into%20the%20volume%20with%20weight%20%241%2F(Jw)%24%3B%20with%20LGL%20nodes%20it%20only%20touches%20the%20two%20endpoint%20nodes.%20For%20a%0A%20%20%20%20continuous%20exact%20solution%20the%20correction%20vanishes%20and%20we%20are%20left%20with%20%24%5Cpartial_t%20u%20%2B%20%5Cpartial_x%20F%20%3D%20S%24%20at%20every%20node.%0A%0A%20%20%20%20The%20rest%20of%20this%20lecture%20builds%20each%20ingredient%3A%20good%20nodes%20(Part%202)%2C%20%24D_%7Bij%7D%24%20(Part%203)%2C%20%24G%24%20(Part%204)%2C%20and%20then%20assembles%0A%20%20%20%20the%20solver%20(Part%205).%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20---%0A%20%20%20%20%23%23%20Part%202%20%E2%80%94%20Polynomial%20interpolation%3A%20where%20to%20put%20the%20points%0A%0A%20%20%20%20**Legendre%20polynomials**%20%24P_n(%5Cxi)%24%20are%20orthogonal%20on%20%24%5B-1%2C1%5D%24%3A%20%24%5Cint_%7B-1%7D%5E1%20P_m%20P_n%5C%2Cd%5Cxi%20%3D%20%5Cfrac%7B2%7D%7B2n%2B1%7D%5Cdelta_%7Bmn%7D%24.%0A%20%20%20%20The%20**Legendre%E2%80%93Gauss%E2%80%93Lobatto%20(LGL)%20points**%20of%20degree%20%24N%24%20are%20the%20endpoints%20%24%5Cpm1%24%20together%20with%20the%20%24N-1%24%20roots%20of%20%24P_N'(%5Cxi)%24%2C%0A%20%20%20%20and%20the%20associated%20quadrature%20weights%20are%20%24w_j%20%3D%20%5Cfrac%7B2%7D%7BN(N%2B1)%5C%2CP_N(%5Cxi_j)%5E2%7D%24.%0A%0A%20%20%20%20**Nodal%20vs.%20modal.**%20The%20same%20polynomial%20can%20be%20written%20in%20two%20bases%3A%0A%20%20%20%20%24u_h%20%3D%20%5Csum_j%20u_j%5C%2C%5Cell_j(%5Cxi)%24%20(**nodal**%2C%20coefficients%20%3D%20values%20at%20the%20nodes)%20or%0A%20%20%20%20%24u_h%20%3D%20%5Csum_n%20a_n%20P_n(%5Cxi)%24%20(**modal**).%20They%20are%20related%20by%20the%20Vandermonde%20matrix%20%24V_%7Bin%7D%20%3D%20P_n(%5Cxi_i)%24%3A%20%24u%20%3D%20V%20a%24.%0A%0A%20%20%20%20The%20functions%20below%20are%20provided.%20The%20interpolation%20matrix%20uses%20the%20numerically%20stable%0A%20%20%20%20**barycentric%20formula**%20%24u_h(%5Cxi)%20%3D%20%5Cdfrac%7B%5Csum_j%20%5Cfrac%7B%5Clambda_j%7D%7B%5Cxi-%5Cxi_j%7Du_j%7D%7B%5Csum_j%20%5Cfrac%7B%5Clambda_j%7D%7B%5Cxi-%5Cxi_j%7D%7D%24%0A%20%20%20%20with%20%24%5Clambda_j%20%3D%201%2F%5Cprod_%7Bm%5Cne%20j%7D(%5Cxi_j-%5Cxi_m)%24.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(legendre%2C%20np)%3A%0A%20%20%20%20def%20lgl_nodes_weights(N)%3A%0A%20%20%20%20%20%20%20%20%22%22%22Legendre-Gauss-Lobatto%20nodes%20and%20weights%20for%20polynomial%20degree%20N%20(N%20%2B%201%20points).%22%22%22%0A%20%20%20%20%20%20%20%20P_N%20%3D%20legendre.Legendre.basis(N)%0A%20%20%20%20%20%20%20%20interior%20%3D%20np.sort(P_N.deriv().roots().real)%20if%20N%20%3E%201%20else%20np.array(%5B%5D)%0A%20%20%20%20%20%20%20%20nodes%20%3D%20np.concatenate(%5B%5B-1.0%5D%2C%20interior%2C%20%5B1.0%5D%5D)%0A%20%20%20%20%20%20%20%20weights%20%3D%202.0%20%2F%20(N%20*%20(N%20%2B%201)%20*%20P_N(nodes)%20**%202)%0A%20%20%20%20%20%20%20%20return%20nodes%2C%20weights%0A%0A%0A%20%20%20%20def%20equispaced_nodes(N)%3A%0A%20%20%20%20%20%20%20%20return%20np.linspace(-1.0%2C%201.0%2C%20N%20%2B%201)%0A%0A%0A%20%20%20%20def%20barycentric_weights(nodes)%3A%0A%20%20%20%20%20%20%20%20return%20np.array(%5B1.0%20%2F%20np.prod(xj%20-%20np.delete(nodes%2C%20j))%20for%20j%2C%20xj%20in%20enumerate(nodes)%5D)%0A%0A%0A%20%20%20%20def%20interp_matrix(nodes%2C%20x)%3A%0A%20%20%20%20%20%20%20%20%22%22%22Matrix%20I%20with%20(I%20%40%20u)%5Bk%5D%20%3D%20u_h(x%5Bk%5D)%20for%20nodal%20values%20u%20at%20%60nodes%60.%22%22%22%0A%20%20%20%20%20%20%20%20lam%20%3D%20barycentric_weights(nodes)%0A%20%20%20%20%20%20%20%20diff%20%3D%20x%5B%3A%2C%20None%5D%20-%20nodes%5BNone%2C%20%3A%5D%0A%20%20%20%20%20%20%20%20on_node%20%3D%20np.abs(diff)%20%3C%201e-14%0A%20%20%20%20%20%20%20%20diff%5Bon_node%5D%20%3D%201.0%0A%20%20%20%20%20%20%20%20I%20%3D%20lam%20%2F%20diff%0A%20%20%20%20%20%20%20%20I%20%2F%3D%20I.sum(axis%3D1%2C%20keepdims%3DTrue)%0A%20%20%20%20%20%20%20%20rows%20%3D%20on_node.any(axis%3D1)%0A%20%20%20%20%20%20%20%20I%5Brows%5D%20%3D%20on_node%5Brows%5D.astype(float)%0A%20%20%20%20%20%20%20%20return%20I%0A%0A%0A%20%20%20%20def%20legendre_modal_coeffs(nodes%2C%20u)%3A%0A%20%20%20%20%20%20%20%20%22%22%22Legendre%20coefficients%20a_n%20of%20the%20interpolating%20polynomial%3A%20u%20%3D%20V%20a.%22%22%22%0A%20%20%20%20%20%20%20%20V%20%3D%20legendre.legvander(nodes%2C%20len(nodes)%20-%201)%0A%20%20%20%20%20%20%20%20return%20np.linalg.solve(V%2C%20u)%0A%0A%20%20%20%20return%20(%0A%20%20%20%20%20%20%20%20barycentric_weights%2C%0A%20%20%20%20%20%20%20%20equispaced_nodes%2C%0A%20%20%20%20%20%20%20%20interp_matrix%2C%0A%20%20%20%20%20%20%20%20legendre_modal_coeffs%2C%0A%20%20%20%20%20%20%20%20lgl_nodes_weights%2C%0A%20%20%20%20)%0A%0A%0A%40app.cell%0Adef%20_(equispaced_nodes%2C%20interp_matrix%2C%20lgl_nodes_weights%2C%20np%2C%20plt)%3A%0A%20%20%20%20def%20runge(x)%3A%0A%20%20%20%20%20%20%20%20return%201.0%20%2F%20(1.0%20%2B%2025.0%20*%20x**2)%0A%0A%0A%20%20%20%20_N%20%3D%2016%0A%20%20%20%20_xf%20%3D%20np.linspace(-1%2C%201%2C%201000)%0A%20%20%20%20_fig%2C%20_axes%20%3D%20plt.subplots(1%2C%202%2C%20figsize%3D(10%2C%203.8)%2C%20sharey%3DTrue)%0A%20%20%20%20for%20_ax%2C%20(_name%2C%20_nodes)%20in%20zip(%0A%20%20%20%20%20%20%20%20_axes%2C%20%5B(%22equally%20spaced%22%2C%20equispaced_nodes(_N))%2C%20(%22LGL%22%2C%20lgl_nodes_weights(_N)%5B0%5D)%5D%0A%20%20%20%20)%3A%0A%20%20%20%20%20%20%20%20_ax.plot(_xf%2C%20runge(_xf)%2C%20%22k-%22%2C%20lw%3D1%2C%20label%3Dr%22%24f(x)%3D1%2F(1%2B25x%5E2)%24%22)%0A%20%20%20%20%20%20%20%20_ax.plot(_xf%2C%20interp_matrix(_nodes%2C%20_xf)%20%40%20runge(_nodes)%2C%20%22C3-%22%2C%20label%3Df%22interpolant%2C%20N%20%3D%20%7B_N%7D%22)%0A%20%20%20%20%20%20%20%20_ax.plot(_nodes%2C%20runge(_nodes)%2C%20%22o%22%2C%20color%3D%22C0%22%2C%20ms%3D4%2C%20label%3D%22nodes%22)%0A%20%20%20%20%20%20%20%20_ax.set_title(_name%20%2B%20%22%20nodes%22)%0A%20%20%20%20%20%20%20%20_ax.set_ylim(-1.5%2C%201.5)%0A%20%20%20%20%20%20%20%20_ax.set_xlabel(%22%24x%24%22)%0A%20%20%20%20%20%20%20%20_ax.legend(loc%3D%22lower%20center%22)%0A%20%20%20%20_fig.tight_layout()%0A%20%20%20%20_fig%0A%20%20%20%20return%20(runge%2C)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20This%20is%20the%20**Runge%20phenomenon**%3A%20with%20equally%20spaced%20points%2C%20raising%20the%20degree%20makes%20the%20interpolant%20oscillate%0A%20%20%20%20wildly%20near%20the%20ends%20of%20the%20interval.%20LGL%20points%20cluster%20near%20the%20endpoints%20(spacing%20%24%5Csim%201%2FN%5E2%24%20there)%2C%20which%20is%0A%20%20%20%20exactly%20what%20is%20needed%20to%20control%20the%20interpolant%20everywhere.%20This%20is%20knob%203%2C%20the%20*location*%20of%20the%20points%2C%20at%20work.%0A%0A%20%20%20%20%23%23%23%20Exercise%204%20%E2%80%94%20interpolation%20error%20versus%20%24N%24%0A%0A%20%20%20%20Complete%20%60interp_errors(f%2C%20nodes_fn%2C%20Ns%2C%20x_fine)%60%3A%20for%20each%20%24N%24%20in%20%60Ns%60%2C%20get%20the%20nodes%20from%20%60nodes_fn(N)%60%2C%20interpolate%0A%20%20%20%20the%20values%20%60f(nodes)%60%20to%20the%20fine%20grid%20%60x_fine%60%20with%20%60interp_matrix%60%2C%20and%20record%20the%20maximum%20absolute%20error.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(interp_matrix%2C%20np)%3A%0A%20%20%20%20def%20interp_errors(f%2C%20nodes_fn%2C%20Ns%2C%20x_fine)%3A%0A%20%20%20%20%20%20%20%20%22%22%22Max-norm%20interpolation%20error%20of%20f%20on%20the%20fine%20grid%2C%20for%20each%20N%20in%20Ns.%22%22%22%0A%20%20%20%20%20%20%20%20errs%20%3D%20%5B%5D%0A%20%20%20%20%20%20%20%20%23%23%23%20BEGIN%20SOLUTION%0A%20%20%20%20%20%20%20%20for%20N%20in%20Ns%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20nodes%20%3D%20nodes_fn(N)%0A%20%20%20%20%20%20%20%20%20%20%20%20errs.append(np.max(np.abs(interp_matrix(nodes%2C%20x_fine)%20%40%20f(nodes)%20-%20f(x_fine))))%0A%20%20%20%20%20%20%20%20%23%23%23%20END%20SOLUTION%0A%20%20%20%20%20%20%20%20return%20np.array(errs)%0A%0A%20%20%20%20return%20(interp_errors%2C)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(check%2C%20equispaced_nodes%2C%20interp_errors%2C%20lgl_nodes_weights%2C%20np%2C%20runge)%3A%0A%20%20%20%20_xf%20%3D%20np.linspace(-1%2C%201%2C%201001)%0A%20%20%20%20_lgl%20%3D%20lambda%20N%3A%20lgl_nodes_weights(N)%5B0%5D%0A%20%20%20%20_e_smooth%20%3D%20interp_errors(lambda%20x%3A%20np.sin(np.pi%20*%20x)%2C%20_lgl%2C%20%5B4%2C%208%2C%2016%5D%2C%20_xf)%0A%20%20%20%20ex4_done%20%3D%20_e_smooth.size%20%3D%3D%203%0A%20%20%20%20_checks%20%3D%20%5B%5D%0A%20%20%20%20if%20ex4_done%3A%0A%20%20%20%20%20%20%20%20_e_eq%20%3D%20interp_errors(runge%2C%20equispaced_nodes%2C%20%5B20%5D%2C%20_xf)%0A%20%20%20%20%20%20%20%20_checks%20%3D%20%5B%0A%20%20%20%20%20%20%20%20%20%20%20%20(np.all(np.diff(_e_smooth)%20%3C%200)%2C%20%22The%20LGL%20error%20for%20sin(pi%20x)%20should%20decrease%20with%20N.%22)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20(_e_smooth%5B-1%5D%20%3C%201e-10%2C%20f%22LGL%20error%20for%20sin(pi%20x)%20at%20N%3D16%20should%20be%20%3C%201e-10%2C%20got%20%7B_e_smooth%5B-1%5D%3A.2e%7D.%22)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20(_e_eq%5B0%5D%20%3E%2010%2C%20%22Equispaced%20interpolation%20of%20the%20Runge%20function%20at%20N%3D20%20should%20have%20a%20large%20error.%22)%2C%0A%20%20%20%20%20%20%20%20%5D%0A%20%20%20%20check(%22Exercise%204%3A%20interpolation%20errors%22%2C%20_checks)%0A%20%20%20%20return%20(ex4_done%2C)%0A%0A%0A%40app.cell%0Adef%20_(%0A%20%20%20%20equispaced_nodes%2C%0A%20%20%20%20ex4_done%2C%0A%20%20%20%20interp_errors%2C%0A%20%20%20%20lgl_nodes_weights%2C%0A%20%20%20%20mo%2C%0A%20%20%20%20np%2C%0A%20%20%20%20plt%2C%0A%20%20%20%20runge%2C%0A)%3A%0A%20%20%20%20mo.stop(not%20ex4_done%2C%20mo.callout(mo.md(%22%E2%AC%86%20Complete%20**Exercise%204**%20to%20see%20the%20convergence%20plot.%22)%2C%20kind%3D%22warn%22))%0A%20%20%20%20_Ns%20%3D%20np.arange(2%2C%2041%2C%202)%0A%20%20%20%20_xf%20%3D%20np.linspace(-1%2C%201%2C%202001)%0A%20%20%20%20_lgl%20%3D%20lambda%20N%3A%20lgl_nodes_weights(N)%5B0%5D%0A%20%20%20%20_fig%2C%20_ax%20%3D%20plt.subplots(figsize%3D(6.5%2C%204.2))%0A%20%20%20%20for%20_fname%2C%20_f%2C%20_c%20in%20%5B%0A%20%20%20%20%20%20%20%20(r%22%24%5Csin(%5Cpi%20x)%24%22%2C%20lambda%20x%3A%20np.sin(np.pi%20*%20x)%2C%20%22C0%22)%2C%0A%20%20%20%20%20%20%20%20(r%22%241%2F(1%2B25x%5E2)%24%22%2C%20runge%2C%20%22C1%22)%2C%0A%20%20%20%20%5D%3A%0A%20%20%20%20%20%20%20%20_ax.semilogy(_Ns%2C%20interp_errors(_f%2C%20_lgl%2C%20_Ns%2C%20_xf)%2C%20%22o-%22%2C%20color%3D_c%2C%20label%3Df%22%7B_fname%7D%2C%20LGL%22)%0A%20%20%20%20%20%20%20%20_ax.semilogy(_Ns%2C%20interp_errors(_f%2C%20equispaced_nodes%2C%20_Ns%2C%20_xf)%2C%20%22x--%22%2C%20color%3D_c%2C%20label%3Df%22%7B_fname%7D%2C%20equispaced%22)%0A%20%20%20%20_ax.set_ylim(1e-16%2C%201e4)%0A%20%20%20%20_ax.set_xlabel(%22polynomial%20degree%20%24N%24%22)%0A%20%20%20%20_ax.set_ylabel(%22max%20interpolation%20error%22)%0A%20%20%20%20_ax.legend()%0A%20%20%20%20_fig.tight_layout()%0A%20%20%20%20_fig%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20On%20a%20semilog%20plot%20**exponential%20convergence%20is%20a%20straight%20line**.%20The%20rate%20depends%20on%20the%20function%3A%20%24%5Csin%5Cpi%20x%24%20is%0A%20%20%20%20entire%20and%20converges%20super-exponentially%20down%20to%20round-off%2C%20while%20the%20Runge%20function%20converges%20more%20slowly%20because%20it%0A%20%20%20%20has%20poles%20at%20%24x%3D%5Cpm%20i%2F5%24%2C%20close%20to%20the%20interval.%20With%20equally%20spaced%20points%2C%20even%20%24%5Csin%5Cpi%20x%24%20eventually%20diverges%0A%20%20%20%20because%20round-off%20errors%20are%20amplified%20exponentially.%0A%0A%20%20%20%20**Modal%20coefficients%20tell%20you%20how%20well%20resolved%20you%20are.**%20For%20a%20smooth%20function%20the%20Legendre%20coefficients%20%24%7Ca_n%7C%24%20decay%0A%20%20%20%20at%20the%20same%20exponential%20rate%20as%20the%20error.%20Adaptive%20DG%20codes%20use%20exactly%20this%20decay%20to%20decide%20where%20to%20refine%0A%20%20%20%20(in%20%24h%24%20or%20in%20%24p%24).%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(legendre_modal_coeffs%2C%20lgl_nodes_weights%2C%20np%2C%20plt%2C%20runge)%3A%0A%20%20%20%20_N%20%3D%2040%0A%20%20%20%20_nodes%20%3D%20lgl_nodes_weights(_N)%5B0%5D%0A%20%20%20%20_fig%2C%20_ax%20%3D%20plt.subplots(figsize%3D(6.5%2C%204))%0A%20%20%20%20_n%20%3D%20np.arange(_N%20%2B%201)%0A%20%20%20%20%23%20Odd%20functions%20have%20only%20odd%20Legendre%20modes%20and%20even%20functions%20only%20even%20ones%2C%20so%20plot%20just%20those.%0A%20%20%20%20for%20_fname%2C%20_f%2C%20_modes%20in%20%5B%0A%20%20%20%20%20%20%20%20(r%22%24%5Csin(%5Cpi%20x)%24%20(odd%20modes)%22%2C%20lambda%20x%3A%20np.sin(np.pi%20*%20x)%2C%20_n%20%25%202%20%3D%3D%201)%2C%0A%20%20%20%20%20%20%20%20(r%22%24e%5E%7B%5Csin(%5Cpi%20x)%7D%24%20(all%20modes)%22%2C%20lambda%20x%3A%20np.exp(np.sin(np.pi%20*%20x))%2C%20_n%20%3E%3D%200)%2C%0A%20%20%20%20%20%20%20%20(r%22%241%2F(1%2B25x%5E2)%24%20(even%20modes)%22%2C%20runge%2C%20_n%20%25%202%20%3D%3D%200)%2C%0A%20%20%20%20%5D%3A%0A%20%20%20%20%20%20%20%20_a%20%3D%20legendre_modal_coeffs(_nodes%2C%20_f(_nodes))%0A%20%20%20%20%20%20%20%20_ax.semilogy(_n%5B_modes%5D%2C%20np.abs(_a%5B_modes%5D)%2C%20%22o-%22%2C%20ms%3D3%2C%20label%3D_fname)%0A%20%20%20%20_ax.set_xlabel(%22mode%20number%20%24n%24%22)%0A%20%20%20%20_ax.set_ylabel(r%22%24%7Ca_n%7C%24%22)%0A%20%20%20%20_ax.set_title(f%22Legendre%20coefficients%20from%20N%20%3D%20%7B_N%7D%20LGL%20nodes%22)%0A%20%20%20%20_ax.legend()%0A%20%20%20%20_fig.tight_layout()%0A%20%20%20%20_fig%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20---%0A%20%20%20%20%23%23%20Part%203%20%E2%80%94%20Differentiating%20the%20interpolant%0A%0A%20%20%20%20The%20derivative%20of%20the%20interpolant%20at%20the%20nodes%20is%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%5Cpartial_%5Cxi%20u_h(%5Cxi_i)%20%3D%20%5Csum_%7Bj%7D%20u_j%5C%2C%5Cell_j'(%5Cxi_i)%20%3D%20%5Csum_j%20D_%7Bij%7D%5C%2Cu_j%2C%0A%20%20%20%20%5Cqquad%20D_%7Bij%7D%20%3D%20%5Cell_j'(%5Cxi_i)%20%3D%0A%20%20%20%20%5Cbegin%7Bcases%7D%0A%20%20%20%20%5Cdfrac%7B%5Clambda_j%2F%5Clambda_i%7D%7B%5Cxi_i-%5Cxi_j%7D%2C%20%26%20i%5Cneq%20j%2C%5C%5C%5B2mm%5D%0A%20%20%20%20-%5Csum_%7Bm%5Cne%20i%7D%20D_%7Bim%7D%2C%20%26%20i%3Dj%2C%0A%20%20%20%20%5Cend%7Bcases%7D%0A%20%20%20%20%24%24%0A%0A%20%20%20%20and%20physical%20derivatives%20follow%20from%20the%20chain%20rule%2C%20%24%5Cpartial_x%20%3D%20J%5E%7B-1%7D%5Cpartial_%5Cxi%24.%20This%20is%20**exact**%20for%20the%0A%20%20%20%20polynomial%20%24u_h%24%3A%20all%20of%20the%20error%20comes%20from%20how%20well%20%24u_h%24%20approximates%20the%20true%20solution.%20Because%20the%0A%20%20%20%20diagonal%20is%20chosen%20so%20that%20each%20row%20sums%20to%20zero%2C%20%24D%24%20differentiates%20constants%20exactly%20even%20in%20floating%20point.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(barycentric_weights%2C%20lgl_nodes_weights%2C%20mo%2C%20np)%3A%0A%20%20%20%20def%20diff_matrix(nodes)%3A%0A%20%20%20%20%20%20%20%20%22%22%22Lagrange%20differentiation%20matrix%20D%5Bi%2C%20j%5D%20%3D%20l_j'(xi_i).%22%22%22%0A%20%20%20%20%20%20%20%20lam%20%3D%20barycentric_weights(nodes)%0A%20%20%20%20%20%20%20%20diff%20%3D%20nodes%5B%3A%2C%20None%5D%20-%20nodes%5BNone%2C%20%3A%5D%0A%20%20%20%20%20%20%20%20np.fill_diagonal(diff%2C%201.0)%0A%20%20%20%20%20%20%20%20D%20%3D%20(lam%5BNone%2C%20%3A%5D%20%2F%20lam%5B%3A%2C%20None%5D)%20%2F%20diff%0A%20%20%20%20%20%20%20%20np.fill_diagonal(D%2C%200.0)%0A%20%20%20%20%20%20%20%20np.fill_diagonal(D%2C%20-D.sum(axis%3D1))%0A%20%20%20%20%20%20%20%20return%20D%0A%0A%0A%20%20%20%20%23%20Sanity%20check%3A%20D%20differentiates%20polynomials%20of%20degree%20%3C%3D%20N%20exactly.%0A%20%20%20%20_nodes%20%3D%20lgl_nodes_weights(8)%5B0%5D%0A%20%20%20%20_max_err%20%3D%20max(%0A%20%20%20%20%20%20%20%20np.max(np.abs(diff_matrix(_nodes)%20%40%20_nodes**p%20-%20p%20*%20_nodes%20**%20max(p%20-%201%2C%200)))%20for%20p%20in%20range(9)%0A%20%20%20%20)%0A%20%20%20%20mo.md(f%22Sanity%20check%3A%20max%20error%20of%20%24D%24%20applied%20to%20%24%5C%5Cxi%5Ep%24%2C%20%24p%5C%5Cle%208%24%2C%20with%20%24N%3D8%24%3A%20**%7B_max_err%3A.1e%7D**%22)%0A%20%20%20%20return%20(diff_matrix%2C)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%23%20Exercise%205%20%E2%80%94%20spectral%20derivative%20of%20an%20analytic%20function%0A%0A%20%20%20%20Take%20%24f(x)%20%3D%20e%5E%7B%5Csin%20x%7D%24%20on%20%24%5B0%2C%202%5Cpi%5D%24%20as%20a%20**single%20element**%20(%24x%20%3D%20%5Cpi(1%2B%5Cxi)%24%2C%20%24J%3D%5Cpi%24).%0A%20%20%20%20Complete%20%60derivative_errors(Ns)%60%3A%20for%20each%20%24N%24%2C%20compute%20%24f%24%20at%20the%20LGL%20nodes%2C%20differentiate%20with%20%60diff_matrix%60%0A%20%20%20%20(remember%20the%20%241%2FJ%24)%2C%20and%20record%20the%20maximum%20error%20at%20the%20nodes%20compared%20with%20%24f'(x)%20%3D%20%5Ccos%20x%5C%2Ce%5E%7B%5Csin%20x%7D%24.%0A%0A%20%20%20%20The%20cell%20below%20then%20compares%20this%20with%20the%204th-order%20finite%20differences%20of%20Lecture%201%2C%20which%20*knows*%20the%20function%20is%0A%20%20%20%20periodic%2C%20using%20the%20same%20number%20of%20points.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(diff_matrix%2C%20lgl_nodes_weights%2C%20np)%3A%0A%20%20%20%20def%20derivative_errors(Ns)%3A%0A%20%20%20%20%20%20%20%20%22%22%22Max%20error%20at%20the%20nodes%20of%20the%20LGL%20derivative%20of%20exp(sin%20x)%20on%20one%20element%20%5B0%2C%202pi%5D.%22%22%22%0A%20%20%20%20%20%20%20%20errs%20%3D%20%5B%5D%0A%20%20%20%20%20%20%20%20%23%23%23%20BEGIN%20SOLUTION%0A%20%20%20%20%20%20%20%20for%20N%20in%20Ns%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20xi%20%3D%20lgl_nodes_weights(N)%5B0%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20x%20%3D%20np.pi%20*%20(1%20%2B%20xi)%0A%20%20%20%20%20%20%20%20%20%20%20%20dfdx%20%3D%20diff_matrix(xi)%20%40%20np.exp(np.sin(x))%20%2F%20np.pi%0A%20%20%20%20%20%20%20%20%20%20%20%20errs.append(np.max(np.abs(dfdx%20-%20np.cos(x)%20*%20np.exp(np.sin(x)))))%0A%20%20%20%20%20%20%20%20%23%23%23%20END%20SOLUTION%0A%20%20%20%20%20%20%20%20return%20np.array(errs)%0A%0A%20%20%20%20return%20(derivative_errors%2C)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(check%2C%20derivative_errors)%3A%0A%20%20%20%20_e%20%3D%20derivative_errors(%5B8%2C%2016%2C%2048%5D)%0A%20%20%20%20ex5_done%20%3D%20_e.size%20%3D%3D%203%0A%20%20%20%20_checks%20%3D%20%5B%5D%0A%20%20%20%20if%20ex5_done%3A%0A%20%20%20%20%20%20%20%20_checks%20%3D%20%5B%0A%20%20%20%20%20%20%20%20%20%20%20%20(_e%5B0%5D%20%3E%20_e%5B1%5D%20%3E%20_e%5B2%5D%2C%20%22The%20error%20should%20decrease%20with%20N.%22)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20(_e%5B2%5D%20%3C%201e-10%2C%20f%22At%20N%20%3D%2048%20the%20error%20should%20be%20%3C%201e-10%2C%20got%20%7B_e%5B2%5D%3A.2e%7D%20(did%20you%20divide%20by%20J%3F).%22)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20(_e%5B0%5D%20%3E%201e-4%2C%20%22At%20N%20%3D%208%20the%20error%20looks%20too%20small%3B%20are%20you%20comparing%20with%20the%20right%20derivative%3F%22)%2C%0A%20%20%20%20%20%20%20%20%5D%0A%20%20%20%20check(%22Exercise%205%3A%20spectral%20derivative%22%2C%20_checks)%0A%20%20%20%20return%20(ex5_done%2C)%0A%0A%0A%40app.cell%0Adef%20_(T%2C%20derivative_errors%2C%20ex5_done%2C%20mo%2C%20np%2C%20plt)%3A%0A%20%20%20%20mo.stop(not%20ex5_done%2C%20mo.callout(mo.md(%22%E2%AC%86%20Complete%20**Exercise%205**%20first.%22)%2C%20kind%3D%22warn%22))%0A%0A%0A%20%20%20%20def%20_fd4_error(n_points)%3A%0A%20%20%20%20%20%20%20%20_dx%20%3D%20T%20%2F%20n_points%0A%20%20%20%20%20%20%20%20_x%20%3D%20np.arange(n_points)%20*%20_dx%0A%20%20%20%20%20%20%20%20_u%20%3D%20np.exp(np.sin(_x))%0A%20%20%20%20%20%20%20%20_du%20%3D%20(np.roll(_u%2C%202)%20-%208%20*%20np.roll(_u%2C%201)%20%2B%208%20*%20np.roll(_u%2C%20-1)%20-%20np.roll(_u%2C%20-2))%20%2F%20(12%20*%20_dx)%0A%20%20%20%20%20%20%20%20return%20np.max(np.abs(_du%20-%20np.cos(_x)%20*%20_u))%0A%0A%0A%20%20%20%20_Ns%20%3D%20np.arange(4%2C%2061%2C%204)%0A%20%20%20%20_fig%2C%20(_ax0%2C%20_ax1)%20%3D%20plt.subplots(1%2C%202%2C%20figsize%3D(10%2C%204)%2C%20sharey%3DTrue)%0A%20%20%20%20for%20_ax%20in%20(_ax0%2C%20_ax1)%3A%0A%20%20%20%20%20%20%20%20_ax.plot(_Ns%20%2B%201%2C%20derivative_errors(_Ns)%2C%20%22o-%22%2C%20label%3D%22LGL%20spectral%2C%20one%20element%22)%0A%20%20%20%20%20%20%20%20_ax.plot(_Ns%20%2B%201%2C%20%5B_fd4_error(n)%20for%20n%20in%20_Ns%20%2B%201%5D%2C%20%22s-%22%2C%20label%3D%22FD4%2C%20periodic%20grid%22)%0A%20%20%20%20%20%20%20%20_ax.set_xlabel(%22number%20of%20points%22)%0A%20%20%20%20%20%20%20%20_ax.legend()%0A%20%20%20%20_ax0.set_yscale(%22log%22)%0A%20%20%20%20_ax0.set_title(%22semilog%3A%20spectral%20is%20a%20straight%20line%22)%0A%20%20%20%20_ax1.set_xscale(%22log%22)%0A%20%20%20%20_ax1.set_yscale(%22log%22)%0A%20%20%20%20_ax1.set_xticks(%5B5%2C%2010%2C%2020%2C%2040%2C%2060%5D%2C%20%5B%225%22%2C%20%2210%22%2C%20%2220%22%2C%20%2240%22%2C%20%2260%22%5D)%0A%20%20%20%20_ax1.minorticks_off()%0A%20%20%20%20_ax1.set_title(%22log-log%3A%20FD4%20is%20a%20straight%20line%22)%0A%20%20%20%20_ax0.set_ylabel(r%22max%20error%20in%20%24f'%24%22)%0A%20%20%20%20_fig.tight_layout()%0A%20%20%20%20_fig%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%23%20Why%20nodal%20methods%20make%20nonlinear%20terms%20easy%0A%0A%20%20%20%20Real%20PDEs%20have%20nonlinear%20fluxes%20and%20sources%3A%20%24u%5E2%24%20in%20Burgers'%20equation%2C%20%24%5Crho%20v%5E2%20%2B%20p(%5Crho%2C%5Cepsilon)%24%20in%20hydrodynamics%2C%0A%20%20%20%20products%20of%20metric%20components%20in%20general%20relativity.%20In%20a%20**nodal**%20method%20we%20evaluate%20them%20**pointwise%20at%20the%20nodes**%3A%0A%20%20%20%20%24F(u)_j%20%3D%20F(u_j)%24.%20That%20costs%20%24%5Cmathcal%20O(N)%24%20operations%20and%20works%20for%20*any*%20function%2C%20including%20%24%5Csqrt%7Bu%7D%24%2C%20%24e%5Eu%24%2C%0A%20%20%20%20tabulated%20equations%20of%20state%2C%20or%20conditionals.%0A%0A%20%20%20%20In%20a%20**modal**%20method%20the%20unknowns%20are%20coefficients%20%24a_n%24.%20The%20product%20of%20two%20degree-%24N%24%20series%20is%20a%20degree-%242N%24%20series%0A%20%20%20%20whose%20coefficients%20require%20%24%5Cmathcal%20O(N%5E2)%24%20work%20(a%20%22convolution%22%20of%20Legendre%20coefficients)%2C%20and%20there%20is%20no%20finite%0A%20%20%20%20formula%20at%20all%20for%20%24e%5Eu%24%20or%20%24%5Csqrt%20u%24%3A%20you%20would%20transform%20to%20point%20values%2C%20evaluate%2C%20and%20transform%20back%2C%20which%20is%20exactly%0A%20%20%20%20what%20the%20nodal%20method%20does%20from%20the%20start.%20(Pointwise%20evaluation%20of%20a%20nonlinear%20term%20does%20introduce%20a%20small%0A%20%20%20%20*aliasing*%20error%2C%20which%20converges%20away%20at%20the%20same%20spectral%20rate%20for%20smooth%20solutions.)%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(%0A%20%20%20%20interp_matrix%2C%0A%20%20%20%20legendre%2C%0A%20%20%20%20legendre_modal_coeffs%2C%0A%20%20%20%20lgl_nodes_weights%2C%0A%20%20%20%20mo%2C%0A%20%20%20%20np%2C%0A)%3A%0A%20%20%20%20import%20timeit%0A%0A%20%20%20%20_N%20%3D%2032%0A%20%20%20%20_xi%20%3D%20lgl_nodes_weights(_N)%5B0%5D%0A%20%20%20%20_u%20%3D%20np.exp(np.sin(np.pi%20*%20_xi))%20%20%23%20nodal%20values%20of%20u%0A%20%20%20%20_a%20%3D%20legendre_modal_coeffs(_xi%2C%20_u)%20%20%23%20modal%20coefficients%20of%20the%20same%20polynomial%0A%20%20%20%20_xf%20%3D%20np.linspace(-1%2C%201%2C%201001)%0A%20%20%20%20_exact%20%3D%20np.exp(2%20*%20np.sin(np.pi%20*%20_xf))%0A%0A%20%20%20%20_nodal_sq%20%3D%20_u%20*%20_u%20%20%23%20nodal%3A%20pointwise%20product%2C%20O(N)%0A%20%20%20%20_modal_sq%20%3D%20legendre.legmul(_a%2C%20_a)%5B%3A%20_N%20%2B%201%5D%20%20%23%20modal%3A%20series%20product%2C%20then%20truncate%20to%20degree%20N%0A%0A%20%20%20%20_err_nodal%20%3D%20np.max(np.abs(interp_matrix(_xi%2C%20_xf)%20%40%20_nodal_sq%20-%20_exact))%0A%20%20%20%20_err_modal%20%3D%20np.max(np.abs(legendre.legval(_xf%2C%20_modal_sq)%20-%20_exact))%0A%20%20%20%20_t_nodal%20%3D%20timeit.timeit(lambda%3A%20_u%20*%20_u%2C%20number%3D2000)%20%2F%202000%0A%20%20%20%20_t_modal%20%3D%20timeit.timeit(lambda%3A%20legendre.legmul(_a%2C%20_a)%5B%3A%20_N%20%2B%201%5D%2C%20number%3D2000)%20%2F%202000%0A%0A%20%20%20%20mo.md(%0A%20%20%20%20%20%20%20%20f%22%22%22%0A%20%20%20%20Computing%20%24u%5E2%24%20for%20%24u%20%3D%20e%5E%7B%7B%5C%5Csin%5C%5Cpi%20x%7D%7D%24%20with%20%24N%20%3D%20%7B_N%7D%24%3A%0A%0A%20%20%20%20%7C%20method%20%7C%20max%20error%20of%20%24u%5E2%24%20%7C%20time%20per%20evaluation%20%7C%0A%20%20%20%20%7C---%7C---%7C---%7C%0A%20%20%20%20%7C%20nodal%20(pointwise%20%60u%20*%20u%60)%20%7C%20%7B_err_nodal%3A.1e%7D%20%7C%20%7B_t_nodal%20*%201e6%3A.2f%7D%20%C2%B5s%20%7C%0A%20%20%20%20%7C%20modal%20(%60legmul%60%20%2B%20truncation)%20%7C%20%7B_err_modal%3A.1e%7D%20%7C%20%7B_t_modal%20*%201e6%3A.2f%7D%20%C2%B5s%20%7C%0A%0A%20%20%20%20Same%20accuracy%2C%20but%20the%20modal%20product%20is%20more%20expensive%20and%20more%20complicated%2C%20and%20it%20doesn't%20generalise%20to%20%24e%5Eu%24.%0A%20%20%20%20%22%22%22%0A%20%20%20%20)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20---%0A%20%20%20%20%23%23%20Part%204%20%E2%80%94%20Boundary%20corrections%20(numerical%20fluxes)%0A%0A%20%20%20%20Inside%20each%20element%20the%20solution%20is%20a%20smooth%20polynomial%3B%20at%20a%20face%20we%20have%20two%20values%2C%20%24u%5E%7B%5Crm%20int%7D%24%20and%20%24u%5E%7B%5Crm%20ext%7D%24.%0A%20%20%20%20The%20numerical%20flux%20%24G(u%5E%7B%5Crm%20int%7D%2C%20u%5E%7B%5Crm%20ext%7D%3B%20n)%24%20decides%20how%20they%20interact.%20It%20must%20have%20these%20properties%3A%0A%0A%20%20%20%201.%20**Consistency%3A**%20%24G(u%2C%20u%3B%20n)%20%3D%20n%20F(u)%24.%20If%20the%20solution%20is%20continuous%20the%20correction%20%24G%20-%20nF%5E%7B%5Crm%20int%7D%24%20vanishes%2C%0A%20%20%20%20%20%20%20and%20the%20scheme%20reduces%20to%20the%20PDE.%0A%20%20%20%202.%20**Conservation%3A**%20%24G(u%5E%7B%5Crm%20int%7D%2C%20u%5E%7B%5Crm%20ext%7D%3B%20n)%20%3D%20-G(u%5E%7B%5Crm%20ext%7D%2C%20u%5E%7B%5Crm%20int%7D%3B%20-n)%24.%20What%20leaves%20one%20element%20enters%20its%0A%20%20%20%20%20%20%20neighbour%2C%20so%20conserved%20quantities%20(mass%2C%20energy%2C%20charge%2C%20%E2%80%A6)%20are%20conserved%20exactly%20by%20the%20discrete%20scheme.%0A%20%20%20%203.%20**Stability%20(dissipation)%3A**%20%24G%24%20should%20depend%20on%20the%20jump%20%24%5Bu%5D%20%3D%20u%5E%7B%5Crm%20ext%7D-u%5E%7B%5Crm%20int%7D%24%20in%20a%20way%20that%20removes%0A%20%20%20%20%20%20%20energy%20from%20the%20discontinuities.%0A%0A%20%20%20%20%23%23%23%20Why%20dissipation%3A%20an%20energy%20estimate%0A%0A%20%20%20%20Take%20advection%2C%20%24%5Cpartial_t%20u%20%2B%20a%5C%2C%5Cpartial_x%20u%20%3D%200%24%2C%20with%20energy%20%24E%20%3D%20%5Ctfrac12%5Csum_k%5Cint_%7B%5COmega_k%7Du_h%5E2%5C%2Cdx%24%2C%20and%20the%20flux%0A%20%20%20%20family%20%24G%20%3D%20n%5Cleft(a%5C%2C%5C%7Bu%5C%7D%20-%20%5Ctfrac%7B%5Calpha%7D%7B2%7D%5C%2Cn%5C%2C%5Bu%5D%5Cright)%24%2C%20where%20%24%5C%7Bu%5C%7D%24%20is%20the%20average%20of%20the%20two%20sides.%0A%20%20%20%20Multiplying%20the%20strong%20form%20by%20%24u_h%24%2C%20integrating%2C%20and%20adding%20the%20contributions%20of%20the%20two%20elements%20that%20share%20each%20face%20gives%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%5Cfrac%7BdE%7D%7Bdt%7D%20%3D%20-%5Cfrac%7B%5Calpha%7D%7B2%7D%5Csum_%7B%5Crm%20faces%7D%5Bu%5D%5E2%20.%0A%20%20%20%20%24%24%0A%0A%20%20%20%20*%20**Central%20flux**%20(%24%5Calpha%3D0%24)%3A%20energy%20is%20conserved%20exactly%2C%20and%20nothing%20damps%20under-resolved%20or%20spurious%20modes.%0A%20%20%20%20*%20**Upwind%20flux**%20(%24%5Calpha%3D%7Ca%7C%24)%3A%20energy%20decreases%20in%20proportion%20to%20the%20jumps%20squared.%0A%0A%20%20%20%20For%20a%20smooth%2C%20well-resolved%20solution%20the%20jumps%20are%20tiny%20(%24%5Csim%5CDelta%20x%5E%7BN%2B1%7D%24)%2C%20so%20this%20dissipation%20does%20not%20hurt%20accuracy.%0A%20%20%20%20It%20acts%20only%20on%20under-resolved%20features%2C%20where%20it%20is%20exactly%20what%20we%20want.%20For%20nonlinear%20problems%20this%20dissipation%20is%0A%20%20%20%20essential%20for%20stability.%0A%0A%20%20%20%20%23%23%23%20Characteristic%20(upwind)%20correction%20for%20the%20scalar%20wave%0A%0A%20%20%20%20For%20our%20system%20the%20**characteristic%20fields**%20along%20the%20outward%20normal%20%24n%24%20are%0A%0A%20%20%20%20%24%24%0A%20%20%20%20w%5E%5CPsi%20%3D%20%5CPsi%5C%20%5C%20(%5Clambda%20%3D%200)%2C%20%5Cqquad%20w%5E%5Cpm%20%3D%20%5CPi%20%5Cpm%20n%5C%2C%5CPhi%20-%20%5Cgamma_2%5CPsi%5C%20%5C%20(%5Clambda%20%3D%20%5Cpm%201)%2C%0A%20%20%20%20%24%24%0A%0A%20%20%20%20with%20inverse%0A%20%20%20%20%24%5CPsi%20%3D%20w%5E%5CPsi%2C%5C%20%5C%20%5CPi%20%3D%20%5Ctfrac12(w%5E%2B%2Bw%5E-)%20%2B%20%5Cgamma_2%20w%5E%5CPsi%2C%5C%20%5C%20%5CPhi%20%3D%20%5Ctfrac%7Bn%7D%7B2%7D(w%5E%2B-w%5E-)%24.%0A%20%20%20%20Only%20fields%20with%20negative%20speed%20(%24%5Clambda%3C0%24)%20**enter**%20the%20element.%20The%20upwind%20correction%20replaces%20exactly%20those%20fields%0A%20%20%20%20by%20the%20neighbour's%20values%3A%0A%0A%20%20%20%20%24%24%0A%20%20%20%20G%20-%20nF%5E%7B%5Crm%20int%7D%20%3D%20T%5C%2C%5CLambda%5E-%5C%2C%5Cbig(w%5E%7B%5Crm%20ext%7D%20-%20w%5E%7B%5Crm%20int%7D%5Cbig)%2C%0A%20%20%20%20%24%24%0A%0A%20%20%20%20where%20%24%5CLambda%5E-%24%20keeps%20only%20the%20negative%20speeds%20and%20%24T%24%20is%20the%20matrix%20of%20the%20inverse%20transformation%20above%0A%20%20%20%20(both%20%24w%5E%7B%5Crm%20ext%7D%24%20and%20%24w%5E%7B%5Crm%20int%7D%24%20are%20computed%20with%20the%20*same*%2C%20interior%2C%20normal%20%24n%24).%0A%0A%20%20%20%20%23%23%23%20Rusanov%20(local%20Lax%E2%80%93Friedrichs)%20correction%0A%0A%20%20%20%20A%20simpler%2C%20more%20dissipative%20choice%20that%20needs%20no%20characteristic%20decomposition%2C%20only%20the%20largest%20speed%20%24%5Clambda_%7B%5Cmax%7D%24%3A%0A%0A%20%20%20%20%24%24%0A%20%20%20%20G%20%3D%20%5Ctfrac12%20n%5Cleft(F%5E%7B%5Crm%20int%7D%2BF%5E%7B%5Crm%20ext%7D%5Cright)%20-%20%5Ctfrac12%5Clambda_%7B%5Cmax%7D%5Cleft(u%5E%7B%5Crm%20ext%7D-u%5E%7B%5Crm%20int%7D%5Cright)%0A%20%20%20%20%5Cquad%5CLongrightarrow%5Cquad%0A%20%20%20%20G%20-%20nF%5E%7B%5Crm%20int%7D%20%3D%20%5Ctfrac12%20n%5Cleft(F%5E%7B%5Crm%20ext%7D-F%5E%7B%5Crm%20int%7D%5Cright)%20-%20%5Ctfrac12%5Clambda_%7B%5Cmax%7D%5Cleft(u%5E%7B%5Crm%20ext%7D-u%5E%7B%5Crm%20int%7D%5Cright).%0A%20%20%20%20%24%24%0A%0A%20%20%20%20Here%20%24%5Clambda_%7B%5Cmax%7D%3D1%24.%20For%20systems%20where%20computing%20the%20characteristic%20fields%20is%20hard%20(e.g.%20relativistic%0A%20%20%20%20magnetohydrodynamics)%2C%20Rusanov%20and%20HLL-type%20fluxes%20are%20the%20workhorses.%0A%0A%20%20%20%20*(In%20curved%20spacetimes%20or%20on%20moving%20meshes%20the%20normals%20and%20speeds%20differ%20on%20the%20two%20sides%20of%20a%20face%3B%20one%20then%20uses%0A%20%20%20%20%24T%5E%7B%5Crm%20ext%7D%5CLambda%5E%7B-%2C%5Crm%20ext%7Dw%5E%7B%5Crm%20ext%7D%20-%20T%5E%7B%5Crm%20int%7D%5CLambda%5E%7B-%2C%5Crm%20int%7Dw%5E%7B%5Crm%20int%7D%24%20to%20keep%20the%20conservation%20property.)*%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20---%0A%20%20%20%20%23%23%20Part%205%20%E2%80%94%20A%20multi-element%20DG%20solver%20for%20the%20scalar%20wave%0A%0A%20%20%20%20**Data%20layout.**%20The%20solution%20is%20stored%20as%20an%20array%20%60u%60%20of%20shape%20%60(3%2C%20K%2C%20N%20%2B%201)%60%3A%20variable%20%24(%5CPsi%2C%5CPi%2C%5CPhi)%24%20%C3%97%20element%20%C3%97%20node.%0A%20%20%20%20Face%20values%20are%20then%20arrays%20of%20shape%20%60(3%2C%20K)%60%3A%20%60u%5B%3A%2C%20%3A%2C%200%5D%60%20is%20every%20element's%20left%20face%20and%20%60u%5B%3A%2C%20%3A%2C%20-1%5D%60%20its%20right%20face.%0A%20%20%20%20Periodicity%20means%20the%20neighbour%20to%20the%20right%20of%20element%20%24k%24%20is%20element%20%24k%2B1%24%20modulo%20%24K%24%2C%20i.e.%20%60np.roll%60%20along%20axis%201.%0A%0A%20%20%20%20Provided%20below%3A%20the%20grid%2C%20%60flux%60%2C%20%60source%60%2C%20the%20**central**%20correction%20(as%20an%20example%20of%20the%20interface)%2C%20and%20%60dg_rhs%60%2C%0A%20%20%20%20which%20assembles%20the%20scheme%20including%20the%20lifting.%20You%20will%20write%0A%0A%20%20%20%20*%20**Exercise%206a**%20%60volume_terms(u%2C%20grid%2C%20gamma2)%60%3A%20%24-%5Cfrac1J%5Csum_j%20D_%7Bij%7DF(u_j)%20%2B%20S(u_i)%24%20for%20all%20variables%2C%20elements%20and%20nodes%3B%0A%20%20%20%20*%20**Exercise%206b**%20%60boundary_correction_upwind(u_int%2C%20u_ext%2C%20n%2C%20gamma2)%60%3A%20%24G%20-%20nF%5E%7B%5Crm%20int%7D%24%20from%20the%20characteristic%20fields%3B%0A%20%20%20%20*%20**Exercise%206c**%20%60boundary_correction_rusanov(u_int%2C%20u_ext%2C%20n%2C%20gamma2)%60.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(T%2C%20diff_matrix%2C%20lgl_nodes_weights%2C%20np)%3A%0A%20%20%20%20from%20types%20import%20SimpleNamespace%0A%0A%0A%20%20%20%20def%20make_grid(K%2C%20N)%3A%0A%20%20%20%20%20%20%20%20%22%22%22K%20equal%20elements%20covering%20%5B0%2C%202pi)%2C%20each%20with%20N%20%2B%201%20LGL%20nodes.%22%22%22%0A%20%20%20%20%20%20%20%20xi%2C%20w%20%3D%20lgl_nodes_weights(N)%0A%20%20%20%20%20%20%20%20dx%20%3D%20T%20%2F%20K%0A%20%20%20%20%20%20%20%20x%20%3D%20(np.arange(K)%20*%20dx)%5B%3A%2C%20None%5D%20%2B%20(1%20%2B%20xi)%5BNone%2C%20%3A%5D%20*%20dx%20%2F%202%20%20%23%20shape%20(K%2C%20N%20%2B%201)%0A%20%20%20%20%20%20%20%20return%20SimpleNamespace(K%3DK%2C%20N%3DN%2C%20xi%3Dxi%2C%20w%3Dw%2C%20D%3Ddiff_matrix(xi)%2C%20J%3Ddx%20%2F%202%2C%20x%3Dx)%0A%0A%0A%20%20%20%20def%20flux(u%2C%20gamma2)%3A%0A%20%20%20%20%20%20%20%20psi%2C%20pi%2C%20phi%20%3D%20u%0A%20%20%20%20%20%20%20%20return%20np.stack(%5Bnp.zeros_like(psi)%2C%20phi%2C%20pi%20-%20gamma2%20*%20psi%5D)%0A%0A%0A%20%20%20%20def%20source(u%2C%20gamma2)%3A%0A%20%20%20%20%20%20%20%20psi%2C%20pi%2C%20phi%20%3D%20u%0A%20%20%20%20%20%20%20%20return%20np.stack(%5B-pi%2C%20np.zeros_like(pi)%2C%20-gamma2%20*%20phi%5D)%0A%0A%0A%20%20%20%20def%20boundary_correction_central(u_int%2C%20u_ext%2C%20n%2C%20gamma2)%3A%0A%20%20%20%20%20%20%20%20%22%22%22G%20-%20n%20F_int%20for%20the%20central%20flux%20G%20%3D%20n%20(F_int%20%2B%20F_ext)%20%2F%202.%20Face%20arrays%20have%20shape%20(3%2C%20K).%22%22%22%0A%20%20%20%20%20%20%20%20return%200.5%20*%20n%20*%20(flux(u_ext%2C%20gamma2)%20-%20flux(u_int%2C%20gamma2))%0A%0A%20%20%20%20return%20boundary_correction_central%2C%20flux%2C%20make_grid%2C%20source%0A%0A%0A%40app.cell%0Adef%20_(hint)%3A%0A%20%20%20%20hint(%0A%20%20%20%20%20%20%20%20r%22%60grid.D%60%20has%20shape%20%60(N%2B1%2C%20N%2B1)%60%20and%20%60F%20%3D%20flux(u%2C%20gamma2)%60%20has%20shape%20%60(3%2C%20K%2C%20N%2B1)%60.%20%22%0A%20%20%20%20%20%20%20%20r%22%60F%20%40%20grid.D.T%60%20computes%20%24%5Csum_j%20D_%7Bij%7DF_j%24%20along%20the%20last%20axis%20for%20every%20variable%20and%20element%20at%20once.%22%2C%0A%20%20%20%20%20%20%20%20r%22For%20the%20upwind%20correction%2C%20compute%20%24w%5E-%24%20on%20both%20sides%2C%20%24w%5E-%20%3D%20%5CPi%20-%20n%5CPhi%20-%20%5Cgamma_2%5CPsi%24%20(same%20%60n%60%20for%20both)%2C%20%22%0A%20%20%20%20%20%20%20%20r%22and%20the%20jump%20%24%5CDelta%20w%5E-%20%3D%20w%5E%7B-%2C%5Crm%20ext%7D%20-%20w%5E%7B-%2C%5Crm%20int%7D%24.%20Only%20%24%5Clambda_-%20%3D%20-1%24%20is%20negative.%22%2C%0A%20%20%20%20%20%20%20%20r%22%24T%5CLambda%5E-%5CDelta%20w%24%20with%20only%20%24w%5E-%24%20incoming%3A%20the%20%24%5CPsi%24%20component%20is%20%240%24%3B%20the%20%24%5CPi%24%20component%20is%20%22%0A%20%20%20%20%20%20%20%20r%22%24%5Ctfrac12%5Clambda_-%5CDelta%20w%5E-%24%3B%20the%20%24%5CPhi%24%20component%20is%20%24-%5Ctfrac%20n2%5Clambda_-%5CDelta%20w%5E-%24.%20%22%0A%20%20%20%20%20%20%20%20r%22Stack%20them%20with%20%60np.stack(%5B...%5D)%60%20to%20get%20shape%20%60(3%2C%20K)%60.%22%2C%0A%20%20%20%20)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(flux%2C%20source)%3A%0A%20%20%20%20def%20volume_terms(u%2C%20grid%2C%20gamma2)%3A%0A%20%20%20%20%20%20%20%20%22%22%22-(1%2FJ)%20D%20F(u)%20%2B%20S(u)%2C%20shape%20(3%2C%20K%2C%20N%20%2B%201).%22%22%22%0A%20%20%20%20%20%20%20%20%23%23%23%20BEGIN%20SOLUTION%0A%20%20%20%20%20%20%20%20dt_u_volume%20%3D%20-(flux(u%2C%20gamma2)%20%40%20grid.D.T)%20%2F%20grid.J%20%2B%20source(u%2C%20gamma2)%0A%20%20%20%20%20%20%20%20%23%23%23%20END%20SOLUTION%0A%20%20%20%20%20%20%20%20return%20dt_u_volume%0A%0A%20%20%20%20return%20(volume_terms%2C)%0A%0A%0A%40app.cell%0Adef%20_(np)%3A%0A%20%20%20%20def%20boundary_correction_upwind(u_int%2C%20u_ext%2C%20n%2C%20gamma2)%3A%0A%20%20%20%20%20%20%20%20%22%22%22Upwind%20G%20-%20n%20F_int%20%3D%20T%20Lambda%5E-%20(w_ext%20-%20w_int).%20Face%20arrays%20have%20shape%20(3%2C%20K)%3B%20n%20%3D%20%2B1%20or%20-1.%22%22%22%0A%20%20%20%20%20%20%20%20%23%23%23%20BEGIN%20SOLUTION%0A%20%20%20%20%20%20%20%20w_minus_int%20%3D%20u_int%5B1%5D%20-%20n%20*%20u_int%5B2%5D%20-%20gamma2%20*%20u_int%5B0%5D%0A%20%20%20%20%20%20%20%20w_minus_ext%20%3D%20u_ext%5B1%5D%20-%20n%20*%20u_ext%5B2%5D%20-%20gamma2%20*%20u_ext%5B0%5D%0A%20%20%20%20%20%20%20%20_lam_minus%20%3D%20-1.0%0A%20%20%20%20%20%20%20%20_dw%20%3D%20w_minus_ext%20-%20w_minus_int%0A%20%20%20%20%20%20%20%20correction%20%3D%20np.stack(%5Bnp.zeros_like(_dw)%2C%200.5%20*%20_lam_minus%20*%20_dw%2C%20-0.5%20*%20n%20*%20_lam_minus%20*%20_dw%5D)%0A%20%20%20%20%20%20%20%20%23%23%23%20END%20SOLUTION%0A%20%20%20%20%20%20%20%20return%20correction%0A%0A%20%20%20%20return%20(boundary_correction_upwind%2C)%0A%0A%0A%40app.cell%0Adef%20_(flux)%3A%0A%20%20%20%20def%20boundary_correction_rusanov(u_int%2C%20u_ext%2C%20n%2C%20gamma2%2C%20lambda_max%3D1.0)%3A%0A%20%20%20%20%20%20%20%20%22%22%22Rusanov%20G%20-%20n%20F_int.%20Face%20arrays%20have%20shape%20(3%2C%20K)%3B%20n%20%3D%20%2B1%20or%20-1.%22%22%22%0A%20%20%20%20%20%20%20%20%23%23%23%20BEGIN%20SOLUTION%0A%20%20%20%20%20%20%20%20correction%20%3D%200.5%20*%20n%20*%20(flux(u_ext%2C%20gamma2)%20-%20flux(u_int%2C%20gamma2))%20-%200.5%20*%20lambda_max%20*%20(u_ext%20-%20u_int)%0A%20%20%20%20%20%20%20%20%23%23%23%20END%20SOLUTION%0A%20%20%20%20%20%20%20%20return%20correction%0A%0A%20%20%20%20return%20(boundary_correction_rusanov%2C)%0A%0A%0A%40app.cell%0Adef%20_(np%2C%20volume_terms)%3A%0A%20%20%20%20def%20dg_rhs(t%2C%20y%2C%20grid%2C%20gamma2%2C%20correction)%3A%0A%20%20%20%20%20%20%20%20%22%22%22Semi-discrete%20strong-form%20nodal%20DG%20right-hand%20side%20for%20the%20scalar%20wave.%22%22%22%0A%20%20%20%20%20%20%20%20u%20%3D%20y.reshape(3%2C%20grid.K%2C%20grid.N%20%2B%201)%0A%20%20%20%20%20%20%20%20dt_u%20%3D%20np.array(volume_terms(u%2C%20grid%2C%20gamma2))%0A%20%20%20%20%20%20%20%20u_right%2C%20u_left%20%3D%20u%5B%3A%2C%20%3A%2C%20-1%5D%2C%20u%5B%3A%2C%20%3A%2C%200%5D%0A%20%20%20%20%20%20%20%20%23%20Exterior%20values%20come%20from%20the%20neighbouring%20elements%20(periodic%20domain).%0A%20%20%20%20%20%20%20%20corr_right%20%3D%20correction(u_right%2C%20np.roll(u_left%2C%20-1%2C%20axis%3D1)%2C%20%2B1.0%2C%20gamma2)%0A%20%20%20%20%20%20%20%20corr_left%20%3D%20correction(u_left%2C%20np.roll(u_right%2C%201%2C%20axis%3D1)%2C%20-1.0%2C%20gamma2)%0A%20%20%20%20%20%20%20%20%23%20Lift%20the%20corrections%3A%20M%5E%7B-1%7D%20%3D%201%20%2F%20(w%20J)%2C%20and%20only%20the%20endpoint%20nodes%20are%20touched.%0A%20%20%20%20%20%20%20%20dt_u%5B%3A%2C%20%3A%2C%20-1%5D%20-%3D%20corr_right%20%2F%20(grid.w%5B-1%5D%20*%20grid.J)%0A%20%20%20%20%20%20%20%20dt_u%5B%3A%2C%20%3A%2C%200%5D%20-%3D%20corr_left%20%2F%20(grid.w%5B0%5D%20*%20grid.J)%0A%20%20%20%20%20%20%20%20return%20dt_u.ravel()%0A%0A%20%20%20%20return%20(dg_rhs%2C)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(%0A%20%20%20%20boundary_correction_central%2C%0A%20%20%20%20boundary_correction_rusanov%2C%0A%20%20%20%20boundary_correction_upwind%2C%0A%20%20%20%20check%2C%0A%20%20%20%20dg_rhs%2C%0A%20%20%20%20flux%2C%0A%20%20%20%20incomplete%2C%0A%20%20%20%20make_grid%2C%0A%20%20%20%20np%2C%0A%20%20%20%20volume_terms%2C%0A)%3A%0A%20%20%20%20_g%20%3D%20make_grid(4%2C%208)%0A%20%20%20%20_x%20%3D%20_g.x%0A%20%20%20%20_u%20%3D%20np.stack(%5Bnp.sin(_x)%2C%20np.cos(_x)%2C%20np.cos(_x)%5D)%0A%20%20%20%20_exact_dt%20%3D%20np.stack(%5B-np.cos(_x)%2C%20np.sin(_x)%2C%20np.sin(_x)%5D)%0A%20%20%20%20_rng%20%3D%20np.random.default_rng(1)%0A%20%20%20%20_ui%2C%20_ue%20%3D%20_rng.normal(size%3D(2%2C%203%2C%205))%0A%20%20%20%20_g2%20%3D%200.7%0A%0A%20%20%20%20vol_done%20%3D%20not%20incomplete(volume_terms(_u%2C%20_g%2C%20_g2))%0A%20%20%20%20upwind_done%20%3D%20not%20incomplete(boundary_correction_upwind(_ui%2C%20_ue%2C%201.0%2C%20_g2))%0A%20%20%20%20rusanov_done%20%3D%20not%20incomplete(boundary_correction_rusanov(_ui%2C%20_ue%2C%201.0%2C%20_g2))%0A%0A%0A%20%20%20%20def%20_conservation(corr)%3A%0A%20%20%20%20%20%20%20%20%23%20G(int%2C%20ext%3B%20n)%20%3D%20-G(ext%2C%20int%3B%20-n)%20with%20G%20%3D%20correction%20%2B%20n%20F_int%0A%20%20%20%20%20%20%20%20_G1%20%3D%20corr(_ui%2C%20_ue%2C%201.0%2C%20_g2)%20%2B%20flux(_ui%2C%20_g2)%0A%20%20%20%20%20%20%20%20_G2%20%3D%20corr(_ue%2C%20_ui%2C%20-1.0%2C%20_g2)%20-%20flux(_ue%2C%20_g2)%0A%20%20%20%20%20%20%20%20return%20np.allclose(_G1%2C%20-_G2)%0A%0A%0A%20%20%20%20_checks%20%3D%20%5B%5D%0A%20%20%20%20if%20vol_done%3A%0A%20%20%20%20%20%20%20%20_v%20%3D%20volume_terms(_u%2C%20_g%2C%201.0)%0A%20%20%20%20%20%20%20%20_checks%20%2B%3D%20%5B%0A%20%20%20%20%20%20%20%20%20%20%20%20(np.shape(_v)%20%3D%3D%20_u.shape%2C%20f%22volume_terms%20should%20return%20shape%20%7B_u.shape%7D%2C%20got%20%7Bnp.shape(_v)%7D.%22)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20(%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20np.shape(_v)%20%3D%3D%20_u.shape%20and%20np.max(np.abs(_v%20-%20_exact_dt))%20%3C%201e-6%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22For%20smooth%20continuous%20data%2C%20volume_terms%20should%20equal%20du%2Fdt%20to%20spectral%20accuracy.%20Check%20signs%20and%201%2FJ.%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20)%2C%0A%20%20%20%20%20%20%20%20%5D%0A%20%20%20%20if%20upwind_done%3A%0A%20%20%20%20%20%20%20%20_c%20%3D%20boundary_correction_upwind%0A%20%20%20%20%20%20%20%20_wplus_jump%20%3D%20np.stack(%5Bnp.zeros(5)%2C%20np.ones(5)%2C%20np.ones(5)%5D)%20%20%23%20jump%20only%20in%20w%5E%2B%20(n%20%3D%20%2B1)%0A%20%20%20%20%20%20%20%20_checks%20%2B%3D%20%5B%0A%20%20%20%20%20%20%20%20%20%20%20%20(np.allclose(_c(_ui%2C%20_ui%2C%201.0%2C%20_g2)%2C%200)%2C%20%22Upwind%3A%20the%20correction%20must%20vanish%20for%20continuous%20data%20(consistency).%22)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20(_conservation(_c)%2C%20%22Upwind%3A%20the%20conservation%20property%20G(int%2C%20ext%3B%20n)%20%3D%20-G(ext%2C%20int%3B%20-n)%20fails.%22)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20(np.allclose(_c(_ui%2C%20_ui%20%2B%20_wplus_jump%2C%201.0%2C%20_g2)%2C%200)%2C%20%22Upwind%3A%20a%20jump%20only%20in%20the%20outgoing%20field%20w%5E%2B%20must%20not%20produce%20a%20correction.%22)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20(%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20np.allclose(_c(_ui%2C%20_ui%20%2B%20np.stack(%5Bnp.zeros(5)%2C%20np.ones(5)%2C%20-np.ones(5)%5D)%2C%201.0%2C%200.0)%5B1%3A%5D%2C%20%5B%5B-1.0%5D%20*%205%2C%20%5B1.0%5D%20*%205%5D)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22Upwind%3A%20a%20jump%20of%202%20in%20w%5E-%20(n%20%3D%20%2B1)%20should%20give%20a%20correction%20(0%2C%20-1%2C%20%2B1).%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20)%2C%0A%20%20%20%20%20%20%20%20%5D%0A%20%20%20%20if%20rusanov_done%3A%0A%20%20%20%20%20%20%20%20_c%20%3D%20boundary_correction_rusanov%0A%20%20%20%20%20%20%20%20_jump_psi%20%3D%20np.stack(%5Bnp.ones(5)%2C%20np.zeros(5)%2C%20np.zeros(5)%5D)%0A%20%20%20%20%20%20%20%20_checks%20%2B%3D%20%5B%0A%20%20%20%20%20%20%20%20%20%20%20%20(np.allclose(_c(_ui%2C%20_ui%2C%20-1.0%2C%20_g2)%2C%200)%2C%20%22Rusanov%3A%20the%20correction%20must%20vanish%20for%20continuous%20data.%22)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20(_conservation(_c)%2C%20%22Rusanov%3A%20the%20conservation%20property%20fails.%22)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20(%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20np.allclose(_c(_ui%2C%20_ui%20%2B%20_jump_psi%2C%201.0%2C%20_g2)%2C%20np.stack(%5B-0.5%20*%20np.ones(5)%2C%20np.zeros(5)%2C%20-0.5%20*%20_g2%20*%20np.ones(5)%5D))%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22Rusanov%3A%20a%20unit%20jump%20in%20Psi%20(n%20%3D%20%2B1)%20should%20give%20(-1%2F2%2C%200%2C%20-gamma2%2F2).%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20)%2C%0A%20%20%20%20%20%20%20%20%5D%0A%20%20%20%20if%20vol_done%20and%20upwind_done%3A%0A%20%20%20%20%20%20%20%20for%20_corr%20in%20%5Bboundary_correction_upwind%2C%20boundary_correction_central%5D%20%2B%20(%5Bboundary_correction_rusanov%5D%20if%20rusanov_done%20else%20%5B%5D)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20_r%20%3D%20dg_rhs(0.0%2C%20_u.ravel()%2C%20_g%2C%201.0%2C%20_corr).reshape(_u.shape)%0A%20%20%20%20%20%20%20%20%20%20%20%20_checks.append((np.max(np.abs(_r%20-%20_exact_dt))%20%3C%201e-6%2C%20f%22The%20full%20DG%20right-hand%20side%20with%20%7B_corr.__name__%7D%20is%20not%20accurate.%22))%0A%20%20%20%20ex6_done%20%3D%20vol_done%20and%20upwind_done%0A%20%20%20%20check(%22Exercise%206%3A%20DG%20volume%20terms%20and%20boundary%20corrections%22%2C%20_checks)%0A%20%20%20%20return%20ex6_done%2C%20rusanov_done%0A%0A%0A%40app.cell%0Adef%20_(dg_rhs%2C%20make_grid%2C%20np%2C%20solve_ivp)%3A%0A%20%20%20%20def%20exact_solution(t%2C%20x)%3A%0A%20%20%20%20%20%20%20%20%22%22%22Exact%20(psi%2C%20pi%2C%20phi)%20for%20the%20right-moving%20wave%20Psi%20%3D%20sin(x%20-%20t)%2C%20stacked.%22%22%22%0A%20%20%20%20%20%20%20%20return%20np.stack(%5Bnp.sin(x%20-%20t)%2C%20np.cos(x%20-%20t)%2C%20np.cos(x%20-%20t)%5D)%0A%0A%0A%20%20%20%20def%20evolve_dg(K%2C%20N%2C%20t_final%2C%20correction%2C%20gamma2%3D1.0%2C%20rtol%3D1e-12%2C%20atol%3D1e-12%2C%20data%3DNone)%3A%0A%20%20%20%20%20%20%20%20%22%22%22Evolve%20the%20DG%20scalar%20wave%3B%20returns%20(sol%2C%20grid).%20%60data(x)%60%20optionally%20sets%20the%20initial%20data.%22%22%22%0A%20%20%20%20%20%20%20%20grid%20%3D%20make_grid(K%2C%20N)%0A%20%20%20%20%20%20%20%20u0%20%3D%20exact_solution(0.0%2C%20grid.x)%20if%20data%20is%20None%20else%20data(grid.x)%0A%20%20%20%20%20%20%20%20sol%20%3D%20solve_ivp(%0A%20%20%20%20%20%20%20%20%20%20%20%20dg_rhs%2C%20(0.0%2C%20t_final)%2C%20u0.ravel()%2C%20method%3D%22RK45%22%2C%20rtol%3Drtol%2C%20atol%3Datol%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20dense_output%3DTrue%2C%20args%3D(grid%2C%20gamma2%2C%20correction)%2C%0A%20%20%20%20%20%20%20%20)%0A%20%20%20%20%20%20%20%20return%20sol%2C%20grid%0A%0A%20%20%20%20return%20evolve_dg%2C%20exact_solution%0A%0A%0A%40app.cell%0Adef%20_(%0A%20%20%20%20T%2C%0A%20%20%20%20boundary_correction_central%2C%0A%20%20%20%20boundary_correction_rusanov%2C%0A%20%20%20%20boundary_correction_upwind%2C%0A%20%20%20%20mo%2C%0A%20%20%20%20rusanov_done%2C%0A)%3A%0A%20%20%20%20corrections%20%3D%20%7B%22upwind%22%3A%20boundary_correction_upwind%2C%20%22central%22%3A%20boundary_correction_central%7D%0A%20%20%20%20if%20rusanov_done%3A%0A%20%20%20%20%20%20%20%20corrections%5B%22Rusanov%22%5D%20%3D%20boundary_correction_rusanov%0A%20%20%20%20K_slider%20%3D%20mo.ui.slider(1%2C%2016%2C%20value%3D4%2C%20label%3D%22elements%20%24K%24%22%2C%20show_value%3DTrue)%0A%20%20%20%20Ndeg_slider%20%3D%20mo.ui.slider(1%2C%2012%2C%20value%3D3%2C%20label%3D%22degree%20%24N%24%22%2C%20show_value%3DTrue)%0A%20%20%20%20flux_dropdown%20%3D%20mo.ui.dropdown(options%3Dlist(corrections)%2C%20value%3D%22upwind%22%2C%20label%3D%22boundary%20correction%22)%0A%20%20%20%20dg_gamma2_slider%20%3D%20mo.ui.slider(0.0%2C%203.0%2C%20step%3D0.1%2C%20value%3D1.0%2C%20label%3Dr%22%24%5Cgamma_2%24%22%2C%20show_value%3DTrue)%0A%20%20%20%20dg_t_slider%20%3D%20mo.ui.slider(0.0%2C%202%20*%20T%2C%20step%3D2%20*%20T%20%2F%20400%2C%20value%3D0.25%20*%20T%2C%20label%3D%22output%20time%20%24t%24%22%2C%20show_value%3DTrue)%0A%20%20%20%20mo.vstack(%0A%20%20%20%20%20%20%20%20%5B%0A%20%20%20%20%20%20%20%20%20%20%20%20mo.hstack(%5BK_slider%2C%20Ndeg_slider%2C%20flux_dropdown%5D%2C%20justify%3D%22start%22%2C%20gap%3D2)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20mo.hstack(%5Bdg_gamma2_slider%2C%20dg_t_slider%5D%2C%20justify%3D%22start%22%2C%20gap%3D2)%2C%0A%20%20%20%20%20%20%20%20%5D%0A%20%20%20%20)%0A%20%20%20%20return%20(%0A%20%20%20%20%20%20%20%20K_slider%2C%0A%20%20%20%20%20%20%20%20Ndeg_slider%2C%0A%20%20%20%20%20%20%20%20corrections%2C%0A%20%20%20%20%20%20%20%20dg_gamma2_slider%2C%0A%20%20%20%20%20%20%20%20dg_t_slider%2C%0A%20%20%20%20%20%20%20%20flux_dropdown%2C%0A%20%20%20%20)%0A%0A%0A%40app.cell%0Adef%20_(%0A%20%20%20%20K_slider%2C%0A%20%20%20%20Ndeg_slider%2C%0A%20%20%20%20T%2C%0A%20%20%20%20corrections%2C%0A%20%20%20%20dg_gamma2_slider%2C%0A%20%20%20%20evolve_dg%2C%0A%20%20%20%20ex6_done%2C%0A%20%20%20%20flux_dropdown%2C%0A%20%20%20%20mo%2C%0A)%3A%0A%20%20%20%20mo.stop(not%20ex6_done%2C%20mo.callout(mo.md(%22%E2%AC%86%20Complete%20**Exercise%206**%20(at%20least%206a%20and%206b)%20to%20run%20the%20DG%20evolution.%22)%2C%20kind%3D%22warn%22))%0A%20%20%20%20dg_sol%2C%20dg_grid%20%3D%20evolve_dg(%0A%20%20%20%20%20%20%20%20K_slider.value%2C%20Ndeg_slider.value%2C%202%20*%20T%2C%20corrections%5Bflux_dropdown.value%5D%2C%0A%20%20%20%20%20%20%20%20gamma2%3Ddg_gamma2_slider.value%2C%20rtol%3D1e-8%2C%20atol%3D1e-8%2C%0A%20%20%20%20)%0A%20%20%20%20mo.md(%0A%20%20%20%20%20%20%20%20f%22%24K%3D%7Bdg_grid.K%7D%24%20elements%20%C3%97%20%24N%2B1%3D%7Bdg_grid.N%20%2B%201%7D%24%20nodes%20%24%3D%20%7Bdg_grid.K%20*%20(dg_grid.N%20%2B%201)%7D%24%20points%20per%20variable%3B%20%22%0A%20%20%20%20%20%20%20%20f%22**%7Bdg_sol.t.size%20-%201%7D%20accepted%20steps**%2C%20%7Bdg_sol.nfev%7D%20RHS%20evaluations%20(%60rtol%20%3D%20atol%20%3D%201e-8%60).%22%0A%20%20%20%20)%0A%20%20%20%20return%20dg_grid%2C%20dg_sol%0A%0A%0A%40app.cell%0Adef%20_(T%2C%20dg_grid%2C%20dg_sol%2C%20dg_t_slider%2C%20exact_solution%2C%20interp_matrix%2C%20np%2C%20plt)%3A%0A%20%20%20%20_t%20%3D%20dg_t_slider.value%0A%20%20%20%20_u%20%3D%20dg_sol.sol(_t).reshape(3%2C%20dg_grid.K%2C%20dg_grid.N%20%2B%201)%0A%20%20%20%20_xi_f%20%3D%20np.linspace(-1%2C%201%2C%2040)%0A%20%20%20%20_I%20%3D%20interp_matrix(dg_grid.xi%2C%20_xi_f)%0A%20%20%20%20_dx%20%3D%20T%20%2F%20dg_grid.K%0A%20%20%20%20_names%20%3D%20%5Br%22%24%5CPsi%24%22%2C%20r%22%24%5CPi%24%22%2C%20r%22%24%5CPhi%24%22%5D%0A%0A%20%20%20%20_fig%2C%20(_ax0%2C%20_ax1)%20%3D%20plt.subplots(2%2C%201%2C%20figsize%3D(7.5%2C%206)%2C%20sharex%3DTrue)%0A%20%20%20%20for%20_k%20in%20range(dg_grid.K)%3A%0A%20%20%20%20%20%20%20%20_xf%20%3D%20_k%20*%20_dx%20%2B%20(1%20%2B%20_xi_f)%20*%20_dx%20%2F%202%0A%20%20%20%20%20%20%20%20_ex%20%3D%20exact_solution(_t%2C%20_xf)%0A%20%20%20%20%20%20%20%20for%20_i%20in%20range(3)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20_uf%20%3D%20_I%20%40%20_u%5B_i%2C%20_k%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20_ax0.plot(_xf%2C%20_uf%2C%20%22-%22%2C%20color%3Df%22C%7B_i%7D%22%2C%20lw%3D1.5%2C%20label%3D_names%5B_i%5D%20if%20_k%20%3D%3D%200%20else%20None)%0A%20%20%20%20%20%20%20%20%20%20%20%20_ax0.plot(dg_grid.x%5B_k%5D%2C%20_u%5B_i%2C%20_k%5D%2C%20%22o%22%2C%20color%3Df%22C%7B_i%7D%22%2C%20ms%3D3)%0A%20%20%20%20%20%20%20%20%20%20%20%20_ax1.semilogy(_xf%2C%20np.abs(_uf%20-%20_ex%5B_i%5D)%20%2B%201e-17%2C%20%22-%22%2C%20color%3Df%22C%7B_i%7D%22%2C%20lw%3D1)%0A%20%20%20%20%20%20%20%20_ax0.axvline(_k%20*%20_dx%2C%20color%3D%22gray%22%2C%20lw%3D0.5)%0A%20%20%20%20%20%20%20%20_ax1.axvline(_k%20*%20_dx%2C%20color%3D%22gray%22%2C%20lw%3D0.5)%0A%20%20%20%20_ax0.set_title(f%22t%20%3D%20%7B_t%20%2F%20T%3A.3f%7D%20T%20%20%20(lines%3A%20DG%20polynomial%20in%20each%20element%2C%20dots%3A%20nodes)%22)%0A%20%20%20%20_ax0.legend(loc%3D%22upper%20right%22)%0A%20%20%20%20_ax1.set_ylabel(%22%7CDG%20%E2%88%92%20exact%7C%22)%0A%20%20%20%20_ax1.set_xlabel(%22%24x%24%22)%0A%20%20%20%20_fig.tight_layout()%0A%20%20%20%20_fig%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20Things%20to%20try%20with%20the%20controls%3A%0A%0A%20%20%20%20*%20Fix%20%24K%24%20and%20raise%20%24N%24%20(knobs%202%20%2B%203)%3A%20the%20error%20drops%20by%20orders%20of%20magnitude%20per%20step.%0A%20%20%20%20*%20Fix%20%24N%24%20and%20raise%20%24K%24%20(knob%201)%3A%20the%20error%20drops%20algebraically.%0A%20%20%20%20*%20Compare%20*upwind*%20and%20*central*%20at%20low%20resolution%20(e.g.%20%24K%3D4%24%2C%20%24N%3D2%24)%3A%20watch%20the%20jumps%20between%20elements%20in%20the%20error%0A%20%20%20%20%20%20panel%2C%20and%20the%20number%20of%20time%20steps.%0A%0A%20%20%20%20%23%23%23%20Demo%20%E2%80%94%20dissipation%20from%20the%20boundary%20correction%0A%0A%20%20%20%20Evolve%20an%20**under-resolved**%20narrow%20Gaussian%20pulse%20(%24K%3D8%24%2C%20%24N%3D3%24%2C%20%24%5Cgamma_2%3D0%24)%20and%20monitor%20the%20energy%0A%20%20%20%20%24E%3D%5Cfrac12%5Cint(%5CPi%5E2%2B%5CPhi%5E2)%5C%2Cdx%24%2C%20which%20is%20exactly%20conserved%20by%20the%20continuum%20equations.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(T%2C%20corrections%2C%20evolve_dg%2C%20ex6_done%2C%20mo%2C%20np%2C%20plt)%3A%0A%20%20%20%20mo.stop(not%20ex6_done%2C%20mo.callout(mo.md(%22%E2%AC%86%20Complete%20**Exercise%206**%20first.%22)%2C%20kind%3D%22warn%22))%0A%0A%0A%20%20%20%20def%20_pulse(x%2C%20width%3D0.3)%3A%0A%20%20%20%20%20%20%20%20_psi%20%3D%20np.exp(-(((x%20-%20np.pi)%20%2F%20width)%20**%202))%0A%20%20%20%20%20%20%20%20_dpsi%20%3D%20-2%20*%20(x%20-%20np.pi)%20%2F%20width**2%20*%20_psi%0A%20%20%20%20%20%20%20%20return%20np.stack(%5B_psi%2C%20_dpsi%2C%20_dpsi%5D)%20%20%23%20right-moving%3A%20Pi%20%3D%20Phi%20%3D%20dPsi%2Fdx%0A%0A%0A%20%20%20%20_ts%20%3D%20np.linspace(0%2C%203%20*%20T%2C%20300)%0A%20%20%20%20_fig%2C%20_ax%20%3D%20plt.subplots(figsize%3D(6.5%2C%204))%0A%20%20%20%20_rows%20%3D%20%5B%5D%0A%20%20%20%20for%20_name%2C%20_corr%20in%20corrections.items()%3A%0A%20%20%20%20%20%20%20%20_sol%2C%20_grid%20%3D%20evolve_dg(8%2C%203%2C%203%20*%20T%2C%20_corr%2C%20gamma2%3D0.0%2C%20rtol%3D1e-10%2C%20atol%3D1e-10%2C%20data%3D_pulse)%0A%20%20%20%20%20%20%20%20_E%20%3D%20%5B%0A%20%20%20%20%20%20%20%20%20%20%20%200.5%20*%20np.sum(_grid.w%20*%20_grid.J%20*%20_sol.sol(_t).reshape(3%2C%20_grid.K%2C%20_grid.N%20%2B%201)%5B1%3A%5D%20**%202)%0A%20%20%20%20%20%20%20%20%20%20%20%20for%20_t%20in%20_ts%0A%20%20%20%20%20%20%20%20%5D%0A%20%20%20%20%20%20%20%20_style%20%3D%20%7B%22upwind%22%3A%20dict(lw%3D4%2C%20alpha%3D0.5)%2C%20%22central%22%3A%20dict(lw%3D2)%2C%20%22Rusanov%22%3A%20dict(lw%3D1.5%2C%20ls%3D%22--%22%2C%20color%3D%22k%22)%7D%0A%20%20%20%20%20%20%20%20_ax.plot(_ts%20%2F%20T%2C%20np.array(_E)%20%2F%20_E%5B0%5D%2C%20label%3D_name%2C%20**_style%5B_name%5D)%0A%20%20%20%20%20%20%20%20_rows.append(f%22%7C%20%7B_name%7D%20%7C%20%7B_sol.t.size%20-%201%7D%20%7C%22)%0A%20%20%20%20_ax.set_xlabel(%22%24t%2FT%24%22)%0A%20%20%20%20_ax.set_ylabel(%22%24E(t)%2FE(0)%24%22)%0A%20%20%20%20_ax.set_title(%22Energy%20of%20an%20under-resolved%20pulse%20(K%20%3D%208%2C%20N%20%3D%203)%22)%0A%20%20%20%20_ax.legend()%0A%20%20%20%20_fig.tight_layout()%0A%20%20%20%20mo.hstack(%0A%20%20%20%20%20%20%20%20%5B_fig%2C%20mo.md(%22%7C%20correction%20%7C%20accepted%20steps%20%7C%5Cn%7C---%7C---%7C%5Cn%22%20%2B%20%22%5Cn%22.join(_rows))%5D%2C%0A%20%20%20%20%20%20%20%20align%3D%22center%22%2C%0A%20%20%20%20)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20*%20The%20**central**%20flux%20conserves%20the%20discrete%20energy%20to%20integration%20accuracy%3A%20the%20under-resolved%2C%20wrong%20parts%20of%20the%20solution%0A%20%20%20%20%20%20keep%20sloshing%20around%20forever.%0A%20%20%20%20*%20**Upwind**%20(and%20**Rusanov**)%20remove%20energy%20through%20the%20jumps%2C%20damping%20precisely%20the%20unresolved%20content.%20For%20this%20system%2C%20where%0A%20%20%20%20%20%20every%20nonzero%20speed%20equals%20%24%5Clambda_%7B%5Cmax%7D%3D1%24%2C%20the%20two%20coincide%20in%20the%20%24%5CPi%24%20and%20%24%5CPhi%24%20equations%20and%20differ%20only%20in%0A%20%20%20%20%20%20%24%5CPsi%24%2C%20where%20Rusanov%20adds%20extra%20dissipation.%0A%20%20%20%20*%20The%20central%20flux%20also%20needs%20**more%20time%20steps**.%20Its%20semi-discrete%20eigenvalues%20lie%20on%20the%20imaginary%20axis%2C%20exactly%20where%0A%20%20%20%20%20%20Dormand%E2%80%93Prince%20is%20weakest%20(Lecture%201)%2C%20whereas%20the%20upwind%20eigenvalues%20have%20negative%20real%20parts.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20---%0A%20%20%20%20%23%23%20Part%206%20%E2%80%94%20Convergence%20of%20DG%3A%20%24p%24-%20and%20%24h%24-refinement%0A%0A%20%20%20%20We%20measure%20the%20%24L_2%24%20error%20using%20the%20LGL%20quadrature%20on%20each%20element%2C%0A%20%20%20%20%24%5C%7Ce%5C%7C_2%5E2%20%5Capprox%20%5Csum_k%5Csum_i%20J%5C%2Cw_i%5C%2Ce_%7Bk%2Ci%7D%5E2%24%20(provided%20below%20as%20%60dg_l2_errors%60).%20As%20in%20Lecture%201%20we%20check%20all%20three%0A%20%20%20%20variables%20at%20%24t%3DT%24%20and%20%24t%3D1.0927%5C%2CT%24%2C%20with%20%60rtol%20%3D%20atol%20%3D%201e-12%60.%0A%0A%20%20%20%20%23%23%23%20Exercise%207%20%E2%80%94%20convergence%20test%0A%0A%20%20%20%20Using%20%60evolve_dg%60%20with%20the%20**upwind**%20correction%2C%20fill%20in%0A%0A%20%20%20%201.%20**%24p%24-convergence**%20(knobs%202%20%2B%203)%3A%20%60p_errors%5Blabel%5D%60%2C%20shape%20%60(len(p_Ns)%2C%203)%60%2C%20for%20%24K%3D4%24%20and%20%24N%24%20in%20%60p_Ns%60%3B%0A%20%20%20%202.%20**%24h%24-convergence**%20(knob%201)%3A%20%60h_errors%5BN%5D%5Blabel%5D%60%2C%20shape%20%60(len(h_Ks)%2C%203)%60%2C%20for%20each%20%24N%24%20in%20%60h_Ns%60%20and%20%24K%24%20in%20%60h_Ks%60.%0A%0A%20%20%20%20Expect%20a%20straight%20line%20on%20a%20**semilog**%20plot%20for%20%24p%24-convergence%2C%20and%20slope%20%24-(N%2B1)%24%20on%20a%20**log-log**%20plot%20for%0A%20%20%20%20%24h%24-convergence.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(exact_solution%2C%20np)%3A%0A%20%20%20%20def%20dg_l2_errors(sol%2C%20grid%2C%20t)%3A%0A%20%20%20%20%20%20%20%20%22%22%22L2%20errors%20of%20(psi%2C%20pi%2C%20phi)%20at%20time%20t%2C%20via%20LGL%20quadrature.%20Returns%20shape%20(3%2C).%22%22%22%0A%20%20%20%20%20%20%20%20_e%20%3D%20sol.sol(t).reshape(3%2C%20grid.K%2C%20grid.N%20%2B%201)%20-%20exact_solution(t%2C%20grid.x)%0A%20%20%20%20%20%20%20%20return%20np.sqrt(np.sum(grid.w%20*%20grid.J%20*%20_e**2%2C%20axis%3D(1%2C%202)))%0A%0A%20%20%20%20return%20(dg_l2_errors%2C)%0A%0A%0A%40app.cell%0Adef%20_(%0A%20%20%20%20T%2C%0A%20%20%20%20boundary_correction_upwind%2C%0A%20%20%20%20dg_l2_errors%2C%0A%20%20%20%20evolve_dg%2C%0A%20%20%20%20ex6_done%2C%0A%20%20%20%20mo%2C%0A%20%20%20%20np%2C%0A)%3A%0A%20%20%20%20mo.stop(not%20ex6_done%2C%20mo.callout(mo.md(%22%E2%AC%86%20Complete%20**Exercise%206**%20first.%22)%2C%20kind%3D%22warn%22))%0A%20%20%20%20dg_t_checks%20%3D%20%7B%22t%20%3D%20T%22%3A%20T%2C%20%22t%20%3D%201.0927%20T%22%3A%201.0927%20*%20T%7D%0A%20%20%20%20p_K%20%3D%204%0A%20%20%20%20p_Ns%20%3D%20np.arange(2%2C%2013)%0A%20%20%20%20h_Ns%20%3D%20%5B2%2C%203%2C%204%5D%0A%20%20%20%20h_Ks%20%3D%20np.array(%5B2%2C%204%2C%208%2C%2016%2C%2032%5D)%0A%20%20%20%20p_errors%20%3D%20%7B%7D%20%20%23%20p_errors%5Blabel%5D%3A%20shape%20(len(p_Ns)%2C%203)%0A%20%20%20%20h_errors%20%3D%20%7B%7D%20%20%23%20h_errors%5BN%5D%5Blabel%5D%3A%20shape%20(len(h_Ks)%2C%203)%0A%20%20%20%20%23%20Tip%3A%20evolve%20each%20(K%2C%20N)%20once%20to%201.1%20T%20and%20evaluate%20dg_l2_errors%20at%20both%20check%20times.%0A%20%20%20%20%23%23%23%20BEGIN%20SOLUTION%0A%20%20%20%20_p_runs%20%3D%20%5Bevolve_dg(p_K%2C%20_N%2C%201.1%20*%20T%2C%20boundary_correction_upwind)%20for%20_N%20in%20p_Ns%5D%0A%20%20%20%20for%20_label%2C%20_t%20in%20dg_t_checks.items()%3A%0A%20%20%20%20%20%20%20%20p_errors%5B_label%5D%20%3D%20np.array(%5Bdg_l2_errors(_s%2C%20_g%2C%20_t)%20for%20_s%2C%20_g%20in%20_p_runs%5D)%0A%20%20%20%20for%20_N%20in%20h_Ns%3A%0A%20%20%20%20%20%20%20%20_h_runs%20%3D%20%5Bevolve_dg(_K%2C%20_N%2C%201.1%20*%20T%2C%20boundary_correction_upwind)%20for%20_K%20in%20h_Ks%5D%0A%20%20%20%20%20%20%20%20h_errors%5B_N%5D%20%3D%20%7B%0A%20%20%20%20%20%20%20%20%20%20%20%20_label%3A%20np.array(%5Bdg_l2_errors(_s%2C%20_g%2C%20_t)%20for%20_s%2C%20_g%20in%20_h_runs%5D)%0A%20%20%20%20%20%20%20%20%20%20%20%20for%20_label%2C%20_t%20in%20dg_t_checks.items()%0A%20%20%20%20%20%20%20%20%7D%0A%20%20%20%20%23%23%23%20END%20SOLUTION%0A%20%20%20%20return%20dg_t_checks%2C%20h_Ks%2C%20h_Ns%2C%20h_errors%2C%20p_K%2C%20p_Ns%2C%20p_errors%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(check%2C%20dg_t_checks%2C%20h_Ks%2C%20h_Ns%2C%20h_errors%2C%20np%2C%20p_Ns%2C%20p_errors)%3A%0A%20%20%20%20ex7_done%20%3D%20len(p_errors)%20%3D%3D%20len(dg_t_checks)%20and%20len(h_errors)%20%3D%3D%20len(h_Ns)%0A%20%20%20%20_checks%20%3D%20%5B%5D%0A%20%20%20%20if%20len(p_errors)%20%3D%3D%20len(dg_t_checks)%3A%0A%20%20%20%20%20%20%20%20for%20_label%20in%20dg_t_checks%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20_e%20%3D%20np.asarray(p_errors%5B_label%5D)%0A%20%20%20%20%20%20%20%20%20%20%20%20_i10%20%3D%20list(p_Ns).index(10)%0A%20%20%20%20%20%20%20%20%20%20%20%20_checks.append(%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20(_e.shape%20%3D%3D%20(len(p_Ns)%2C%203)%20and%20np.all(_e%5B_i10%5D%20%3C%201e-9)%2C%20f%22p-convergence%20at%20%7B_label%7D%3A%20errors%20at%20N%20%3D%2010%20should%20be%20%3C%201e-9.%22)%0A%20%20%20%20%20%20%20%20%20%20%20%20)%0A%20%20%20%20if%20len(h_errors)%20%3D%3D%20len(h_Ns)%3A%0A%20%20%20%20%20%20%20%20for%20_N%20in%20h_Ns%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20for%20_label%20in%20dg_t_checks%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20_e%20%3D%20np.asarray(h_errors%5B_N%5D%5B_label%5D)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20_ok%20%3D%20_e.shape%20%3D%3D%20(len(h_Ks)%2C%203)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20_order%20%3D%20np.log2(_e%5B-2%5D%20%2F%20_e%5B-1%5D)%20if%20_ok%20else%20np.nan%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20_checks.append(%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20(_ok%20and%20np.all(np.abs(_order%20-%20(_N%20%2B%201))%20%3C%200.35)%2C%20f%22h-convergence%20N%20%3D%20%7B_N%7D%20at%20%7B_label%7D%3A%20finest%20orders%20%7Bnp.round(_order%2C%202)%7D%2C%20expected%20about%20%7B_N%20%2B%201%7D.%22)%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20)%0A%20%20%20%20check(%22Exercise%207%3A%20DG%20convergence%22%2C%20_checks)%0A%20%20%20%20return%20(ex7_done%2C)%0A%0A%0A%40app.cell%0Adef%20_(%0A%20%20%20%20dg_t_checks%2C%0A%20%20%20%20ex7_done%2C%0A%20%20%20%20h_Ks%2C%0A%20%20%20%20h_Ns%2C%0A%20%20%20%20h_errors%2C%0A%20%20%20%20mo%2C%0A%20%20%20%20np%2C%0A%20%20%20%20p_K%2C%0A%20%20%20%20p_Ns%2C%0A%20%20%20%20p_errors%2C%0A%20%20%20%20plt%2C%0A)%3A%0A%20%20%20%20mo.stop(not%20ex7_done%2C%20mo.callout(mo.md(%22%E2%AC%86%20Complete%20**Exercise%207**%20to%20see%20the%20convergence%20plots.%22)%2C%20kind%3D%22warn%22))%0A%20%20%20%20_names%20%3D%20%5Br%22%24%5CPsi%24%22%2C%20r%22%24%5CPi%24%22%2C%20r%22%24%5CPhi%24%22%5D%0A%20%20%20%20_styles%20%3D%20%7B%22t%20%3D%20T%22%3A%20%22-%22%2C%20%22t%20%3D%201.0927%20T%22%3A%20%22%3A%22%7D%0A%20%20%20%20_fig%2C%20(_ax0%2C%20_ax1)%20%3D%20plt.subplots(1%2C%202%2C%20figsize%3D(11%2C%204.3))%0A%20%20%20%20for%20_label%20in%20dg_t_checks%3A%0A%20%20%20%20%20%20%20%20for%20_i%20in%20range(3)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20_ax0.semilogy(p_Ns%2C%20p_errors%5B_label%5D%5B%3A%2C%20_i%5D%2C%20%22o%22%20%2B%20_styles%5B_label%5D%2C%20color%3Df%22C%7B_i%7D%22%2C%20ms%3D4%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20label%3Df%22%7B_names%5B_i%5D%7D%2C%20%7B_label%7D%22)%0A%20%20%20%20_ax0.set_xlabel(%22polynomial%20degree%20%24N%24%22)%0A%20%20%20%20_ax0.set_ylabel(r%22%24L_2%24%20error%22)%0A%20%20%20%20_ax0.set_title(f%22p-convergence%20(K%20%3D%20%7Bp_K%7D)%3A%20exponential%22)%0A%20%20%20%20_ax0.legend(ncol%3D2%2C%20fontsize%3D8)%0A%0A%20%20%20%20_dxs%20%3D%202%20*%20np.pi%20%2F%20h_Ks%0A%20%20%20%20for%20_j%2C%20_N%20in%20enumerate(h_Ns)%3A%0A%20%20%20%20%20%20%20%20for%20_i%20in%20range(3)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20_ax1.loglog(_dxs%2C%20h_errors%5B_N%5D%5B%22t%20%3D%201.0927%20T%22%5D%5B%3A%2C%20_i%5D%2C%20%22o-%22%2C%20color%3Df%22C%7B_j%7D%22%2C%20ms%3D3%2C%20alpha%3D0.4%20%2B%200.2%20*%20_i%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20label%3Df%22N%20%3D%20%7B_N%7D%22%20if%20_i%20%3D%3D%200%20else%20None)%0A%20%20%20%20%20%20%20%20_e0%20%3D%20h_errors%5B_N%5D%5B%22t%20%3D%201.0927%20T%22%5D%5B-1%2C%200%5D%0A%20%20%20%20%20%20%20%20_ax1.loglog(_dxs%2C%20_e0%20*%20(_dxs%20%2F%20_dxs%5B-1%5D)%20**%20(_N%20%2B%201)%2C%20%22k--%22%2C%20lw%3D0.8)%0A%20%20%20%20_ax1.plot(%5B%5D%2C%20%5B%5D%2C%20%22k--%22%2C%20lw%3D0.8%2C%20label%3Dr%22%24%5Cpropto%5CDelta%20x%5E%7BN%2B1%7D%24%22)%0A%20%20%20%20_ax1.set_xticks(_dxs%2C%20%5Bf%22%7B_d%3A.2f%7D%22%20for%20_d%20in%20_dxs%5D)%0A%20%20%20%20_ax1.minorticks_off()%0A%20%20%20%20_ax1.set_xlabel(r%22element%20size%20%24%5CDelta%20x%24%22)%0A%20%20%20%20_ax1.set_title(%22h-convergence%20at%20t%20%3D%201.0927%20T%20(all%203%20variables)%22)%0A%20%20%20%20_ax1.legend(fontsize%3D8)%0A%20%20%20%20_fig.tight_layout()%0A%0A%20%20%20%20_rows%20%3D%20%22%5Cn%22.join(%0A%20%20%20%20%20%20%20%20f%22%7C%20%7B_N%7D%20%7C%20%7B_N%20%2B%201%7D%20%7C%20%22%0A%20%20%20%20%20%20%20%20%2B%20%22%20%7C%20%22.join(%0A%20%20%20%20%20%20%20%20%20%20%20%20%22%20%2F%20%22.join(f%22%7Bnp.log2(h_errors%5B_N%5D%5B_l%5D%5B-2%2C%20_i%5D%20%2F%20h_errors%5B_N%5D%5B_l%5D%5B-1%2C%20_i%5D)%3A.2f%7D%22%20for%20_l%20in%20dg_t_checks)%0A%20%20%20%20%20%20%20%20%20%20%20%20for%20_i%20in%20range(3)%0A%20%20%20%20%20%20%20%20)%0A%20%20%20%20%20%20%20%20%2B%20%22%20%7C%22%0A%20%20%20%20%20%20%20%20for%20_N%20in%20h_Ns%0A%20%20%20%20)%0A%20%20%20%20mo.vstack(%0A%20%20%20%20%20%20%20%20%5B%0A%20%20%20%20%20%20%20%20%20%20%20%20_fig%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20mo.md(%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22**Observed%20h-convergence%20orders**%20between%20the%20two%20finest%20grids%20(%24t%3DT%24%20%2F%20%24t%3D1.0927%5C%5C%2CT%24)%3A%5Cn%5Cn%22%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22%7C%20%24N%24%20%7C%20expected%20%7C%20%24%5C%5CPsi%24%20%7C%20%24%5C%5CPi%24%20%7C%20%24%5C%5CPhi%24%20%7C%5Cn%7C---%7C---%7C---%7C---%7C---%7C%5Cn%22%20%2B%20_rows%0A%20%20%20%20%20%20%20%20%20%20%20%20)%2C%0A%20%20%20%20%20%20%20%20%5D%0A%20%20%20%20)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20**Reading%20the%20plots.**%0A%0A%20%20%20%20*%20%24p%24-convergence%20is%20a%20straight%20line%20on%20the%20semilog%20plot%20until%20it%20reaches%20the%20**time-integration%20floor**%20set%20by%0A%20%20%20%20%20%20%60rtol%20%3D%20atol%20%3D%201e-12%60.%20With%20%24K%3D4%24%2C%20%24N%3D10%24%20(%2444%24%20points%20per%20variable)%20the%20error%20is%20%24%5Csim10%5E%7B-11%7D%24.%20Compare%20Lecture%201%2C%20where%0A%20%20%20%20%20%204th-order%20FD%20with%20%24256%24%20points%20only%20reached%20%24%5Csim10%5E%7B-7%7D%24.%0A%20%20%20%20*%20%24h%24-refinement%20at%20fixed%20%24N%24%20converges%20as%20%24%5CDelta%20x%5E%7BN%2B1%7D%24%20(with%20upwind%20fluxes)%2C%20like%20a%20high-order%20finite-difference%20method.%0A%20%20%20%20*%20In%20practice%20one%20chooses%20%24N%24%20as%20high%20as%20the%20smoothness%20of%20the%20solution%20warrants%20and%20uses%20elements%20(and%20%24h%24-refinement)%20to%20handle%0A%20%20%20%20%20%20geometry%20and%20to%20isolate%20non-smooth%20regions.%0A%0A%20%20%20%20---%0A%20%20%20%20%23%23%20Summary%0A%0A%20%20%20%20*%20Nodal%20DG%20approximates%20the%20**solution**%20by%20a%20polynomial%20in%20each%20element%20and%20differentiates%20that%20approximation%20**exactly**.%0A%20%20%20%20*%20**Three%20knobs%3A**%20element%20size%2C%20number%20of%20points%2C%20point%20location.%20Points%20%2B%20location%20(LGL)%20give%20exponential%20convergence%3B%0A%20%20%20%20%20%20equally%20spaced%20points%20fail%20(Runge).%0A%20%20%20%20*%20Nodal%20representation%3A%20nonlinear%20fluxes%20and%20sources%20are%20evaluated%20**pointwise**%2C%20which%20is%20cheap%2C%20simple%20and%20general.%0A%20%20%20%20*%20Elements%20couple%20only%20through%20**boundary%20corrections**%2C%20which%20must%20be%20consistent%20and%20conservative%20and%20should%20be%20dissipative%0A%20%20%20%20%20%20(upwind%2C%20Rusanov%2C%20HLL%2C%20%E2%80%A6).%20The%20dissipation%20acts%20on%20the%20jumps%2C%20i.e.%20only%20where%20the%20solution%20is%20under-resolved.%0A%20%20%20%20*%20With%20LGL%20nodes%20and%20a%20lumped%20(diagonal)%20mass%20matrix%20the%20scheme%20is%0A%20%20%20%20%20%20%24%5Cdot%20u%20%3D%20-J%5E%7B-1%7DDF%20%2B%20S%20-%20%5Ctext%7Blift%7D(G%20-%20nF%5E%7B%5Crm%20int%7D)%24%2C%20which%20is%20only%20a%20few%20lines%20of%20NumPy.%0A%0A%20%20%20%20**Where%20next%3F**%20Non-smooth%20solutions%20(shocks)%20break%20exponential%20convergence%20and%20produce%20Gibbs%20oscillations.%20That%20is%20where%0A%20%20%20%20limiters%2C%20finite-volume%2Ffinite-difference%20subcells%2C%20and%20the%20SSP%20time%20integrators%20mentioned%20in%20Lecture%201%20come%20in.%20Curved%0A%20%20%20%20elements%2C%20multiple%20dimensions%20(tensor%20products%20of%20the%201d%20operators)%2C%20and%20moving%20meshes%20follow%20the%20same%20template.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0Aif%20__name__%20%3D%3D%20%22__main__%22%3A%0A%20%20%20%20app.run()%0A