Week 4a — ODEs: direction fields and the discrete origin

import%20marimo%0A%0A__generated_with%20%3D%20%220.13.0%22%0Aapp%20%3D%20marimo.App(%0A%20%20%20%20width%3D%22medium%22%2C%0A%20%20%20%20app_title%3D%22Week%204a%20%E2%80%94%20ODEs%3A%20direction%20fields%20and%20the%20discrete%20origin%22%2C%0A)%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%20Week%204a%3A%20Ordinary%20differential%20equations%20%E2%80%94%20geometry%20and%20the%20discrete%20origin%0A%20%20%20%20%23%23%20Computation%20lecture%0A%0A%20%20%20%20**MATH4120%20%E2%80%94%20Mathematical%20Modelling%20and%20Programming**%0A%20%20%20%20Lancaster%20University%2C%202026%E2%80%9327%0A%0A%20%20%20%20---%0A%0A%20%20%20%20Lectures%20this%20week%20focused%20on%20solving%20ODEs%20by%20hand%20%E2%80%94%20separation%20of%20variables%2C%0A%20%20%20%20integrating%20factors%2C%20and%20the%20auxiliary%20equation.%20%20Here%20we%20step%20back%20and%20ask%0A%20%20%20%20a%20more%20fundamental%20question%3A%20*what%20is%20an%20ODE%2C%20geometrically%3F*%0A%0A%20%20%20%20Two%20ideas.%20%20First%3A%20a%20direction%20field.%20%20Second%3A%20where%20the%20ODE%20comes%20from.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_()%3A%0A%20%20%20%20import%20numpy%20as%20np%0A%20%20%20%20import%20matplotlib.pyplot%20as%20plt%0A%20%20%20%20return%20np%2C%20plt%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%201.%20Direction%20fields%0A%0A%20%20%20%20The%20ODE%20%24y'%20%3D%20f(t%2C%20y)%24%20assigns%20a%20slope%20%24f(t%2C%20y)%24%20to%20every%20point%20%24(t%2C%20y)%24%20in%0A%20%20%20%20the%20plane.%20%20Drawing%20those%20slopes%20as%20unit%20arrows%20gives%20a%20**direction%20field**%20%E2%80%94%0A%20%20%20%20the%20geometric%20portrait%20of%20the%20equation%20before%20any%20solving.%0A%0A%20%20%20%20A%20solution%20curve%20is%20simply%20a%20curve%20that%20is%20**tangent%20to%20the%20field%20everywhere**.%0A%20%20%20%20The%20analytic%20formula%20%24y%20%3D%20y_0%20e%5E%7Bkt%7D%24%20is%20just%20a%20compact%20way%20to%20describe%20such%20a%0A%20%20%20%20curve%3B%20the%20direction%20field%20is%20the%20more%20fundamental%20object.%0A%0A%20%20%20%20Use%20the%20sliders%20below%3A%20first%20fix%20%24k%24%20and%20watch%20the%20arrows%20change%20direction%2C%20then%0A%20%20%20%20drag%20the%20**trace**%20slider%20to%20paint%20the%20solution%20curve%20one%20point%20at%20a%20time.%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%20k_field_slider%20%3D%20mo.ui.slider(-1.5%2C%201.5%2C%20step%3D0.1%2C%20value%3D0.5%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%20%20%20%20%20%20%20%20%20label%3Dr%22%24k%24%22)%0A%20%20%20%20y0_field_slider%20%3D%20mo.ui.slider(0.2%2C%203.0%2C%20step%3D0.1%2C%20value%3D1.0%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%20%20%20%20%20%20%20%20%20%20label%3Dr%22%24y_0%24%22)%0A%20%20%20%20t_paint_slider%20%3D%20mo.ui.slider(0.0%2C%206.0%2C%20step%3D0.1%2C%20value%3D1.0%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%20%20%20%20%20%20%20%20%20label%3Dr%22Trace%20to%20%24t%24%20%3D%22)%0A%20%20%20%20mo.vstack(%5Bk_field_slider%2C%20y0_field_slider%2C%20t_paint_slider%5D)%0A%20%20%20%20return%20k_field_slider%2C%20t_paint_slider%2C%20y0_field_slider%0A%0A%0A%40app.cell%0Adef%20_(k_field_slider%2C%20np%2C%20plt%2C%20t_paint_slider%2C%20y0_field_slider)%3A%0A%20%20%20%20_k_f%20%20%20%3D%20k_field_slider.value%0A%20%20%20%20_y0_f%20%20%3D%20y0_field_slider.value%0A%20%20%20%20_t_end%20%3D%20t_paint_slider.value%0A%0A%20%20%20%20%23%20Direction%20field%0A%20%20%20%20_t_g%20%3D%20np.linspace(0%2C%206%2C%2024)%0A%20%20%20%20_y_g%20%3D%20np.linspace(0.0%2C%205.0%2C%2024)%0A%20%20%20%20_T_g%2C%20_Y_g%20%3D%20np.meshgrid(_t_g%2C%20_y_g)%0A%20%20%20%20_U_g%20%3D%20np.ones_like(_T_g)%0A%20%20%20%20_V_g%20%3D%20_k_f%20*%20_Y_g%0A%20%20%20%20_norm%20%3D%20np.hypot(_U_g%2C%20_V_g)%0A%20%20%20%20_norm%5B_norm%20%3D%3D%200%5D%20%3D%201%0A%0A%20%20%20%20%23%20Faded%20family%20of%20solution%20curves%0A%20%20%20%20_t_full%20%20%3D%20np.linspace(0%2C%206%2C%20400)%0A%20%20%20%20_y0_fam%20%20%3D%20%5B0.3%2C%200.8%2C%201.5%2C%202.5%2C%204.0%5D%0A%20%20%20%20_col_fam%20%3D%20%5B%22%231a5276%22%2C%20%22%231e8449%22%2C%20%22%23922b21%22%2C%20%22%237d3c98%22%2C%20%22%23b7950b%22%5D%0A%0A%20%20%20%20%23%20Solution%20traced%20up%20to%20t_end%0A%20%20%20%20_t_trace%20%3D%20np.linspace(0%2C%20_t_end%2C%20400)%0A%20%20%20%20_y_trace%20%3D%20_y0_f%20*%20np.exp(_k_f%20*%20_t_trace)%0A%20%20%20%20_y_now%20%20%20%3D%20_y0_f%20*%20np.exp(_k_f%20*%20_t_end)%0A%0A%20%20%20%20_fig%2C%20_ax%20%3D%20plt.subplots(figsize%3D(9%2C%205))%0A%0A%20%20%20%20_ax.quiver(%0A%20%20%20%20%20%20%20%20_T_g%2C%20_Y_g%2C%20_U_g%20%2F%20_norm%2C%20_V_g%20%2F%20_norm%2C%0A%20%20%20%20%20%20%20%20color%3D%22grey%22%2C%20alpha%3D0.4%2C%20scale%3D28%2C%0A%20%20%20%20%20%20%20%20headwidth%3D3%2C%20headlength%3D4%2C%20width%3D0.003%2C%0A%20%20%20%20)%0A%0A%20%20%20%20for%20_y0i%2C%20_ci%20in%20zip(_y0_fam%2C%20_col_fam)%3A%0A%20%20%20%20%20%20%20%20_ax.plot(_t_full%2C%20_y0i%20*%20np.exp(_k_f%20*%20_t_full)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20color%3D_ci%2C%20linewidth%3D1.0%2C%20alpha%3D0.2)%0A%0A%20%20%20%20if%20_t_end%20%3E%200%3A%0A%20%20%20%20%20%20%20%20_ax.plot(_t_trace%2C%20_y_trace%2C%20color%3D%22%231a5276%22%2C%20linewidth%3D2.5%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20label%3Drf%22%24y_0%20%3D%20%7B_y0_f%3A.1f%7D%24%2C%20%20traced%20to%20%24t%20%3D%20%7B_t_end%3A.1f%7D%24%22)%0A%20%20%20%20_ax.plot(_t_end%2C%20_y_now%2C%20%22o%22%2C%20color%3D%22%23922b21%22%2C%20markersize%3D9%2C%20zorder%3D5)%0A%0A%20%20%20%20_ax.set_xlim(0%2C%206)%0A%20%20%20%20_ax.set_ylim(0%2C%205)%0A%20%20%20%20_ax.set_xlabel(r%22%24t%24%22%2C%20fontsize%3D13)%0A%20%20%20%20_ax.set_ylabel(r%22%24y%24%22%2C%20fontsize%3D13)%0A%20%20%20%20_ax.set_title(%0A%20%20%20%20%20%20%20%20rf%22Direction%20field%20for%20%24y'%20%3D%20ky%24%20%20%20(%24k%20%3D%20%7B_k_f%3A.1f%7D%24)%22%2C%0A%20%20%20%20%20%20%20%20fontsize%3D13%2C%0A%20%20%20%20)%0A%20%20%20%20_ax.legend(fontsize%3D11)%0A%20%20%20%20_ax.grid(True%2C%20alpha%3D0.2)%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%20The%20faded%20curves%20are%20other%20solution%20curves%20%E2%80%94%20each%20one%20threads%20through%20the%20arrows%0A%20%20%20%20without%20ever%20crossing%20another.%20%20(Two%20solution%20curves%20cannot%20cross%3A%20the%20field%0A%20%20%20%20assigns%20exactly%20one%20slope%20at%20each%20point%2C%20so%20a%20crossing%20would%20mean%20two%20different%0A%20%20%20%20slopes%20at%20the%20same%20point.)%0A%0A%20%20%20%20Change%20%24k%24%20through%20zero%20and%20the%20whole%20field%20flips%3A%20growth%20becomes%20decay.%0A%20%20%20%20The%20formula%20%24y%20%3D%20y_0%20e%5E%7Bkt%7D%24%20tells%20us%20the%20same%20thing%2C%20but%20the%20direction%20field%0A%20%20%20%20*shows*%20it.%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%202.%20The%20discrete%20origin%20of%20the%20ODE%0A%0A%20%20%20%20In%20Week%201%20we%20wrote%20the%20discrete%20model%20%24N_%7Bn%2B1%7D%20%3D%20r%5C%2CN_n%24%20for%20geometric%20growth.%0A%20%20%20%20Let%20%24h%24%20be%20the%20time%20step%20between%20generations%20and%20%24r%20%3D%201%20%2B%20kh%24.%20%20Then%0A%20%20%20%20%5C%5BN_%7Bn%2B1%7D%20%3D%20N_n%20%2B%20h%5C%2Ck%5C%2CN_n%20%5Cqquad%5CLongrightarrow%5Cqquad%20%5Cfrac%7BN_%7Bn%2B1%7D%20-%20N_n%7D%7Bh%7D%20%3D%20k%5C%2CN_n.%5C%5D%0A%20%20%20%20As%20%24h%20%5Cto%200%24%20the%20left%20side%20is%20the%20definition%20of%20%24%5Cmathrm%7Bd%7DN%2F%5Cmathrm%7Bd%7Dt%24%2C%20and%0A%20%20%20%20we%20recover%20%24N'%20%3D%20kN%24.%20%20**The%20ODE%20is%20the%20continuous%20limit%20of%20the%20difference%0A%20%20%20%20equation.**%0A%0A%20%20%20%20The%20update%20rule%20%24y_%7Bn%2B1%7D%20%3D%20y_n%20%2B%20h%5C%2Cf(t_n%2C%20y_n)%24%20is%20**Euler's%20method**%20%E2%80%94%20the%0A%20%20%20%20simplest%20numerical%20solver.%20%20Here%20we%20just%20watch%20it%20converge.%20%20Chapter%206%20asks%20how%0A%20%20%20%20to%20make%20it%20accurate%20and%20efficient.%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%20h_euler_slider%20%3D%20mo.ui.slider(0.05%2C%202.0%2C%20step%3D0.05%2C%20value%3D1.5%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%20%20%20%20%20%20%20%20%20label%3Dr%22Step%20size%20%24h%24%22)%0A%20%20%20%20h_euler_slider%0A%20%20%20%20return%20(h_euler_slider%2C)%0A%0A%0A%40app.cell%0Adef%20_(h_euler_slider%2C%20np%2C%20plt)%3A%0A%20%20%20%20_h%20%20%20%3D%20h_euler_slider.value%0A%20%20%20%20_k_e%20%3D%200.5%0A%20%20%20%20_y0e%20%3D%201.0%0A%20%20%20%20_T_e%20%3D%206.0%0A%0A%20%20%20%20%23%20Euler%20staircase%0A%20%20%20%20_t_eu%20%3D%20np.arange(0%2C%20_T_e%20%2B%20_h%20%2F%202%2C%20_h)%0A%20%20%20%20_y_eu%20%3D%20np.zeros(len(_t_eu))%0A%20%20%20%20_y_eu%5B0%5D%20%3D%20_y0e%0A%20%20%20%20for%20_i%20in%20range(1%2C%20len(_t_eu))%3A%0A%20%20%20%20%20%20%20%20_y_eu%5B_i%5D%20%3D%20_y_eu%5B_i%20-%201%5D%20%2B%20_h%20*%20_k_e%20*%20_y_eu%5B_i%20-%201%5D%0A%0A%20%20%20%20%23%20Analytic%0A%20%20%20%20_t_ex%20%3D%20np.linspace(0%2C%20_T_e%2C%20500)%0A%20%20%20%20_y_ex%20%3D%20_y0e%20*%20np.exp(_k_e%20*%20_t_ex)%0A%0A%20%20%20%20%23%20Error%20at%20t%20%3D%20T_e%0A%20%20%20%20_err%20%3D%20abs(_y_eu%5B-1%5D%20-%20_y0e%20*%20np.exp(_k_e%20*%20_t_eu%5B-1%5D))%0A%0A%20%20%20%20_fig_e%2C%20_ax_e%20%3D%20plt.subplots(figsize%3D(9%2C%204))%0A%20%20%20%20_ax_e.step(_t_eu%2C%20_y_eu%2C%20where%3D%22post%22%2C%20color%3D%22%231a5276%22%2C%20linewidth%3D2.0%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20label%3Drf%22Euler%20staircase%20(%24h%20%3D%20%7B_h%3A.2f%7D%24)%22)%0A%20%20%20%20_ax_e.plot(_t_ex%2C%20_y_ex%2C%20%22k--%22%2C%20linewidth%3D1.5%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20label%3Dr%22analytic%20%24y_0%5C%2Ce%5E%7Bkt%7D%24%22)%0A%20%20%20%20_ax_e.set_xlabel(r%22%24t%24%22%2C%20fontsize%3D13)%0A%20%20%20%20_ax_e.set_ylabel(r%22%24y(t)%24%22%2C%20fontsize%3D13)%0A%20%20%20%20_ax_e.set_title(%0A%20%20%20%20%20%20%20%20rf%22Euler%20staircase%20vs%20exact%20%20%20(error%20at%20%24t%20%3D%20%7B_T_e%3A.0f%7D%24%3A%20%7B_err%3A.3f%7D)%22%2C%0A%20%20%20%20%20%20%20%20fontsize%3D13%2C%0A%20%20%20%20)%0A%20%20%20%20_ax_e.legend(fontsize%3D11)%0A%20%20%20%20_ax_e.grid(True%2C%20alpha%3D0.3)%0A%20%20%20%20_fig_e.tight_layout()%0A%20%20%20%20_fig_e%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%20Drag%20%24h%24%20down%20to%20zero%3A%20the%20staircase%20converges%20to%20the%20smooth%20curve.%0A%20%20%20%20This%20is%20not%20approximation%20magic%20%E2%80%94%20it%20is%20literally%20the%20definition%20of%20the%20derivative%2C%0A%20%20%20%20reversed.%0A%0A%20%20%20%20Notice%20that%20halving%20%24h%24%20roughly%20halves%20the%20error.%20%20Euler's%20method%20has%0A%20%20%20%20**first-order%20convergence**.%20%20Chapter%206%20introduces%20methods%20where%20halving%20%24h%24%0A%20%20%20%20reduces%20the%20error%20by%20a%20factor%20of%20four%20(second-order)%20or%20sixteen%20(fourth-order).%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%20Summary%0A%0A%20%20%20%20-%20A%20first-order%20ODE%20%24y'%20%3D%20f(t%2C%20y)%24%20is%20a%20**direction%20field**%3A%20a%20slope%20at%20every%0A%20%20%20%20%20%20point.%20%20Solution%20curves%20are%20tangent%20to%20the%20field%3B%20they%20cannot%20cross%20each%20other.%0A%20%20%20%20-%20The%20ODE%20is%20the%20continuous%20limit%20of%20the%20difference%20equation%0A%20%20%20%20%20%20%24y_%7Bn%2B1%7D%20%3D%20y_n%20%2B%20h%5C%2Cf(t_n%2C%20y_n)%24%20as%20%24h%20%5Cto%200%24.%0A%20%20%20%20-%20That%20difference%20equation%20is%20Euler's%20method.%20%20Chapter%206%20is%20the%20story%20of%20how%0A%20%20%20%20%20%20to%20do%20this%20better.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_()%3A%0A%20%20%20%20import%20marimo%20as%20mo%0A%20%20%20%20return%20(mo%2C)%0A%0A%0Aif%20__name__%20%3D%3D%20%22__main__%22%3A%0A%20%20%20%20app.run()%0A