lecture1 fd

%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%20%20%20%20%22sympy%22%2C%0A%23%20%5D%0A%23%20mograder-assignment%20%3D%20%22lecture1_fd%22%0A%23%20mograder-cell-hashes%20%3D%20%22a2e046dd%2C1fa92fa9%2Cb1534735%2Ccdcece45%2C199a8b35%2C7c08c6bb%2C02374ea0%2C87dc323e%2Cc0a44626%2C2ed6ccb9%2C51c70226%2C1b9560e2%2Ce240144f%2Cbe943441%2C66a05548%2C27309b7a%2Cd0f29b09%2Cc774109a%2C76a50567%2C5753641f%2Cee2388cb%2C4bcc56c5%2Ca59286ec%2C1db3e7c7%2Ca8033d1c%2C3efafa04%2C5277b8fe%2Cc2234757%2Ce3b092b1%2Ccaa892c2%2C5d7eaa21%2C58d0e18d%2C77520892%2C15e88334%22%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%20import%20sympy%20as%20sp%0A%20%20%20%20from%20mograder.runtime%20import%20check%2C%20hint%0A%20%20%20%20from%20scipy.integrate%20import%20RK45%2C%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%20RK45%2C%20T%2C%20check%2C%20hint%2C%20incomplete%2C%20mo%2C%20np%2C%20plt%2C%20solve_ivp%2C%20sp%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%201%20%E2%80%94%20Solving%20PDEs%3A%20finite%20differences%20and%20the%20method%20of%20lines%0A%0A%20%20%20%20You%20already%20know%20how%20to%20integrate%20ODEs%20%24%5Cfrac%7Bd%5Cmathbf%7By%7D%7D%7Bdt%7D%20%3D%20%5Cmathbf%7Bf%7D(t%2C%20%5Cmathbf%7By%7D)%24.%0A%20%20%20%20Today%20we%20turn%20a%20**partial**%20differential%20equation%20into%20a%20(large)%20system%20of%20ODEs%20and%20hand%20it%20to%20an%0A%20%20%20%20ODE%20integrator.%20This%20is%20called%20the%20**method%20of%20lines**%3A%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%5Cpartial_t%20u(t%2Cx)%20%3D%20%5Cmathcal%7BL%7D%5Bu%5D(t%2Cx)%0A%20%20%20%20%5Cquad%5Cxrightarrow%7B%5C%3B%5Ctext%7Bdiscretize%20%7D%20x%5C%3B%7D%5Cquad%0A%20%20%20%20%5Cfrac%7Bd%20u_j%7D%7Bdt%7D%20%3D%20%5Cmathcal%7BL%7D_h%5Bu%5D_j%20%2C%20%5Cqquad%20u_j(t)%5Capprox%20u(t%2C%20x_j).%0A%20%20%20%20%24%24%0A%0A%20%20%20%20To%20make%20this%20work%20we%20need%20four%20ingredients%2C%20which%20define%20the%20plan%20for%20today%3A%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%20Finite-difference%20approximations%20of%20%24%5Cpartial_x%24%20(Exercise%201%3A%20derive%20a%204th-order%20stencil%20with%20%60sympy%60)%20%7C%2020%20min%20%7C%0A%20%20%20%20%7C%202%20%7C%20The%20scalar%20wave%20equation%20as%20a%20first-order%20system%3A%20constraints%20and%20initial%20data%20%7C%2015%20min%20%7C%0A%20%20%20%20%7C%203%20%7C%20Evolving%20the%20wave%20with%204th-order%20finite%20differences%20%2B%20Dormand%E2%80%93Prince%205(4)%20(Exercise%202)%20%7C%2015%20min%20%7C%0A%20%20%20%20%7C%204%20%7C%20Time%20steps%3A%20the%20CFL%20condition%2C%20error-based%20step%20control%2C%20dense%20output%20%7C%2015%20min%20%7C%0A%20%20%20%20%7C%205%20%7C%20Convergence%20testing%20(Exercise%203)%20%7C%2010%20min%20%7C%0A%0A%20%20%20%20**How%20the%20exercises%20work.**%20Cells%20containing%20%60%23%20YOUR%20CODE%20HERE%60%20are%20for%20you%20to%20complete.%0A%20%20%20%20A%20coloured%20box%20below%20each%20exercise%20tells%20you%20whether%20your%20code%20passes%20the%20checks%0A%20%20%20%20(amber%20%3D%20waiting%20for%20your%20code%2C%20red%20%3D%20something%20is%20wrong%2C%20green%20%3D%20all%20good).%0A%20%20%20%20Cells%20further%20down%20that%20depend%20on%20an%20exercise%20will%20tell%20you%20to%20finish%20it%20first.%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%20Finite-difference%20derivatives%0A%0A%20%20%20%20Place%20grid%20points%20%24x_j%20%3D%20j%5C%2Ch%24%20with%20spacing%20%24h%3D%5CDelta%20x%24%20and%20write%20%24u_j%20%3D%20u(x_j)%24.%0A%20%20%20%20Taylor%20expanding%20about%20%24x_j%24%2C%0A%0A%20%20%20%20%24%24%0A%20%20%20%20u_%7Bj%5Cpm1%7D%20%3D%20u_j%20%5Cpm%20h%5C%2Cu'_j%20%2B%20%5Cfrac%7Bh%5E2%7D%7B2%7Du''_j%20%5Cpm%20%5Cfrac%7Bh%5E3%7D%7B6%7Du'''_j%20%2B%20%5Cfrac%7Bh%5E4%7D%7B24%7Du''''_j%20%5Cpm%20%5Cfrac%7Bh%5E5%7D%7B120%7Du%5E%7B(5)%7D_j%20%2B%20%5Cmathcal%7BO%7D(h%5E6).%0A%20%20%20%20%24%24%0A%0A%20%20%20%20**First%20derivative.**%20Subtracting%20the%20two%20expansions%2C%20all%20*even*%20derivatives%20cancel%3A%0A%0A%20%20%20%20%24%24%0A%20%20%20%20u_%7Bj%2B1%7D%20-%20u_%7Bj-1%7D%20%3D%202h%5C%2Cu'_j%20%2B%20%5Cfrac%7Bh%5E3%7D%7B3%7Du'''_j%20%2B%20%5Cmathcal%7BO%7D(h%5E5)%0A%20%20%20%20%5Cquad%5CLongrightarrow%5Cquad%0A%20%20%20%20%5Cboxed%7Bu'_j%20%3D%20%5Cfrac%7Bu_%7Bj%2B1%7D-u_%7Bj-1%7D%7D%7B2h%7D%20%5C%3B-%5C%3B%20%5Cfrac%7Bh%5E2%7D%7B6%7Du'''_j%20%2B%20%5Cmathcal%7BO%7D(h%5E4).%7D%0A%20%20%20%20%24%24%0A%0A%20%20%20%20**Second%20derivative.**%20Adding%20the%20two%20expansions%2C%20all%20*odd*%20derivatives%20cancel%3A%0A%0A%20%20%20%20%24%24%0A%20%20%20%20u_%7Bj%2B1%7D%20%2B%20u_%7Bj-1%7D%20%3D%202u_j%20%2B%20h%5E2%20u''_j%20%2B%20%5Cfrac%7Bh%5E4%7D%7B12%7Du''''_j%20%2B%20%5Cmathcal%7BO%7D(h%5E6)%0A%20%20%20%20%5Cquad%5CLongrightarrow%5Cquad%0A%20%20%20%20%5Cboxed%7Bu''_j%20%3D%20%5Cfrac%7Bu_%7Bj%2B1%7D-2u_j%2Bu_%7Bj-1%7D%7D%7Bh%5E2%7D%20%5C%3B-%5C%3B%20%5Cfrac%7Bh%5E2%7D%7B12%7Du''''_j%20%2B%20%5Cmathcal%7BO%7D(h%5E4).%7D%0A%20%20%20%20%24%24%0A%0A%20%20%20%20Both%20are%20**second%20order**%3A%20the%20error%20scales%20as%20%24h%5E2%24%2C%20so%20halving%20%24h%24%20reduces%20the%20error%20by%20%244%24.%0A%20%20%20%20The%20symmetry%20of%20the%20centered%20stencil%20is%20what%20buys%20us%20the%20extra%20order%20%E2%80%94%20the%20one-sided%20(forward)%20difference%0A%20%20%20%20%24(u_%7Bj%2B1%7D-u_j)%2Fh%20%3D%20u'_j%20%2B%20%5Ctfrac%7Bh%7D%7B2%7Du''_j%20%2B%20%5Cdots%24%20is%20only%20first%20order.%0A%0A%20%20%20%20Let's%20verify%20the%20orders%20numerically%20for%20%24u%20%3D%20%5Csin%20x%24%20at%20%24x%3D1%24.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(np%2C%20plt)%3A%0A%20%20%20%20_x0%20%3D%201.0%0A%20%20%20%20_hs%20%3D%202.0%20**%20-np.arange(1%2C%2021)%0A%20%20%20%20_err_d1%20%3D%20np.abs((np.sin(_x0%20%2B%20_hs)%20-%20np.sin(_x0%20-%20_hs))%20%2F%20(2%20*%20_hs)%20-%20np.cos(_x0))%0A%20%20%20%20_err_d2%20%3D%20np.abs(%0A%20%20%20%20%20%20%20%20(np.sin(_x0%20%2B%20_hs)%20-%202%20*%20np.sin(_x0)%20%2B%20np.sin(_x0%20-%20_hs))%20%2F%20_hs**2%20%2B%20np.sin(_x0)%0A%20%20%20%20)%0A%0A%20%20%20%20_fig%2C%20_ax%20%3D%20plt.subplots(figsize%3D(6%2C%204))%0A%20%20%20%20_ax.loglog(_hs%2C%20_err_d1%2C%20%22o-%22%2C%20label%3Dr%22%24u'%24%3A%20%24(u_%7Bj%2B1%7D-u_%7Bj-1%7D)%2F2h%24%22)%0A%20%20%20%20_ax.loglog(_hs%2C%20_err_d2%2C%20%22s-%22%2C%20label%3Dr%22%24u''%24%3A%20%24(u_%7Bj%2B1%7D-2u_j%2Bu_%7Bj-1%7D)%2Fh%5E2%24%22)%0A%20%20%20%20_ax.loglog(_hs%2C%200.1%20*%20_hs**2%2C%20%22k--%22%2C%20label%3Dr%22%24%5Cpropto%20h%5E2%24%22)%0A%20%20%20%20_ax.set_xlabel(%22%24h%24%22)%0A%20%20%20%20_ax.set_ylabel(%22absolute%20error%22)%0A%20%20%20%20_ax.set_title(r%22Second-order%20centered%20stencils%2C%20%24u%3D%5Csin%20x%24%20at%20%24x%3D1%24%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%20Notice%20that%20for%20very%20small%20%24h%24%20the%20errors%20stop%20decreasing%20and%20start%20growing%20again%3A%0A%20%20%20%20the%20stencils%20subtract%20nearly%20equal%20numbers%20and%20divide%20by%20%24h%24%20or%20%24h%5E2%24%2C%20so%20**round-off%20error**%0A%20%20%20%20(%24%5Csim%20%5Cepsilon_%7B%5Crm%20mach%7D%2Fh%24%20for%20%24u'%24%2C%20%24%5Csim%20%5Cepsilon_%7B%5Crm%20mach%7D%2Fh%5E2%24%20for%20%24u''%24)%20eventually%20wins.%20Higher-order%20stencils%20reach%20a%20small%0A%20%20%20%20truncation%20error%20at%20much%20larger%20%24h%24%2C%20before%20round-off%20matters.%0A%0A%20%20%20%20%23%23%23%20Exercise%201%20%E2%80%94%20Derive%20the%204th-order%20centered%20first%20derivative%20with%20%60sympy%60%0A%0A%20%20%20%20Wider%20stencils%20let%20us%20cancel%20more%20terms%20of%20the%20Taylor%20series.%20We%20seek%20weights%20%24a_k%24%20such%20that%0A%0A%20%20%20%20%24%24%0A%20%20%20%20u'_j%20%5Capprox%20%5Cfrac%7B1%7D%7Bh%7D%5Csum_%7Bk%3D-2%7D%5E%7B2%7D%20a_k%5C%2C%20u_%7Bj%2Bk%7D%0A%20%20%20%20%24%24%0A%0A%20%20%20%20is%20as%20accurate%20as%20possible.%20This%20is%20the%20**method%20of%20undetermined%20coefficients**%3A%0A%0A%20%20%20%201.%20Replace%20each%20%24u_%7Bj%2Bk%7D%20%3D%20u(x_j%20%2B%20kh)%24%20by%20its%20Taylor%20series%20%24%5Csum_n%20%5Cfrac%7B(kh)%5En%7D%7Bn!%7Du%5E%7B(n)%7D_j%24.%0A%20%20%20%202.%20Collect%20the%20coefficient%20multiplying%20each%20derivative%20%24u%5E%7B(n)%7D_j%24.%0A%20%20%20%203.%20Demand%20that%20the%20combination%20equals%20%24h%5C%2Cu'_j%24%3A%20the%20coefficients%20of%0A%20%20%20%20%20%20%20%24u_j%2C%20u'_j%2C%20u''_j%2C%20u'''_j%2C%20u''''_j%24%20must%20be%20%240%2C%20h%2C%200%2C%200%2C%200%24.%20Five%20equations%2C%20five%20unknowns.%0A%20%20%20%204.%20Whatever%20is%20left%20over%20is%20the%20truncation%20error.%0A%0A%20%20%20%20Below%2C%20the%20symbol%20%60du%5Bn%5D%60%20stands%20for%20%24u%5E%7B(n)%7D(x_j)%24.%20Fill%20in%20%60taylor(k)%60%20and%20then%20compute%0A%20%20%20%20%60coeffs%60%20%24%3D%5Ba_%7B-2%7D%2Ca_%7B-1%7D%2Ca_0%2Ca_1%2Ca_2%5D%24%20and%20%60leading_error%60%2C%20the%20leading%20term%20of%0A%20%20%20%20%24%5Cfrac%7B1%7D%7Bh%7D%5Csum_k%20a_k%20u_%7Bj%2Bk%7D%20-%20u'_j%24.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(hint)%3A%0A%20%20%20%20hint(%0A%20%20%20%20%20%20%20%20r%22Write%20%24u(x%2Bkh)%3D%5Csum_%7Bn%3D0%7D%5E%7B%5Crm%20order%7D%5Cfrac%7B(kh)%5En%7D%7Bn!%7Du%5E%7B(n)%7D(x)%24.%20With%20%60du%5Bn%5D%60%20standing%20for%20%22%0A%20%20%20%20%20%20%20%20r%22%24u%5E%7B(n)%7D(x)%24%20this%20is%20a%20one-line%20%60sum(...)%60%20over%20%60range(order%20%2B%201)%60%3B%20%60sp.factorial(n)%60%20gives%20%24n!%24.%22%2C%0A%20%20%20%20%20%20%20%20r%22Build%20%60combo%20%3D%20sp.expand(sum(a_k%20*%20taylor(k)%20for%20a_k%2C%20k%20in%20zip(a%2C%20offsets)))%60%2C%20which%20should%20%22%0A%20%20%20%20%20%20%20%20r%22approximate%20%24h%5C%2Cu'(x)%24.%20Then%20%60combo.coeff(du%5Bn%5D)%60%20extracts%20the%20coefficient%20multiplying%20%24u%5E%7B(n)%7D%24.%22%2C%0A%20%20%20%20%20%20%20%20r%22Make%20a%20list%20of%20equations%20%60sp.Eq(combo.coeff(du%5Bn%5D)%2C%20h%20if%20n%20%3D%3D%201%20else%200)%60%20for%20%60n%60%20in%20%60range(5)%60%20and%20%22%0A%20%20%20%20%20%20%20%20r%22solve%20them%20with%20%60sp.solve(eqs%2C%20a)%60%2C%20which%20returns%20a%20dictionary%20%60%7Ba_k%3A%20value%7D%60.%22%2C%0A%20%20%20%20%20%20%20%20r%22The%20error%20is%20what%20is%20left%20of%20%60combo.subs(solution)%20%2F%20h%20-%20du%5B1%5D%60.%20After%20%60sp.simplify%60%20it%20should%20be%20%22%0A%20%20%20%20%20%20%20%20r%22proportional%20to%20%24h%5E4%20u%5E%7B(5)%7D%24%20(why%20does%20the%20%24h%5E5%20u%5E%7B(6)%7D%24%20term%20vanish%3F).%22%2C%0A%20%20%20%20)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(sp)%3A%0A%20%20%20%20h%20%3D%20sp.symbols(%22h%22%2C%20positive%3DTrue)%0A%20%20%20%20offsets%20%3D%20%5B-2%2C%20-1%2C%200%2C%201%2C%202%5D%0A%20%20%20%20a%20%3D%20sp.symbols(%22a_m2%20a_m1%20a_0%20a_p1%20a_p2%22)%20%20%23%20the%20unknown%20weights%20a_k%0A%20%20%20%20du%20%3D%20sp.symbols(%22u0%3A7%22)%20%20%23%20du%5Bn%5D%20represents%20the%20n-th%20derivative%20of%20u%20at%20x_j%0A%0A%0A%20%20%20%20def%20taylor(k%2C%20order%3D6)%3A%0A%20%20%20%20%20%20%20%20%22%22%22Taylor%20expansion%20of%20u(x%20%2B%20k%20h)%20about%20x%2C%20through%20the%20h**order%20term.%22%22%22%0A%20%20%20%20%20%20%20%20expansion%20%3D%20...%0A%20%20%20%20%20%20%20%20%23%20YOUR%20CODE%20HERE%0A%20%20%20%20%20%20%20%20pass%0A%20%20%20%20%20%20%20%20return%20expansion%0A%0A%20%20%20%20return%20du%2C%20h%0A%0A%0A%40app.cell%0Adef%20_()%3A%0A%20%20%20%20%23%20Compute%20%60coeffs%60%20%3D%20%5Ba_%7B-2%7D%2C%20a_%7B-1%7D%2C%20a_0%2C%20a_1%2C%20a_2%5D%20and%20%60leading_error%60.%0A%20%20%20%20%23%20Use%20underscore-prefixed%20names%20(e.g.%20_combo)%20for%20scratch%20variables%3A%20marimo%20keeps%20them%20local%20to%20this%20cell.%0A%20%20%20%20coeffs%20%3D%20...%0A%20%20%20%20leading_error%20%3D%20...%0A%20%20%20%20%23%20YOUR%20CODE%20HERE%0A%20%20%20%20pass%0A%20%20%20%20return%20coeffs%2C%20leading_error%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(check%2C%20coeffs%2C%20du%2C%20h%2C%20incomplete%2C%20leading_error%2C%20sp)%3A%0A%20%20%20%20_expected%20%3D%20%5Bsp.Rational(1%2C%2012)%2C%20sp.Rational(-2%2C%203)%2C%200%2C%20sp.Rational(2%2C%203)%2C%20sp.Rational(-1%2C%2012)%5D%0A%20%20%20%20_checks%20%3D%20%5B%5D%0A%20%20%20%20if%20not%20incomplete(coeffs%2C%20leading_error)%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(%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20len(coeffs)%20%3D%3D%205%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20and%20all(sp.simplify(c%20-%20e)%20%3D%3D%200%20for%20c%2C%20e%20in%20zip(coeffs%2C%20_expected))%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22The%20weights%20should%20be%20(1%2F12%2C%20-2%2F3%2C%200%2C%202%2F3%2C%20-1%2F12).%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%20%20%20%20(%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20sp.simplify(sp.expand(leading_error)%20%2B%20h**4%20*%20du%5B5%5D%20%2F%2030)%20%3D%3D%200%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22The%20leading%20error%20should%20be%20-h%5E4%20u%5E(5)%2F30%20(did%20you%20expand%20to%20high%20enough%20order%3F).%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%20check(%22Exercise%201%3A%20fourth-order%20centered%20stencil%22%2C%20_checks)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(coeffs%2C%20incomplete%2C%20leading_error%2C%20mo%2C%20sp)%3A%0A%20%20%20%20mo.stop(incomplete(coeffs%2C%20leading_error))%0A%20%20%20%20_terms%20%3D%20%22%20%22.join(%0A%20%20%20%20%20%20%20%20f%22%7B'%2B'%20if%20c%20%3E%3D%200%20else%20'-'%7D%20%7Bsp.latex(abs(c))%7D%5C%5C%2C%20u_%7B%7Bj%7Bk%3A%2Bd%7D%7D%7D%22.replace(%22%7Bj%2B0%7D%22%2C%20%22%7Bj%7D%22)%0A%20%20%20%20%20%20%20%20for%20c%2C%20k%20in%20zip(coeffs%2C%20%5B-2%2C%20-1%2C%200%2C%201%2C%202%5D)%0A%20%20%20%20%20%20%20%20if%20c%20!%3D%200%0A%20%20%20%20)%0A%20%20%20%20mo.md(%0A%20%20%20%20%20%20%20%20r%22%22%22%0A%20%20%20%20**Result.**%20%24%5Cdisplaystyle%20u'_j%20%3D%20%5Cfrac%7B1%7D%7Bh%7D%5Cleft(TERMS%5Cright)%20%2B%20%5Cmathcal%7BO%7D(h%5E4)%2C%5Cqquad%24%0A%20%20%20%20truncation%20error%20%24%3D%20ERR%24.%0A%0A%20%20%20%20Every%20centered%20stencil%20on%20a%20periodic%20grid%20follows%20the%20same%20recipe%3B%20the%20weights%20for%20any%20order%20and%20any%0A%20%20%20%20derivative%20can%20be%20generated%20this%20way%20(Fornberg's%20algorithm%20does%20it%20efficiently).%0A%20%20%20%20%22%22%22.replace(%22TERMS%22%2C%20_terms).replace(%22ERR%22%2C%20sp.latex(leading_error))%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%202%20%E2%80%94%20The%20scalar%20wave%20equation%20as%20a%20first-order%20system%0A%0A%20%20%20%20The%20(flat-space%2C%201d%2C%20unit%20speed)%20scalar%20wave%20equation%20is%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%5Cpartial_t%5E2%20%5CPsi%20%3D%20%5Cpartial_x%5E2%20%5CPsi%20.%0A%20%20%20%20%24%24%0A%0A%20%20%20%20It%20is%20**second%20order%20in%20time**%2C%20but%20ODE%20integrators%20(like%20Dormand%E2%80%93Prince)%20solve%20first-order%0A%20%20%20%20systems%20%24%5Cdot%7B%5Cmathbf%20y%7D%20%3D%20%5Cmathbf%20f(t%2C%5Cmathbf%20y)%24.%20You%20have%20seen%20the%20fix%20before%3A%20Newton's%20law%0A%20%20%20%20%24%5Cddot%20x%20%3D%20F%2Fm%24%20becomes%20a%20first-order%20system%20for%20%24(x%2C%20v)%24.%20Here%20we%20introduce%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%5CPi%20%5Cequiv%20-%5Cpartial_t%20%5CPsi%2C%20%5Cqquad%20%5CPhi%20%5Cequiv%20%5Cpartial_x%20%5CPsi%20%2C%0A%20%20%20%20%24%24%0A%0A%20%20%20%20(the%20minus%20sign%20in%20%24%5CPi%24%20is%20a%20convention%20common%20in%20numerical%20relativity%20codes%20such%20as%20SpECTRE)%20and%20obtain%0A%20%20%20%20the%20**first-order%20system**%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%5Cbegin%7Baligned%7D%0A%20%20%20%20%5Cpartial_t%20%5CPsi%20%26%3D%20-%5CPi%2C%20%5C%5C%0A%20%20%20%20%5Cpartial_t%20%5CPi%20%20%26%3D%20-%5Cpartial_x%20%5CPhi%2C%20%5C%5C%0A%20%20%20%20%5Cpartial_t%20%5CPhi%20%26%3D%20-%5Cpartial_x%20%5CPi%20%2B%20%5Cgamma_2%5Cleft(%5Cpartial_x%5CPsi%20-%20%5CPhi%5Cright).%0A%20%20%20%20%5Cend%7Baligned%7D%0A%20%20%20%20%24%24%0A%0A%20%20%20%20Check%3A%20%24%5Cpartial_t%5E2%5CPsi%20%3D%20-%5Cpartial_t%5CPi%20%3D%20%5Cpartial_x%5CPhi%20%3D%20%5Cpartial_x%5E2%5CPsi%24.%20%E2%9C%93%0A%0A%20%20%20%20Why%20also%20introduce%20%24%5CPhi%24%20(instead%20of%20just%20%24%5Cpartial_t%5CPi%20%3D%20-%5Cpartial_x%5E2%5CPsi%24)%3F%20A%20**fully%20first-order**%0A%20%20%20%20system%20only%20ever%20needs%20a%20first-derivative%20operator%2C%20and%20it%20is%20the%20natural%20setting%20for%20the%0A%20%20%20%20characteristic%20analysis%20and%20discontinuous%20Galerkin%20methods%20of%20Lecture%202.%20Signals%20in%20this%20system%0A%20%20%20%20travel%20with%20speeds%20%24%5Cpm1%24%20(and%20%240%24)%2C%20which%20will%20set%20our%20time%20step.%0A%0A%20%20%20%20%23%23%23%20The%20reduction%20constraint%0A%0A%20%20%20%20Introducing%20%24%5CPhi%24%20added%20information%20that%20is%20not%20independent%3A%20a%20solution%20of%20the%20first-order%20system%20is%20only%20a%0A%20%20%20%20solution%20of%20the%20wave%20equation%20if%20the%20**constraint**%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%5Cmathcal%7BC%7D%20%5Cequiv%20%5Cpartial_x%5CPsi%20-%20%5CPhi%20%3D%200%0A%20%20%20%20%24%24%0A%0A%20%20%20%20holds.%20How%20does%20%24%5Cmathcal%20C%24%20evolve%3F%20Using%20the%20evolution%20equations%2C%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%5Cpartial_t%20%5Cmathcal%7BC%7D%20%3D%20%5Cpartial_x%5Cpartial_t%5CPsi%20-%20%5Cpartial_t%5CPhi%0A%20%20%20%20%3D%20-%5Cpartial_x%5CPi%20-%20%5Cleft(-%5Cpartial_x%5CPi%20%2B%20%5Cgamma_2%5C%2C%5Cmathcal%7BC%7D%5Cright)%20%3D%20-%5Cgamma_2%5C%2C%5Cmathcal%7BC%7D%0A%20%20%20%20%5Cquad%5CLongrightarrow%5Cquad%0A%20%20%20%20%5Cmathcal%7BC%7D(t%2Cx)%20%3D%20%5Cmathcal%7BC%7D(0%2Cx)%5C%2Ce%5E%7B-%5Cgamma_2%20t%7D.%0A%20%20%20%20%24%24%0A%0A%20%20%20%20*%20Constraint%20violations%20**do%20not%20propagate**%20(their%20speed%20is%20zero).%0A%20%20%20%20*%20With%20%24%5Cgamma_2%3D0%24%20any%20violation%20%E2%80%94%20from%20imperfect%20initial%20data%20or%20accumulated%20truncation%20error%20%E2%80%94%20**stays%20forever**%0A%20%20%20%20%20%20(and%20in%20more%20complicated%20systems%2C%20e.g.%20Einstein's%20equations%2C%20it%20can%20grow%20exponentially).%0A%20%20%20%20*%20The%20term%20%24%5Cgamma_2(%5Cpartial_x%5CPsi-%5CPhi)%24%20vanishes%20for%20true%20solutions%2C%20so%20it%20does%20not%20change%20the%20physics%2C%20but%20it%0A%20%20%20%20%20%20**damps**%20violations%20on%20a%20time%20scale%20%241%2F%5Cgamma_2%24.%20This%20is%20**constraint%20damping**.%0A%20%20%20%20%20%20(In%202d%2F3d%20there%20is%20a%20second%20constraint%20%24%5Cpartial_%7B%5Bi%7D%5CPhi_%7Bj%5D%7D%3D0%24%2C%20which%20is%20damped%20the%20same%20way.)%0A%20%20%20%20*%20Don't%20make%20%24%5Cgamma_2%24%20huge%3A%20the%20damping%20term%20has%20eigenvalue%20%24-%5Cgamma_2%24%2C%20which%20an%20explicit%20integrator%20must%20also%0A%20%20%20%20%20%20resolve%20(%24%5Cgamma_2%5CDelta%20t%20%5Clesssim%203%24).%0A%0A%20%20%20%20%23%23%23%20Choosing%20initial%20data%0A%0A%20%20%20%20A%20second-order-in-time%20equation%20needs%20**two**%20pieces%20of%20initial%20data%3A%20%24%5CPsi(0%2Cx)%24%20and%20%24%5Cpartial_t%5CPsi(0%2Cx)%24.%0A%20%20%20%20In%20the%20first-order%20system%3A%0A%0A%20%20%20%20*%20%24%5CPsi(0%2Cx)%24%20is%20free%3B%0A%20%20%20%20*%20%24%5CPi(0%2Cx)%20%3D%20-%5Cpartial_t%5CPsi(0%2Cx)%24%20is%20the%20second%20free%20function%3B%0A%20%20%20%20*%20%24%5CPhi(0%2Cx)%24%20is%20**not**%20free%3A%20it%20must%20satisfy%20the%20constraint%2C%20%24%5CPhi(0%2Cx)%20%3D%20%5Cpartial_x%5CPsi(0%2Cx)%24.%0A%0A%20%20%20%20The%20choice%20of%20%24%5CPi%24%20determines%20where%20the%20wave%20goes.%20By%20d'Alembert%2C%20%24%5CPsi%20%3D%20f(x-t)%20%2B%20g(x%2Bt)%24%2C%20so%0A%20%20%20%20%24%5CPhi%20%3D%20f'%20%2B%20g'%24%20and%20%24%5CPi%20%3D%20f'%20-%20g'%24%3A%0A%20%20%20%20%24%5CPi%3D%5CPhi%24%20gives%20a%20purely%20**right-moving**%20wave%2C%20%24%5CPi%3D-%5CPhi%24%20a%20**left-moving**%20one%2C%20and%20%24%5CPi%3D0%24%20a%20standing%20wave.%0A%0A%20%20%20%20We%20will%20evolve%20the%20right-moving%20traveling%20wave%20on%20the%20periodic%20domain%20%24x%5Cin%5B0%2C2%5Cpi)%24%2C%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%5CPsi%20%3D%20%5Csin(x-t)%2C%5Cqquad%20%5CPi%20%3D%20%5Ccos(x-t)%2C%5Cqquad%20%5CPhi%20%3D%20%5Ccos(x-t)%2C%0A%20%20%20%20%24%24%0A%0A%20%20%20%20whose%20period%20is%20%24T%3D2%5Cpi%24.%20Knowing%20the%20exact%20solution%20at%20all%20times%20makes%20it%20a%20perfect%20test%20problem.%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%203%20%E2%80%94%20Evolving%20the%20wave%20with%204th-order%20finite%20differences%0A%0A%20%20%20%20%23%23%23%20Exercise%202%20%E2%80%94%20The%20semi-discrete%20system%0A%0A%20%20%20%20Use%20%24N%24%20grid%20points%20%24x_j%20%3D%20j%5C%2C%5CDelta%20x%24%2C%20%24j%3D0%2C%5Cdots%2CN-1%24%2C%20with%20%24%5CDelta%20x%20%3D%202%5Cpi%2FN%24.%0A%20%20%20%20Periodicity%20means%20%24u_%7BN%7D%20%3D%20u_0%24%2C%20%24u_%7B-1%7D%3Du_%7BN-1%7D%24%2C%20etc.%2C%20which%20%60np.roll%60%20handles%20for%20us%3A%0A%20%20%20%20%60np.roll(u%2C%20-1)%5Bj%5D%20%3D%3D%20u%5Bj%2B1%5D%60%20and%20%60np.roll(u%2C%201)%5Bj%5D%20%3D%3D%20u%5Bj-1%5D%60.%0A%0A%20%20%20%20The%20ODE%20state%20vector%20stacks%20the%20three%20fields%2C%20%24%5Cmathbf%7By%7D%20%3D%20%5B%5CPsi_0%5Cdots%5CPsi_%7BN-1%7D%2C%5C%20%5CPi_0%5Cdots%5CPi_%7BN-1%7D%2C%5C%20%5CPhi_0%5Cdots%5CPhi_%7BN-1%7D%5D%24%2C%0A%20%20%20%20and%20%60np.split(y%2C%203)%60%20recovers%20them.%0A%0A%20%20%20%20Fill%20in%3A%0A%0A%20%20%20%201.%20%60d_dx_fd4(u%2C%20dx)%60%3A%20the%204th-order%20centered%20derivative%20from%20Exercise%201%20(periodic)%3B%0A%20%20%20%202.%20%60initial_data(x)%60%3A%20%24(%5CPsi%2C%20%5CPi%2C%20%5CPhi)%24%20at%20%24t%3D0%24%20for%20the%20right-moving%20wave%3B%0A%20%20%20%203.%20%60rhs(t%2C%20y%2C%20dx%2C%20gamma2)%60%3A%20%24d%5Cmathbf%7By%7D%2Fdt%24%20for%20the%20first-order%20system%2C%20including%20the%20constraint-damping%20term.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.function%0Adef%20d_dx_fd4(u%2C%20dx)%3A%0A%20%20%20%20%22%22%22Fourth-order%20centered%20first%20derivative%20of%20periodic%20data%20u%20with%20grid%20spacing%20dx.%0A%0A%20%20%20%20Hint%3A%20np.roll(u%2C%20-1)%5Bj%5D%20%3D%3D%20u%5Bj%20%2B%201%5D%20and%20np.roll(u%2C%201)%5Bj%5D%20%3D%3D%20u%5Bj%20-%201%5D.%0A%20%20%20%20%22%22%22%0A%20%20%20%20du_dx%20%3D%20...%0A%20%20%20%20%23%20YOUR%20CODE%20HERE%0A%20%20%20%20pass%0A%20%20%20%20return%20du_dx%0A%0A%0A%40app.function%0Adef%20initial_data(x)%3A%0A%20%20%20%20%22%22%22(psi%2C%20pi%2C%20phi)%20at%20t%20%3D%200%20for%20the%20right-moving%20wave%20Psi%20%3D%20sin(x%20-%20t).%0A%0A%20%20%20%20Remember%3A%20Pi%20%3D%20-dPsi%2Fdt%20and%20Phi%20%3D%20dPsi%2Fdx.%0A%20%20%20%20%22%22%22%0A%20%20%20%20psi%20%3D%20...%0A%20%20%20%20pi%20%3D%20...%0A%20%20%20%20phi%20%3D%20...%0A%20%20%20%20%23%20YOUR%20CODE%20HERE%0A%20%20%20%20pass%0A%20%20%20%20return%20psi%2C%20pi%2C%20phi%0A%0A%0A%40app.cell%0Adef%20_(np)%3A%0A%20%20%20%20def%20rhs(t%2C%20y%2C%20dx%2C%20gamma2)%3A%0A%20%20%20%20%20%20%20%20%22%22%22dy%2Fdt%20for%20the%20first-order%20scalar%20wave%20system%20on%20a%20periodic%20grid.%22%22%22%0A%20%20%20%20%20%20%20%20psi%2C%20pi%2C%20phi%20%3D%20np.split(y%2C%203)%0A%20%20%20%20%20%20%20%20%23%20Compute%20dt_psi%2C%20dt_pi%20and%20dt_phi%20using%20d_dx_fd4%2C%20then%20stack%20them%3A%0A%20%20%20%20%20%20%20%20%23%20%20%20%20%20dt_y%20%3D%20np.concatenate(%5Bdt_psi%2C%20dt_pi%2C%20dt_phi%5D)%0A%20%20%20%20%20%20%20%20dt_y%20%3D%20...%0A%20%20%20%20%20%20%20%20%23%20YOUR%20CODE%20HERE%0A%20%20%20%20%20%20%20%20pass%0A%20%20%20%20%20%20%20%20return%20dt_y%0A%0A%20%20%20%20return%20(rhs%2C)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(T%2C%20check%2C%20incomplete%2C%20np%2C%20rhs)%3A%0A%20%20%20%20def%20_grid(n)%3A%0A%20%20%20%20%20%20%20%20return%20np.arange(n)%20*%20T%20%2F%20n%2C%20T%20%2F%20n%0A%0A%20%20%20%20_checks%20%3D%20%5B%5D%0A%20%20%20%20_x32%2C%20_dx32%20%3D%20_grid(32)%0A%20%20%20%20_x64%2C%20_dx64%20%3D%20_grid(64)%0A%20%20%20%20_d_ok%20%3D%20not%20incomplete(d_dx_fd4(np.sin(_x32)%2C%20_dx32))%0A%20%20%20%20_id_ok%20%3D%20not%20incomplete(initial_data(_x64))%0A%20%20%20%20_rhs_ok%20%3D%20(%0A%20%20%20%20%20%20%20%20_d_ok%0A%20%20%20%20%20%20%20%20and%20_id_ok%0A%20%20%20%20%20%20%20%20and%20not%20incomplete(rhs(0.0%2C%20np.concatenate(initial_data(_x64))%2C%20_dx64%2C%201.0))%0A%20%20%20%20)%0A%0A%20%20%20%20if%20_d_ok%3A%0A%20%20%20%20%20%20%20%20_e32%20%3D%20np.max(np.abs(d_dx_fd4(np.sin(_x32)%2C%20_dx32)%20-%20np.cos(_x32)))%0A%20%20%20%20%20%20%20%20_e64%20%3D%20np.max(np.abs(d_dx_fd4(np.sin(_x64)%2C%20_dx64)%20-%20np.cos(_x64)))%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(_e64%20%3C%201e-5%2C%20f%22d_dx_fd4%20is%20not%20accurate%20enough%3A%20max%20error%20%7B_e64%3A.2e%7D%20for%20N%3D64%20on%20sin(x).%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%2013%20%3C%20_e32%20%2F%20_e64%20%3C%2019%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20f%22Halving%20dx%20should%20reduce%20the%20error%20by%202%5E4%20%3D%2016%2C%20but%20got%20%7B_e32%20%2F%20_e64%3A.2f%7D.%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%20_id_ok%3A%0A%20%20%20%20%20%20%20%20_psi%2C%20_pi%2C%20_phi%20%3D%20initial_data(_x64)%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(_psi%2C%20np.sin(_x64))%2C%20%22Psi(0%2C%20x)%20should%20be%20sin(x).%22)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20(np.allclose(_pi%2C%20np.cos(_x64))%2C%20%22Pi(0%2C%20x)%20%3D%20-dPsi%2Fdt%20should%20be%20cos(x)%20for%20sin(x%20-%20t).%22)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20(np.allclose(_phi%2C%20np.cos(_x64))%2C%20%22Phi(0%2C%20x)%20%3D%20dPsi%2Fdx%20should%20be%20cos(x).%22)%2C%0A%20%20%20%20%20%20%20%20%5D%0A%20%20%20%20if%20_d_ok%20and%20_rhs_ok%3A%0A%20%20%20%20%20%20%20%20_y0%20%3D%20np.concatenate(initial_data(_x64))%0A%20%20%20%20%20%20%20%20_exact_dt%20%3D%20np.concatenate(%5B-np.cos(_x64)%2C%20np.sin(_x64)%2C%20np.sin(_x64)%5D)%0A%20%20%20%20%20%20%20%20_shift%20%3D%20np.concatenate(%5Bnp.zeros(128)%2C%201e-3%20*%20np.ones(64)%5D)%20%20%23%20constant%20offset%20in%20Phi%0A%20%20%20%20%20%20%20%20_damp%20%3D%20(rhs(0.0%2C%20_y0%20%2B%20_shift%2C%20_dx64%2C%202.0)%20-%20rhs(0.0%2C%20_y0%2C%20_dx64%2C%202.0))%5B128%3A%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(%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20np.max(np.abs(rhs(0.0%2C%20_y0%2C%20_dx64%2C%201.0)%20-%20_exact_dt))%20%3C%201e-5%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22rhs(0%2C%20y0)%20does%20not%20match%20the%20exact%20time%20derivative%20(-cos%20x%2C%20sin%20x%2C%20sin%20x).%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%20%20%20%20(%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20np.allclose(_damp%2C%20-2.0%20*%201e-3)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22The%20constraint-damping%20term%20in%20dt_phi%20is%20wrong%3A%20a%20constant%20offset%20delta%20in%20Phi%20%22%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22should%20change%20dt_phi%20by%20-gamma2%20*%20delta.%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%20ex2_done%20%3D%20_d_ok%20and%20_id_ok%20and%20_rhs_ok%0A%20%20%20%20check(%22Exercise%202%3A%20finite-difference%20scalar%20wave%22%2C%20_checks)%0A%20%20%20%20return%20(ex2_done%2C)%0A%0A%0A%40app.cell%0Adef%20_(np)%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).%22%22%22%0A%20%20%20%20%20%20%20%20return%20np.sin(x%20-%20t)%2C%20np.cos(x%20-%20t)%2C%20np.cos(x%20-%20t)%0A%0A%20%20%20%20return%20(exact_solution%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%20Time%20integration%20with%20Dormand%E2%80%93Prince%205(4)%0A%0A%20%20%20%20SciPy's%20%60solve_ivp(method%3D%22RK45%22)%60%20*is*%20the%20Dormand%E2%80%93Prince%205(4)%20embedded%20Runge%E2%80%93Kutta%20pair%20with%20adaptive%0A%20%20%20%20step-size%20control.%20%60dense_output%3DTrue%60%20asks%20it%20to%20also%20build%20a%20continuous%20interpolant%2C%20%60sol.sol(t)%60%2C%0A%20%20%20%20valid%20for%20any%20%24t%24%20in%20the%20integration%20interval%20(more%20on%20this%20in%20Part%204).%0A%0A%20%20%20%20Use%20the%20controls%20to%20change%20the%20resolution%20and%20the%20constraint-damping%20parameter%2C%20and%20scrub%20through%20time.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(T%2C%20mo)%3A%0A%20%20%20%20N_slider%20%3D%20mo.ui.slider(steps%3D%5B8%2C%2016%2C%2032%2C%2064%2C%20128%2C%20256%5D%2C%20value%3D16%2C%20label%3D%22grid%20points%20%24N%24%22)%0A%20%20%20%20gamma2_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)%0A%20%20%20%20t_slider%20%3D%20mo.ui.slider(%0A%20%20%20%20%20%20%20%200.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%20)%0A%20%20%20%20mo.hstack(%5BN_slider%2C%20gamma2_slider%2C%20t_slider%5D%2C%20justify%3D%22start%22%2C%20gap%3D2)%0A%20%20%20%20return%20N_slider%2C%20gamma2_slider%2C%20t_slider%0A%0A%0A%40app.cell%0Adef%20_(N_slider%2C%20T%2C%20ex2_done%2C%20gamma2_slider%2C%20mo%2C%20np%2C%20rhs%2C%20solve_ivp)%3A%0A%20%20%20%20mo.stop(%0A%20%20%20%20%20%20%20%20not%20ex2_done%2C%0A%20%20%20%20%20%20%20%20mo.callout(mo.md(%22%E2%AC%86%20Complete%20**Exercise%202**%20to%20run%20the%20evolution.%22)%2C%20kind%3D%22warn%22)%2C%0A%20%20%20%20)%0A%20%20%20%20evo_N%20%3D%20N_slider.value%0A%20%20%20%20evo_dx%20%3D%20T%20%2F%20evo_N%0A%20%20%20%20evo_x%20%3D%20np.arange(evo_N)%20*%20evo_dx%0A%20%20%20%20evo_sol%20%3D%20solve_ivp(%0A%20%20%20%20%20%20%20%20rhs%2C%0A%20%20%20%20%20%20%20%20(0.0%2C%202%20*%20T)%2C%0A%20%20%20%20%20%20%20%20np.concatenate(initial_data(evo_x))%2C%0A%20%20%20%20%20%20%20%20method%3D%22RK45%22%2C%0A%20%20%20%20%20%20%20%20rtol%3D1e-8%2C%0A%20%20%20%20%20%20%20%20atol%3D1e-8%2C%0A%20%20%20%20%20%20%20%20dense_output%3DTrue%2C%0A%20%20%20%20%20%20%20%20args%3D(evo_dx%2C%20gamma2_slider.value)%2C%0A%20%20%20%20)%0A%20%20%20%20mo.md(%0A%20%20%20%20%20%20%20%20f%22Evolved%20to%20%24t%3D2T%24%20with%20**%7Bevo_sol.t.size%20-%201%7D%20accepted%20steps**%20and%20%22%0A%20%20%20%20%20%20%20%20f%22**%7Bevo_sol.nfev%7D%20right-hand-side%20evaluations**%20(%60rtol%20%3D%20atol%20%3D%201e-8%60).%22%0A%20%20%20%20)%0A%20%20%20%20return%20evo_sol%2C%20evo_x%0A%0A%0A%40app.cell%0Adef%20_(T%2C%20evo_sol%2C%20evo_x%2C%20exact_solution%2C%20np%2C%20plt%2C%20t_slider)%3A%0A%20%20%20%20_t%20%3D%20t_slider.value%0A%20%20%20%20_num%20%3D%20np.split(evo_sol.sol(_t)%2C%203)%0A%20%20%20%20_ex%20%3D%20exact_solution(_t%2C%20evo_x)%0A%20%20%20%20_xf%20%3D%20np.linspace(0%2C%20T%2C%20400)%0A%20%20%20%20_exf%20%3D%20exact_solution(_t%2C%20_xf)%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%2C%205.5)%2C%20sharex%3DTrue)%0A%20%20%20%20for%20_i%2C%20(_u%2C%20_name)%20in%20enumerate(zip(_num%2C%20_names))%3A%0A%20%20%20%20%20%20%20%20_ax0.plot(_xf%2C%20_exf%5B_i%5D%2C%20color%3Df%22C%7B_i%7D%22%2C%20lw%3D1%2C%20alpha%3D0.6)%0A%20%20%20%20%20%20%20%20_ax0.plot(evo_x%2C%20_u%2C%20%22o%22%2C%20color%3Df%22C%7B_i%7D%22%2C%20ms%3D4%2C%20label%3D_name)%0A%20%20%20%20%20%20%20%20_ax1.semilogy(evo_x%2C%20np.abs(_u%20-%20_ex%5B_i%5D)%20%2B%201e-17%2C%20%22o-%22%2C%20color%3Df%22C%7B_i%7D%22%2C%20ms%3D3%2C%20label%3D_name)%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(markers%3A%20numerical%2C%20lines%3A%20exact)%22)%0A%20%20%20%20_ax0.legend(loc%3D%22upper%20right%22)%0A%20%20%20%20_ax1.set_ylabel(%22%7Cnumerical%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%20%23%23%23%20Demo%20%E2%80%94%20constraint%20damping%20in%20action%0A%0A%20%20%20%20Now%20start%20with%20initial%20data%20that%20deliberately%20**violates**%20the%20constraint%2C%0A%20%20%20%20%24%5CPhi(0%2Cx)%20%3D%20%5Ccos%20x%20%2B%200.1%5Csin%203x%20%5Cneq%20%5Cpartial_x%5CPsi(0%2Cx)%24%2C%20and%20monitor%0A%20%20%20%20%24%5C%7C%5Cmathcal%7BC%7D%5C%7C%20%3D%20%5C%7CD%5CPsi%20-%20%5CPhi%5C%7C%24%20where%20%24D%24%20is%20our%20discrete%20derivative.%0A%20%20%20%20Because%20the%20same%20%24D%24%20appears%20in%20both%20the%20constraint%20and%20the%20evolution%20equations%2C%20the%20semi-discrete%0A%20%20%20%20constraint%20obeys%20%24%5Cdot%7B%5Cmathcal%20C%7D_h%20%3D%20-%5Cgamma_2%5C%2C%5Cmathcal%7BC%7D_h%24%20*exactly*%2C%20so%20we%20should%20see%20a%0A%20%20%20%20perfect%20exponential.%20Change%20%24%5Cgamma_2%24%20with%20the%20slider%20above.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(T%2C%20ex2_done%2C%20gamma2_slider%2C%20mo%2C%20np%2C%20plt%2C%20rhs%2C%20solve_ivp)%3A%0A%20%20%20%20mo.stop(not%20ex2_done%2C%20mo.callout(mo.md(%22%E2%AC%86%20Complete%20**Exercise%202**%20first.%22)%2C%20kind%3D%22warn%22))%0A%20%20%20%20_N%20%3D%2064%0A%20%20%20%20_dx%20%3D%20T%20%2F%20_N%0A%20%20%20%20_x%20%3D%20np.arange(_N)%20*%20_dx%0A%20%20%20%20_y0%20%3D%20np.concatenate(%5Bnp.sin(_x)%2C%20np.cos(_x)%2C%20np.cos(_x)%20%2B%200.1%20*%20np.sin(3%20*%20_x)%5D)%0A%20%20%20%20_ts%20%3D%20np.linspace(0%2C%20T%2C%20200)%0A%0A%20%20%20%20_fig%2C%20_ax%20%3D%20plt.subplots(figsize%3D(6.5%2C%204))%0A%20%20%20%20for%20_g2%20in%20sorted(%7B0.0%2C%20gamma2_slider.value%7D)%3A%0A%20%20%20%20%20%20%20%20_sol%20%3D%20solve_ivp(%0A%20%20%20%20%20%20%20%20%20%20%20%20rhs%2C%20(0%2C%20T)%2C%20_y0%2C%20method%3D%22RK45%22%2C%20rtol%3D1e-10%2C%20atol%3D1e-10%2C%20dense_output%3DTrue%2C%20args%3D(_dx%2C%20_g2)%0A%20%20%20%20%20%20%20%20)%0A%20%20%20%20%20%20%20%20_C%20%3D%20%5B%5D%0A%20%20%20%20%20%20%20%20for%20_t%20in%20_ts%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20_psi%2C%20_pi%2C%20_phi%20%3D%20np.split(_sol.sol(_t)%2C%203)%0A%20%20%20%20%20%20%20%20%20%20%20%20_C.append(np.sqrt(_dx%20*%20np.sum((d_dx_fd4(_psi%2C%20_dx)%20-%20_phi)%20**%202)))%0A%20%20%20%20%20%20%20%20_ax.semilogy(_ts%20%2F%20T%2C%20_C%2C%20lw%3D2%2C%20label%3Drf%22numerical%2C%20%24%5Cgamma_2%3D%7B_g2%3A.1f%7D%24%22)%0A%20%20%20%20%20%20%20%20_ax.semilogy(_ts%20%2F%20T%2C%20_C%5B0%5D%20*%20np.exp(-_g2%20*%20_ts)%2C%20%22k--%22%2C%20lw%3D1)%0A%20%20%20%20_ax.plot(%5B%5D%2C%20%5B%5D%2C%20%22k--%22%2C%20lw%3D1%2C%20label%3Dr%22%24%5C%7C%5Cmathcal%7BC%7D(0)%5C%7C%5C%2Ce%5E%7B-%5Cgamma_2%20t%7D%24%22)%0A%20%20%20%20_ax.set_xlabel(%22%24t%2FT%24%22)%0A%20%20%20%20_ax.set_ylabel(r%22%24%5C%7C%5Cmathcal%7BC%7D%5C%7C_2%24%22)%0A%20%20%20%20_ax.set_title(%22Constraint%20violation%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%204%20%E2%80%94%20Choosing%20the%20time%20step%0A%0A%20%20%20%20%23%23%23%20The%20CFL%20condition%0A%0A%20%20%20%20Insert%20a%20Fourier%20mode%20%24u_j%20%3D%20e%5E%7Bikx_j%7D%24%20into%20the%204th-order%20stencil%3A%0A%0A%20%20%20%20%24%24%0A%20%20%20%20D%20u_j%20%3D%20%5Cfrac%7Bi%7D%7B%5CDelta%20x%7D%5C%2C%5Ckappa(k%5CDelta%20x)%5C%2Cu_j%2C%5Cqquad%0A%20%20%20%20%5Ckappa(%5Ctheta)%20%3D%20%5Cfrac%7B8%5Csin%5Ctheta%20-%20%5Csin%202%5Ctheta%7D%7B6%7D%2C%20%5Cqquad%20%5Cmax_%5Ctheta%20%5Ckappa%20%5Capprox%201.372%20.%0A%20%20%20%20%24%24%0A%0A%20%20%20%20For%20the%20wave%20system%2C%20each%20Fourier%20mode%20then%20behaves%20like%20the%20ODE%20%24%5Cdot%20y%20%3D%20%5Clambda%20y%24%20with%0A%20%20%20%20%24%5Clambda%20%3D%20%5Cpm%20i%5Ckappa%2F%5CDelta%20x%24%20(plus%20%24-%5Cgamma_2%24%20from%20the%20damping%20term).%20A%20Runge%E2%80%93Kutta%20method%20applied%20to%0A%20%20%20%20%24%5Cdot%20y%20%3D%20%5Clambda%20y%24%20gives%20%24y_%7Bn%2B1%7D%20%3D%20R(%5Clambda%5CDelta%20t)%5C%2Cy_n%24%2C%20so%20the%20scheme%20is%20stable%20only%20if%20every%0A%20%20%20%20%24%5Clambda%5CDelta%20t%24%20lies%20in%20the%20**stability%20region**%20%24%7CR(z)%7C%5Cle%201%24.%20Because%20%24%7C%5Clambda%7C%20%5Cpropto%201%2F%5CDelta%20x%24%2C%0A%20%20%20%20this%20requires%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%5CDelta%20t%20%5C%3B%5Cle%5C%3B%20C%5C%2C%5Cfrac%7B%5CDelta%20x%7D%7B%7C%5Clambda_%7B%5Crm%20char%7D%7C%7D%0A%20%20%20%20%5Cqquad%5Ctext%7B(Courant%E2%80%93Friedrichs%E2%80%93Lewy%20condition)%7D%2C%0A%20%20%20%20%24%24%0A%0A%20%20%20%20where%20%24%5Clambda_%7B%5Crm%20char%7D%3D%5Cpm1%24%20are%20the%20characteristic%20speeds%20and%20%24C%24%20depends%20on%20both%20the%20spatial%20stencil%0A%20%20%20%20and%20the%20time%20integrator.%20Physically%3A%20in%20one%20step%2C%20information%20must%20not%20travel%20farther%20than%20the%20stencil%20can%20%22see%22.%0A%20%20%20%20Consequence%3A%20doubling%20the%20resolution%20halves%20the%20time%20step%2C%20so%20the%20cost%20grows%20like%20%24N%5E%7B2%7D%24%20in%201d%0A%20%20%20%20(%24N%5E%7Bd%2B1%7D%24%20in%20%24d%24%20dimensions).%0A%0A%20%20%20%20Below%2C%20the%20stability%20region%20of%20Dormand%E2%80%93Prince%20is%20computed%20directly%20from%20SciPy's%20Butcher%20tableau%0A%20%20%20%20(%60RK45.A%60%2C%20%60RK45.B%60)%20via%20%24R(z)%20%3D%201%20%2B%20z%5C%2C%5Cmathbf%20b%5ET(%5Cmathbb%201%20-%20z%20A)%5E%7B-1%7D%5Cmathbf%201%24%2C%0A%20%20%20%20and%20overlaid%20with%20the%20scaled%20eigenvalues%20%24%5Clambda%5CDelta%20t%24%20of%20our%20discretization.%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%20courant_slider%20%3D%20mo.ui.slider(%0A%20%20%20%20%20%20%20%200.1%2C%203.0%2C%20step%3D0.05%2C%20value%3D0.7%2C%20label%3Dr%22Courant%20number%20%24%5CDelta%20t%2F%5CDelta%20x%24%22%2C%20show_value%3DTrue%0A%20%20%20%20)%0A%20%20%20%20courant_slider%0A%20%20%20%20return%20(courant_slider%2C)%0A%0A%0A%40app.cell%0Adef%20_(RK45%2C%20courant_slider%2C%20np%2C%20plt)%3A%0A%20%20%20%20_b%20%3D%20RK45.B%0A%20%20%20%20_A%20%3D%20np.zeros((_b.size%2C%20_b.size))%0A%20%20%20%20_A%5B%3A%2C%20%3A%20RK45.A.shape%5B1%5D%5D%20%3D%20RK45.A%0A%20%20%20%20_ones%20%3D%20np.ones(_b.size)%0A%0A%0A%20%20%20%20def%20stability_function(z)%3A%0A%20%20%20%20%20%20%20%20return%201%20%2B%20z%20*%20_b%20%40%20np.linalg.solve(np.eye(_b.size)%20-%20z%20*%20_A%2C%20_ones)%0A%0A%0A%20%20%20%20_X%2C%20_Y%20%3D%20np.meshgrid(np.linspace(-4.5%2C%201.5%2C%20301)%2C%20np.linspace(-4%2C%204%2C%20301))%0A%20%20%20%20_absR%20%3D%20np.vectorize(lambda%20z%3A%20abs(stability_function(z)))(_X%20%2B%201j%20*%20_Y)%0A%0A%20%20%20%20%23%20Eigenvalues%20of%20the%20FD4%20wave%20system%20for%20all%20Fourier%20modes%20(gamma2%20%3D%200)%2C%20scaled%20by%20dt%0A%20%20%20%20_theta%20%3D%20np.linspace(-np.pi%2C%20np.pi%2C%20401)%0A%20%20%20%20_kappa%20%3D%20(8%20*%20np.sin(_theta)%20-%20np.sin(2%20*%20_theta))%20%2F%206%0A%20%20%20%20_zs%20%3D%201j%20*%20_kappa%20*%20courant_slider.value%0A%20%20%20%20_amp%20%3D%20max(abs(stability_function(z))%20for%20z%20in%20_zs)%0A%0A%20%20%20%20_fig%2C%20_ax%20%3D%20plt.subplots(figsize%3D(6%2C%205))%0A%20%20%20%20_ax.contourf(_X%2C%20_Y%2C%20_absR%2C%20levels%3D%5B0%2C%201%5D%2C%20colors%3D%5B%22C0%22%5D%2C%20alpha%3D0.25)%0A%20%20%20%20_ax.contour(_X%2C%20_Y%2C%20_absR%2C%20levels%3D%5B1%5D%2C%20colors%3D%5B%22C0%22%5D)%0A%20%20%20%20_ax.plot(_zs.real%2C%20_zs.imag%2C%20%22r.%22%2C%20ms%3D3%2C%20label%3Dr%22%24%5Clambda%5CDelta%20t%24%20(all%20modes)%22)%0A%20%20%20%20_ax.axhline(0%2C%20color%3D%22k%22%2C%20lw%3D0.5)%0A%20%20%20%20_ax.axvline(0%2C%20color%3D%22k%22%2C%20lw%3D0.5)%0A%20%20%20%20_ax.set_xlabel(r%22Re%24(z)%24%22)%0A%20%20%20%20_ax.set_ylabel(r%22Im%24(z)%24%22)%0A%20%20%20%20_ax.set_aspect(%22equal%22)%0A%20%20%20%20_ax.set_title(f%22Dormand%E2%80%93Prince%20stability%20region%3B%20%20max%20%7CR%7C%20%3D%20%7B_amp%3A.6f%7D%22)%0A%20%20%20%20_ax.legend(loc%3D%22lower%20left%22)%0A%20%20%20%20_fig.tight_layout()%0A%20%20%20%20_fig%0A%20%20%20%20return%20(stability_function%2C)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo%2C%20np%2C%20stability_function)%3A%0A%20%20%20%20_ys%20%3D%20np.linspace(0.01%2C%204%2C%204000)%0A%20%20%20%20_y_stable%20%3D%20_ys%5Bnp.argmax(%5Babs(stability_function(1j%20*%20y))%20%3E%201%20%2B%201e-12%20for%20y%20in%20_ys%5D)%5D%0A%20%20%20%20mo.md(rf%22%22%22%0A%20%20%20%20The%20wave%20eigenvalues%20are%20purely%20imaginary%2C%20so%20what%20matters%20is%20how%20far%20the%20stability%20region%20extends%0A%20%20%20%20**along%20the%20imaginary%20axis**.%20For%20Dormand%E2%80%93Prince%2C%20%24%7CR(iy)%7C%5Cle%201%24%20only%20up%20to%20%24y%5Capprox%20%7B_y_stable%3A.2f%7D%24%3B%0A%20%20%20%20beyond%20that%20%24%7CR(iy)%7C%24%20exceeds%20%241%24%20only%20very%20slightly%20at%20first%20(by%20%24%5Csim%2010%5E%7B%7B-3%7D%7D%24%20at%20%24y%3D1.5%24)%2C%20which%20means%20the%0A%20%20%20%20highest-frequency%20modes%20grow%20slowly%20rather%20than%20blowing%20up%20immediately.%20The%20strict%20limit%20is%20therefore%0A%20%20%20%20%24%5CDelta%20t%2F%5CDelta%20x%20%5Cle%20%7B_y_stable%3A.2f%7D%2F1.372%20%5Capprox%20%7B_y_stable%20%2F%201.372%3A.2f%7D%24.%0A%20%20%20%20Try%20Courant%20numbers%20around%20%240.7%24%2C%20%241.5%24%20and%20%242.5%24%20above.%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%23%23%23%20Why%20error-based%20step%20control%20is%20the%20right%20tool%20for%20smooth%20solutions%0A%0A%20%20%20%20Every%20Dormand%E2%80%93Prince%20step%20computes%20both%20a%205th-%20and%20a%204th-order%20solution%20from%20the%20same%20stages.%20Their%20difference%0A%20%20%20%20estimates%20the%20local%20error%20%24%5Cvarepsilon%24%20of%20the%20step%2C%20which%20is%20compared%20with%20the%20requested%20tolerance%0A%20%20%20%20%24%5Ctext%7Batol%7D%20%2B%20%5Ctext%7Brtol%7D%5C%2C%7Cy%7C%24.%20The%20step%20is%20accepted%20if%20%24%5Cvarepsilon%24%20is%20small%20enough%2C%20and%20the%20next%20step%20size%20is%0A%20%20%20%20chosen%20as%20%24%5CDelta%20t_%7B%5Crm%20new%7D%20%5Capprox%200.9%5C%2C%5CDelta%20t%5C%2C(%5Ctext%7Btol%7D%2F%5Cvarepsilon)%5E%7B1%2F5%7D%24.%0A%0A%20%20%20%20*%20You%20specify%20the%20**accuracy%20you%20want**%2C%20not%20a%20step%20size.%20There%20is%20no%20need%20to%20know%20%24%5Clambda_%7B%5Cmax%7D%24%20or%20%24C%24.%0A%20%20%20%20*%20The%20step%20adapts%20automatically%20as%20the%20solution%20changes.%0A%20%20%20%20*%20**The%20CFL%20limit%20is%20enforced%20automatically.**%20If%20%24%5CDelta%20t%24%20exceeds%20the%20stability%20limit%2C%20the%20highest-frequency%0A%20%20%20%20%20%20modes%20start%20to%20grow%2C%20the%20error%20estimate%20sees%20this%2C%20and%20the%20step%20is%20rejected%20or%20shrunk.%0A%20%20%20%20%20%20The%20controller%20ends%20up%20hovering%20near%20the%20stability%20boundary.%0A%20%20%20%20*%20Time%20errors%20are%20controlled%2C%20so%20we%20can%20make%20them%20negligible%20compared%20with%20the%20spatial%20error%2C%20which%20is%0A%20%20%20%20%20%20essential%20for%20convergence%20tests%20(Part%205).%0A%0A%20%20%20%20**Caveat%20%E2%80%94%20shocks.**%20The%20error%20estimate%20assumes%20the%20solution%20is%20smooth%20(Taylor-expandable)%20in%20time.%20Near%20a%20shock%20or%0A%20%20%20%20discontinuity%20the%20estimate%20is%20large%20no%20matter%20how%20small%20the%20step%20is%2C%20so%20the%20controller%20shrinks%20the%20step%0A%20%20%20%20drastically%20and%20rejects%20many%20steps.%20More%20importantly%2C%20what%20matters%20there%20is%20*nonlinear*%20stability%0A%20%20%20%20(no%20new%20oscillations%2C%20positivity)%2C%20not%20local%20truncation%20error.%20Shock-capturing%20codes%20therefore%20usually%20use%0A%20%20%20%20strong-stability-preserving%20(SSP)%20Runge%E2%80%93Kutta%20methods%20at%20a%20fixed%20CFL%20number.%20We%20will%20come%20back%20to%20this.%0A%0A%20%20%20%20**Demo.**%20For%20a%20fixed%20tolerance%2C%20measure%20the%20typical%20step%20size%20the%20controller%20chooses%20as%20we%20refine%20the%20grid.%0A%20%20%20%20On%20coarse%20grids%20the%20step%20is%20limited%20by%20**accuracy**%20(independent%20of%20%24%5CDelta%20x%24)%3B%20on%20fine%20grids%20it%20becomes%20limited%20by%0A%20%20%20%20**stability**%20(%24%5CDelta%20t%5Cpropto%5CDelta%20x%24).%20Notice%20that%20the%20controller%20settles%20at%20%24%5CDelta%20t%5Capprox%5CDelta%20x%24%2C%20slightly%0A%20%20%20%20*above*%20the%20strict%20limit%3A%20there%20the%20highest-frequency%20modes%20are%20amplified%20by%20only%20%24%7CR%7C-1%5Csim10%5E%7B-3%7D%24%20per%20step%2C%0A%20%20%20%20starting%20from%20round-off-sized%20amplitudes%2C%20and%20the%20error%20estimator%20keeps%20that%20growth%20in%20check.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(T%2C%20ex2_done%2C%20mo%2C%20np%2C%20plt%2C%20rhs%2C%20solve_ivp)%3A%0A%20%20%20%20mo.stop(not%20ex2_done%2C%20mo.callout(mo.md(%22%E2%AC%86%20Complete%20**Exercise%202**%20first.%22)%2C%20kind%3D%22warn%22))%0A%20%20%20%20_Ns%20%3D%20np.array(%5B16%2C%2032%2C%2064%2C%20128%2C%20256%2C%20512%2C%201024%2C%202048%5D)%0A%20%20%20%20_median_dt%20%3D%20%5B%5D%0A%20%20%20%20_sols%20%3D%20%7B%7D%0A%20%20%20%20for%20_N%20in%20_Ns%3A%0A%20%20%20%20%20%20%20%20_dx%20%3D%20T%20%2F%20_N%0A%20%20%20%20%20%20%20%20_x%20%3D%20np.arange(_N)%20*%20_dx%0A%20%20%20%20%20%20%20%20_sols%5B_N%5D%20%3D%20solve_ivp(%0A%20%20%20%20%20%20%20%20%20%20%20%20rhs%2C%20(0%2C%20T)%2C%20np.concatenate(initial_data(_x))%2C%20method%3D%22RK45%22%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20rtol%3D1e-8%2C%20atol%3D1e-8%2C%20args%3D(_dx%2C%201.0)%2C%0A%20%20%20%20%20%20%20%20)%0A%20%20%20%20%20%20%20%20_median_dt.append(np.median(np.diff(_sols%5B_N%5D.t)))%0A%20%20%20%20_dxs%20%3D%20T%20%2F%20_Ns%0A%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))%0A%20%20%20%20_ax0.loglog(_dxs%2C%20_median_dt%2C%20%22o-%22%2C%20label%3D%22median%20accepted%20%24%5C%5CDelta%20t%24%20(rtol%20%3D%201e-8)%22)%0A%20%20%20%20_ax0.loglog(_dxs%2C%200.73%20*%20_dxs%2C%20%22k--%22%2C%20label%3Dr%22strict%20stability%20limit%20%24%5CDelta%20t%20%3D%200.73%5C%2C%5CDelta%20x%24%22)%0A%20%20%20%20_ax0.set_xlabel(r%22%24%5CDelta%20x%24%22)%0A%20%20%20%20_ax0.set_ylabel(r%22%24%5CDelta%20t%24%22)%0A%20%20%20%20_ax0.legend()%0A%20%20%20%20_ax0.set_title(%22small%20%24%5C%5CDelta%20x%24%3A%20stability-limited%3B%20%20large%20%24%5C%5CDelta%20x%24%3A%20accuracy-limited%22%2C%20fontsize%3D10)%0A%20%20%20%20for%20_N%20in%20%5B64%2C%20512%2C%202048%5D%3A%0A%20%20%20%20%20%20%20%20_ax1.semilogy(_sols%5B_N%5D.t%5B1%3A%5D%20%2F%20T%2C%20np.diff(_sols%5B_N%5D.t)%20%2F%20(T%20%2F%20_N)%2C%20%22.-%22%2C%20ms%3D3%2C%20label%3Df%22N%20%3D%20%7B_N%7D%22)%0A%20%20%20%20_ax1.set_xlabel(%22%24t%2FT%24%22)%0A%20%20%20%20_ax1.set_ylabel(r%22%24%5CDelta%20t%2F%5CDelta%20x%24%22)%0A%20%20%20%20_ax1.set_title(%22accepted%20steps%22)%0A%20%20%20%20_ax1.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%23%23%23%20Dense%20output%3A%20getting%20the%20solution%20at%20the%20times%20*you*%20want%0A%0A%20%20%20%20The%20accepted%20step%20times%20%60sol.t%60%20are%20non-uniform%20and%20chosen%20by%20the%20controller%2C%20not%20by%20you.%20If%20you%20need%20output%20on%20a%0A%20%20%20%20uniform%20time%20grid%20(for%20plots%2C%20movies%2C%20Fourier%20analysis)%20or%20at%20specific%20times%20(comparisons%2C%20convergence%20tests)%2C%0A%20%20%20%20**do%20not**%20force%20the%20integrator%20to%20land%20on%20them%20(e.g.%20with%20a%20tiny%20%60max_step%60%2C%20or%20by%20restarting%20%60solve_ivp%60%20between%0A%20%20%20%20output%20times).%20That%20wastes%20work%20and%20restarts%20the%20step-size%20controller.%0A%0A%20%20%20%20Instead%20use%20**dense%20output**%3A%20Dormand%E2%80%93Prince%20comes%20with%20a%20free%204th-order%20continuous%20interpolant%20built%20from%20the%20stages%0A%20%20%20%20it%20already%20computed.%20%60sol.sol(t)%60%20evaluates%20it%20at%20any%20%24t%24%2C%20as%20the%20time%20slider%20above%20does%2C%20and%20%60t_eval%3D%60%20uses%20the%20same%0A%20%20%20%20interpolant%20internally.%20Its%20error%20is%20comparable%20to%20the%20error%20of%20the%20steps%20themselves.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(T%2C%20ex2_done%2C%20exact_solution%2C%20mo%2C%20np%2C%20plt%2C%20rhs%2C%20solve_ivp)%3A%0A%20%20%20%20mo.stop(not%20ex2_done%2C%20mo.callout(mo.md(%22%E2%AC%86%20Complete%20**Exercise%202**%20first.%22)%2C%20kind%3D%22warn%22))%0A%20%20%20%20_N%20%3D%2064%0A%20%20%20%20_dx%20%3D%20T%20%2F%20_N%0A%20%20%20%20_x%20%3D%20np.arange(_N)%20*%20_dx%0A%20%20%20%20_sol%20%3D%20solve_ivp(%0A%20%20%20%20%20%20%20%20rhs%2C%20(0%2C%20T)%2C%20np.concatenate(initial_data(_x))%2C%20method%3D%22RK45%22%2C%0A%20%20%20%20%20%20%20%20rtol%3D1e-4%2C%20atol%3D1e-4%2C%20dense_output%3DTrue%2C%20args%3D(_dx%2C%201.0)%2C%0A%20%20%20%20)%0A%20%20%20%20_tu%20%3D%20np.linspace(0%2C%20T%2C%20500)%0A%20%20%20%20_psi0_dense%20%3D%20_sol.sol(_tu)%5B0%5D%20%20%23%20Psi%20at%20x_0%20%3D%200%20on%20a%20uniform%20time%20grid%0A%20%20%20%20_psi0_steps%20%3D%20_sol.y%5B0%5D%20%20%23%20Psi%20at%20x_0%20%3D%200%20at%20the%20accepted%20steps%0A%0A%20%20%20%20_fig%2C%20(_ax0%2C%20_ax1)%20%3D%20plt.subplots(2%2C%201%2C%20figsize%3D(7%2C%205)%2C%20sharex%3DTrue)%0A%20%20%20%20_ax0.plot(_tu%20%2F%20T%2C%20_psi0_dense%2C%20%22-%22%2C%20label%3D%22dense%20output%20%60sol.sol(t)%60%22)%0A%20%20%20%20_ax0.plot(_sol.t%20%2F%20T%2C%20_psi0_steps%2C%20%22o%22%2C%20label%3Df%22accepted%20steps%20%60sol.t%60%20(%7B_sol.t.size%7D)%22)%0A%20%20%20%20_ax0.set_ylabel(r%22%24%5CPsi(t%2C%20x%3D0)%24%22)%0A%20%20%20%20_ax0.legend(loc%3D%22upper%20right%22)%0A%20%20%20%20_ax0.set_title(%22rtol%20%3D%20atol%20%3D%201e-4%2C%20N%20%3D%2064%22)%0A%20%20%20%20_ax1.semilogy(_tu%20%2F%20T%2C%20np.abs(_psi0_dense%20-%20exact_solution(_tu%2C%200.0)%5B0%5D)%2C%20%22-%22%2C%20label%3D%22dense%20output%22)%0A%20%20%20%20_ax1.semilogy(_sol.t%20%2F%20T%2C%20np.abs(_psi0_steps%20-%20exact_solution(_sol.t%2C%200.0)%5B0%5D)%20%2B%201e-17%2C%20%22o%22%2C%20label%3D%22steps%22)%0A%20%20%20%20_ax1.set_xlabel(%22%24t%2FT%24%22)%0A%20%20%20%20_ax1.set_ylabel(%22%7Cerror%7C%22)%0A%20%20%20%20_ax1.legend(loc%3D%22lower%20right%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---%0A%20%20%20%20%23%23%20Part%205%20%E2%80%94%20Convergence%20testing%0A%0A%20%20%20%20A%20code%20that%20produces%20pretty%20pictures%20is%20not%20necessarily%20correct.%20The%20single%20most%20important%20test%20of%20a%20numerical%0A%20%20%20%20PDE%20solver%20is%20a%20**convergence%20test**%3A%20if%20the%20error%20behaves%20like%20%24E(%5CDelta%20x)%5Capprox%20C%5C%2C%5CDelta%20x%5Ep%24%2C%20then%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%5Cfrac%7BE(%5CDelta%20x)%7D%7BE(%5CDelta%20x%2F2)%7D%20%3D%202%5Ep%0A%20%20%20%20%5Cquad%5CLongrightarrow%5Cquad%0A%20%20%20%20p%20%3D%20%5Clog_2%5Cfrac%7BE_N%7D%7BE_%7B2N%7D%7D%2C%0A%20%20%20%20%24%24%0A%0A%20%20%20%20and%20%24p%24%20must%20approach%20the%20design%20order%20(%244%24%20here)%20as%20%24N%24%20grows.%20We%20measure%20errors%20in%20the%20discrete%20%24L_2%24%20norm%0A%20%20%20%20%24%5C%7Ce%5C%7C_2%20%3D%20%5Cbig(%5CDelta%20x%5Csum_j%20e_j%5E2%5Cbig)%5E%7B1%2F2%7D%24.%0A%0A%20%20%20%20Things%20to%20get%20right%3A%0A%0A%20%20%20%201.%20**Time%20error%20must%20be%20negligible.**%20The%20spatial%20error%20is%20what%20we%20are%20testing%2C%20so%20we%20use%20%60rtol%20%3D%20atol%20%3D%201e-12%60.%0A%20%20%20%202.%20**Check%20every%20evolved%20variable.**%20A%20bug%20in%20the%20%24%5CPhi%24%20equation%20may%20barely%20show%20up%20in%20%24%5CPsi%24%20for%20a%20while.%0A%20%20%20%203.%20**Check%20at%20a%20generic%20time%2C%20not%20just%20a%20special%20one.**%20We%20check%20at%20exactly%20one%20period%20%24t%3DT%24%20*and*%20at%0A%20%20%20%20%20%20%20%24t%20%3D%201.0927%5C%2CT%24%3A%0A%20%20%20%20%20%20%20*%20at%20%24t%3DT%24%20the%20exact%20solution%20equals%20the%20initial%20data%2C%20so%20a%20code%20that%20does%20not%20evolve%20at%20all%2C%20or%20sends%20the%20wave%0A%20%20%20%20%20%20%20%20%20in%20the%20**wrong%20direction**%20(a%20sign%20error%20in%20%24%5CPi%24!)%2C%20would%20pass%3B%0A%20%20%20%20%20%20%20*%20%24t%20%3D%201.0927%5C%2CT%24%20is%20not%20a%20step%20endpoint%2C%20so%20we%20are%20also%20testing%20the%20**dense-output**%20interpolant%3B%0A%20%20%20%20%20%20%20*%20special%20times%20can%20hide%20errors%20through%20cancellation.%20Always%20also%20test%20at%20generic%20parameters.%0A%0A%20%20%20%20%23%23%23%20Exercise%203%20%E2%80%94%20convergence%20test%0A%0A%20%20%20%20%60evolve(N%2C%20t_final%2C%20...)%60%20(provided%20below)%20runs%20the%20solver%20and%20returns%20%60(sol%2C%20x)%60.%20Complete%0A%0A%20%20%20%201.%20%60l2_errors(sol%2C%20x%2C%20t)%60%3A%20the%20%24L_2%24%20errors%20of%20%24(%5CPsi%2C%5CPi%2C%5CPhi)%24%20at%20time%20%24t%24%20(use%20dense%20output%20and%20%60exact_solution%60)%3B%0A%20%20%20%202.%20the%20convergence%20loop%3A%20fill%20%60errors%5Blabel%5D%60%20(shape%20%60(len(Ns)%2C%203)%60)%20and%20%60orders%5Blabel%5D%60%0A%20%20%20%20%20%20%20(shape%20%60(len(Ns)%20-%201%2C%203)%60)%20for%20both%20check%20times.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(T%2C%20np%2C%20rhs%2C%20solve_ivp)%3A%0A%20%20%20%20def%20evolve(N%2C%20t_final%2C%20rtol%3D1e-12%2C%20atol%3D1e-12%2C%20gamma2%3D1.0%2C%20data%3DNone)%3A%0A%20%20%20%20%20%20%20%20%22%22%22Evolve%20the%20FD4%20scalar%20wave%20on%20N%20points%20to%20t_final%3B%20returns%20(sol%2C%20x).%0A%0A%20%20%20%20%20%20%20%20%60data%60%20optionally%20replaces%20%60initial_data%60%20(a%20function%20x%20-%3E%20(psi%2C%20pi%2C%20phi)).%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20%20%20%20%20dx%20%3D%20T%20%2F%20N%0A%20%20%20%20%20%20%20%20x%20%3D%20np.arange(N)%20*%20dx%0A%20%20%20%20%20%20%20%20y0%20%3D%20np.concatenate((initial_data%20if%20data%20is%20None%20else%20data)(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%20rhs%2C%20(0.0%2C%20t_final)%2C%20y0%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(dx%2C%20gamma2)%2C%0A%20%20%20%20%20%20%20%20)%0A%20%20%20%20%20%20%20%20return%20sol%2C%20x%0A%0A%20%20%20%20return%20(evolve%2C)%0A%0A%0A%40app.function%0Adef%20l2_errors(sol%2C%20x%2C%20t)%3A%0A%20%20%20%20%22%22%22L2%20norms%20of%20the%20errors%20in%20(psi%2C%20pi%2C%20phi)%20at%20time%20t.%20Returns%20an%20array%20of%20shape%20(3%2C).%22%22%22%0A%20%20%20%20dx%20%3D%20x%5B1%5D%20-%20x%5B0%5D%0A%20%20%20%20errors_t%20%3D%20...%0A%20%20%20%20%23%20YOUR%20CODE%20HERE%0A%20%20%20%20pass%0A%20%20%20%20return%20errors_t%0A%0A%0A%40app.cell%0Adef%20_(T%2C%20ex2_done%2C%20mo%2C%20np)%3A%0A%20%20%20%20mo.stop(not%20ex2_done%2C%20mo.callout(mo.md(%22%E2%AC%86%20Complete%20**Exercise%202**%20first.%22)%2C%20kind%3D%22warn%22))%0A%20%20%20%20Ns%20%3D%20np.array(%5B16%2C%2032%2C%2064%2C%20128%2C%20256%5D)%0A%20%20%20%20t_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%20errors%20%3D%20%7B%7D%20%20%23%20errors%5Blabel%5D%3A%20array%20of%20shape%20(len(Ns)%2C%203)%0A%20%20%20%20orders%20%3D%20%7B%7D%20%20%23%20orders%5Blabel%5D%3A%20array%20of%20shape%20(len(Ns)%20-%201%2C%203)%0A%20%20%20%20%23%20Tip%3A%20evolve%20each%20resolution%20once%20to%20a%20time%20past%20both%20check%20times%20and%20reuse%20it.%0A%20%20%20%20%23%20YOUR%20CODE%20HERE%0A%20%20%20%20pass%0A%20%20%20%20return%20Ns%2C%20errors%2C%20orders%2C%20t_checks%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(%0A%20%20%20%20T%2C%0A%20%20%20%20check%2C%0A%20%20%20%20evolve%2C%0A%20%20%20%20ex2_done%2C%0A%20%20%20%20exact_solution%2C%0A%20%20%20%20incomplete%2C%0A%20%20%20%20np%2C%0A%20%20%20%20orders%2C%0A%20%20%20%20t_checks%2C%0A)%3A%0A%20%20%20%20_checks%20%3D%20%5B%5D%0A%20%20%20%20_l2_ok%20%3D%20False%0A%20%20%20%20if%20ex2_done%3A%0A%20%20%20%20%20%20%20%20_sol%2C%20_x%20%3D%20evolve(32%2C%201.2%20*%20T%2C%20rtol%3D1e-10%2C%20atol%3D1e-10)%0A%20%20%20%20%20%20%20%20_e%20%3D%20l2_errors(_sol%2C%20_x%2C%201.2%20*%20T)%0A%20%20%20%20%20%20%20%20_l2_ok%20%3D%20not%20incomplete(_e)%0A%20%20%20%20%20%20%20%20if%20_l2_ok%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20_dx%20%3D%20_x%5B1%5D%20-%20_x%5B0%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20_ref%20%3D%20%5B%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20np.sqrt(_dx%20*%20np.sum((u%20-%20v)%20**%202))%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20for%20u%2C%20v%20in%20zip(np.split(_sol.sol(1.2%20*%20T)%2C%203)%2C%20exact_solution(1.2%20*%20T%2C%20_x))%0A%20%20%20%20%20%20%20%20%20%20%20%20%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20_checks.append((np.allclose(_e%2C%20_ref%2C%20rtol%3D1e-6)%2C%20%22l2_errors%20does%20not%20return%20the%20L2%20errors%20of%20(psi%2C%20pi%2C%20phi).%22))%0A%20%20%20%20if%20_l2_ok%20and%20len(orders)%20%3D%3D%20len(t_checks)%3A%0A%20%20%20%20%20%20%20%20for%20_label%20in%20t_checks%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20_p%20%3D%20np.asarray(orders%5B_label%5D)%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(%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20_p.shape%5B1%3A%5D%20%3D%3D%20(3%2C)%20and%20np.all(np.abs(_p%5B-1%5D%20-%204)%20%3C%200.2)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20f%22At%20%7B_label%7D%20the%20finest-resolution%20orders%20should%20be%20close%20to%204%2C%20got%20%7Bnp.round(_p%5B-1%5D%2C%202)%7D.%22%2C%0A%20%20%20%20%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)%0A%20%20%20%20ex3_done%20%3D%20_l2_ok%20and%20len(orders)%20%3D%3D%20len(t_checks)%0A%20%20%20%20check(%22Exercise%203%3A%20convergence%20test%22%2C%20_checks%20if%20ex3_done%20or%20_checks%20else%20%5B%5D)%0A%20%20%20%20return%20(ex3_done%2C)%0A%0A%0A%40app.cell%0Adef%20_(Ns%2C%20errors%2C%20ex3_done%2C%20mo%2C%20orders%2C%20plt%2C%20t_checks)%3A%0A%20%20%20%20mo.stop(not%20ex3_done%2C%20mo.callout(mo.md(%22%E2%AC%86%20Complete%20**Exercise%203**%20to%20see%20the%20convergence%20plot.%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_fig%2C%20_axes%20%3D%20plt.subplots(1%2C%202%2C%20figsize%3D(10%2C%204)%2C%20sharey%3DTrue)%0A%20%20%20%20for%20_ax%2C%20_label%20in%20zip(_axes%2C%20t_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_ax.loglog(Ns%2C%20errors%5B_label%5D%5B%3A%2C%20_i%5D%2C%20%22o-%22%2C%20label%3D_names%5B_i%5D)%0A%20%20%20%20%20%20%20%20_ax.loglog(Ns%2C%20errors%5B_label%5D%5B0%2C%200%5D%20*%20(Ns%20%2F%20Ns%5B0%5D)%20**%20-4.0%2C%20%22k--%22%2C%20label%3Dr%22%24%5Cpropto%20N%5E%7B-4%7D%24%22)%0A%20%20%20%20%20%20%20%20_ax.set_title(_label)%0A%20%20%20%20%20%20%20%20_ax.set_xlabel(%22%24N%24%22)%0A%20%20%20%20%20%20%20%20_ax.legend()%0A%20%20%20%20_axes%5B0%5D.set_ylabel(r%22%24L_2%24%20error%22)%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%7BNs%5B_k%5D%7D%20%E2%86%92%20%7BNs%5B_k%20%2B%201%5D%7D%20%7C%20%22%0A%20%20%20%20%20%20%20%20%2B%20%22%20%7C%20%22.join(%22%20%2F%20%22.join(f%22%7Borders%5B_l%5D%5B_k%2C%20_i%5D%3A.2f%7D%22%20for%20_l%20in%20t_checks)%20for%20_i%20in%20range(3))%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_k%20in%20range(len(Ns)%20-%201)%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%20orders**%20%24%5C%5Clog_2(E_N%2FE_%7B2N%7D)%24%2C%20shown%20as%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%20resolutions%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%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%23%23%23%20Demo%20%E2%80%94%20what%20a%20bug%20looks%20like%20at%20a%20special%20time%0A%0A%20%20%20%20Suppose%20we%20got%20the%20sign%20of%20%24%5CPi%24%20wrong%20in%20the%20initial%20data%2C%20%24%5CPi(0%2Cx)%20%3D%20-%5Ccos%20x%24.%20This%20produces%20the%20**left-moving**%0A%20%20%20%20wave%20%24%5Csin(x%2Bt)%24.%20Compare%20it%20against%20the%20(right-moving)%20exact%20solution%3A%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(T%2C%20evolve%2C%20ex3_done%2C%20mo%2C%20np)%3A%0A%20%20%20%20mo.stop(not%20ex3_done%2C%20mo.callout(mo.md(%22%E2%AC%86%20Complete%20**Exercise%203**%20first.%22)%2C%20kind%3D%22warn%22))%0A%20%20%20%20_sol%2C%20_x%20%3D%20evolve(64%2C%201.1%20*%20T%2C%20data%3Dlambda%20x%3A%20(np.sin(x)%2C%20-np.cos(x)%2C%20np.cos(x)))%0A%20%20%20%20_eT%20%3D%20l2_errors(_sol%2C%20_x%2C%20T)%0A%20%20%20%20_eG%20%3D%20l2_errors(_sol%2C%20_x%2C%201.0927%20*%20T)%0A%20%20%20%20mo.md(%0A%20%20%20%20%20%20%20%20%22%7C%20check%20time%20%7C%20%24%5C%5CPsi%24%20error%20%7C%20%24%5C%5CPi%24%20error%20%7C%20%24%5C%5CPhi%24%20error%20%7C%5Cn%7C---%7C---%7C---%7C---%7C%5Cn%22%0A%20%20%20%20%20%20%20%20f%22%7C%20%24t%3DT%24%20%7C%20%7B_eT%5B0%5D%3A.1e%7D%20%7C%20%7B_eT%5B1%5D%3A.1e%7D%20%7C%20%7B_eT%5B2%5D%3A.1e%7D%20%7C%5Cn%22%0A%20%20%20%20%20%20%20%20f%22%7C%20%24t%3D1.0927%5C%5C%2CT%24%20%7C%20%7B_eG%5B0%5D%3A.1e%7D%20%7C%20%7B_eG%5B1%5D%3A.1e%7D%20%7C%20%7B_eG%5B2%5D%3A.1e%7D%20%7C%5Cn%5Cn%22%0A%20%20%20%20%20%20%20%20%22At%20%24t%3DT%24%20the%20%24%5C%5CPsi%24%20and%20%24%5C%5CPhi%24%20errors%20look%20as%20good%20as%20a%20correct%20code's%3B%20only%20a%20generic%20time%20reveals%20the%20%22%0A%20%20%20%20%20%20%20%20%22bug%20unambiguously.%20(Here%20%24%5C%5CPi%24%20happens%20to%20catch%20it%20at%20%24t%3DT%24%20too%2C%20but%20for%20many%20bugs%20no%20variable%20would.)%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%23%23%23%20Demo%20%E2%80%94%20the%20time-integration%20error%20floor%0A%0A%20%20%20%20Repeat%20the%20convergence%20test%20with%20a%20loose%20tolerance%20%60rtol%20%3D%20atol%20%3D%201e-6%60.%20Once%20the%20spatial%20error%20drops%20below%20the%0A%20%20%20%20time-integration%20error%2C%20**convergence%20stalls**%2C%20and%20it%20can%20even%20look%20like%20the%20code%20gets%20*worse*%20with%20resolution%0A%20%20%20%20(finer%20grids%20force%20smaller%2C%20more%20numerous%20steps%2C%20each%20contributing%20error).%20Whenever%20a%20convergence%20test%0A%20%20%20%20misbehaves%2C%20first%20make%20sure%20you%20are%20not%20hitting%20a%20floor%20from%20another%20error%20source%20(time%20stepping%2C%20round-off%2C%0A%20%20%20%20boundary%20treatment%2C%20%E2%80%A6).%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(Ns%2C%20T%2C%20errors%2C%20evolve%2C%20ex3_done%2C%20mo%2C%20np%2C%20plt)%3A%0A%20%20%20%20mo.stop(not%20ex3_done%2C%20mo.callout(mo.md(%22%E2%AC%86%20Complete%20**Exercise%203**%20first.%22)%2C%20kind%3D%22warn%22))%0A%20%20%20%20_loose%20%3D%20np.array(%5Bl2_errors(*evolve(_N%2C%201.1%20*%20T%2C%20rtol%3D1e-6%2C%20atol%3D1e-6)%2C%201.0927%20*%20T)%20for%20_N%20in%20Ns%5D)%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_fig%2C%20_ax%20%3D%20plt.subplots(figsize%3D(6.5%2C%204))%0A%20%20%20%20for%20_i%20in%20range(3)%3A%0A%20%20%20%20%20%20%20%20_ax.loglog(Ns%2C%20errors%5B%22t%20%3D%201.0927%20T%22%5D%5B%3A%2C%20_i%5D%2C%20%22o-%22%2C%20color%3Df%22C%7B_i%7D%22%2C%20label%3Df%22%7B_names%5B_i%5D%7D%2C%20tol%20%3D%201e-12%22)%0A%20%20%20%20%20%20%20%20_ax.loglog(Ns%2C%20_loose%5B%3A%2C%20_i%5D%2C%20%22x--%22%2C%20color%3Df%22C%7B_i%7D%22%2C%20label%3Df%22%7B_names%5B_i%5D%7D%2C%20tol%20%3D%201e-6%22)%0A%20%20%20%20_ax.set_xlabel(%22%24N%24%22)%0A%20%20%20%20_ax.set_ylabel(r%22%24L_2%24%20error%20at%20%24t%3D1.0927%5C%2CT%24%22)%0A%20%20%20%20_ax.legend(ncol%3D2)%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%20Summary%0A%0A%20%20%20%20*%20**Method%20of%20lines%3A**%20discretize%20space%20%E2%86%92%20large%20ODE%20system%20%E2%86%92%20use%20a%20good%20ODE%20integrator.%0A%20%20%20%20*%20**Finite%20differences**%20approximate%20the%20*derivative*%20from%20nearby%20point%20values%3B%20a%20stencil%20of%20width%20%242m%2B1%24%0A%20%20%20%20%20%20gives%20order%20%242m%24%2C%20derived%20by%20matching%20Taylor%20series%20(Exercise%201).%0A%20%20%20%20*%20Second-order-in-time%20equations%20are%20reduced%20to%20**first%20order**%2C%20which%20introduces%20**constraints**%3B%20we%20damp%20their%0A%20%20%20%20%20%20violations%20with%20%24%5Cgamma_2%24.%20Initial%20data%20must%20satisfy%20the%20constraints%3B%20the%20free%20data%20choose%20the%20physics%0A%20%20%20%20%20%20(e.g.%20the%20direction%20of%20propagation).%0A%20%20%20%20*%20Explicit%20time%20stepping%20obeys%20a%20**CFL%20condition**%20%24%5CDelta%20t%5Clesssim%20C%5C%2C%5CDelta%20x%24.%20Error-based%20adaptive%20stepping%0A%20%20%20%20%20%20(Dormand%E2%80%93Prince%205(4))%20picks%20the%20step%20for%20you%20and%20respects%20stability%20automatically%20for%20smooth%20solutions.%0A%20%20%20%20%20%20Use%20**dense%20output**%20for%20output%20at%20chosen%20times.%0A%20%20%20%20*%20**Always**%20do%20convergence%20tests%3A%20all%20variables%2C%20generic%20times%2C%20time%20error%20%E2%89%AA%20space%20error.%0A%0A%20%20%20%20**Next%20time%3A**%20our%20errors%20fall%20like%20%24N%5E%7B-4%7D%24.%20Can%20we%20do%20better%3F%20By%20approximating%20the%20**solution**%20itself%20with%0A%20%20%20%20high-degree%20polynomials%20instead%20of%20approximating%20derivatives%2C%20we%20will%20get%20errors%20that%20fall%20*exponentially*%20with%20%24N%24%3A%0A%20%20%20%20the%20discontinuous%20Galerkin%20method.%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