Week 4a Lab — Direction fields and Euler's method

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%20Lab%20%E2%80%94%20Direction%20fields%20and%20Euler's%20method%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%20lab%20%E2%80%94%20Direction%20fields%20and%20Euler's%20method%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%20The%20scenario%20is%20Newton's%20law%20of%20cooling%20%E2%80%94%20an%20ODE%20you%20already%20solved%20analytically%0A%20%20%20%20in%20Week%202.%20%20Here%20you%20draw%20its%20direction%20field%2C%20overlay%20your%20analytic%20solution%2C%0A%20%20%20%20and%20then%20approximate%20it%20with%20an%20Euler%20staircase.%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%20Parameters%0A%0A%20%20%20%20A%20cup%20of%20coffee%20cools%20in%20a%20room.%20%20You%20can%20change%20these%20at%20any%20point.%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%20k_cool%20%20%3D%200.1%20%20%20%20%20%23%20cooling%20rate%20(min%E2%81%BB%C2%B9)%0A%20%20%20%20T0_cool%20%3D%2085.0%20%20%20%20%23%20initial%20temperature%20(%C2%B0C)%0A%20%20%20%20Tinf%20%20%20%20%3D%2020.0%20%20%20%20%23%20ambient%20temperature%20(%C2%B0C)%0A%20%20%20%20return%20T0_cool%2C%20Tinf%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%20The%20direction%20field%0A%0A%20%20%20%20Newton's%20law%20of%20cooling%20is%0A%20%20%20%20%5C%5B%5Cfrac%7B%5Cmathrm%7Bd%7DT%7D%7B%5Cmathrm%7Bd%7Dt%7D%20%3D%20-k%5C%2C(T%20-%20T_%5Cinfty).%5C%5D%0A%20%20%20%20This%20assigns%20a%20slope%20%24f(t%2C%20T)%20%3D%20-k(T%20-%20T_%5Cinfty)%24%20to%20every%20point%20%24(t%2C%20T)%24%0A%20%20%20%20in%20the%20plane.%0A%0A%20%20%20%20The%20direction%20field%20makes%20something%20visible%20immediately%2C%20before%20any%20solving%3A%0A%20%20%20%20when%20%24T%20%3E%20T_%5Cinfty%24%20the%20slope%20is%20negative%20(coffee%20cools)%2C%20and%20when%0A%20%20%20%20%24T%20%3C%20T_%5Cinfty%24%20it%20is%20positive.%20%20Every%20solution%20must%20converge%20to%20%24T_%5Cinfty%24.%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%20Exercise%201%20%E2%80%94%20draw%20the%20direction%20field%0A%0A%20%20%20%20**Task.**%20Replace%20%60np.zeros_like(_T_g)%60%20with%20the%20correct%20expression%20for%0A%20%20%20%20%24f(t%2C%20T)%20%3D%20-k(T%20-%20T_%5Cinfty)%24.%20%20Use%20%60k_cool%60%2C%20%60Tinf%60%2C%20and%20%60_Y_g%60%0A%20%20%20%20(which%20holds%20the%20%24T%24%20values%20on%20the%20grid).%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(T0_cool%2C%20Tinf%2C%20np)%3A%0A%20%20%20%20_t_g%20%3D%20np.linspace(0%2C%2040%2C%2022)%0A%20%20%20%20_y_g%20%3D%20np.linspace(Tinf%20-%205%2C%20T0_cool%20%2B%205%2C%2022)%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%20np.zeros_like(_T_g)%20%20%20%20%23%20%E2%86%90%20replace%20with%20%20-k_cool%20*%20(_Y_g%20-%20Tinf)%0A%20%20%20%20_norm_g%20%3D%20np.hypot(_U_g%2C%20_V_g)%0A%20%20%20%20_norm_g%5B_norm_g%20%3D%3D%200%5D%20%3D%201%0A%20%20%20%20V_field%20%3D%20_V_g%20%20%20%20%23%20returned%20so%20Exercise%202%20plot%20can%20check%20it%0A%20%20%20%20T_grid%20%20%3D%20_T_g%0A%20%20%20%20Y_grid%20%20%3D%20_Y_g%0A%20%20%20%20U_field%20%3D%20_U_g%20%2F%20_norm_g%0A%20%20%20%20V_norm_field%20%3D%20_V_g%20%2F%20_norm_g%0A%20%20%20%20return%20T_grid%2C%20U_field%2C%20V_norm_field%2C%20Y_grid%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(T_grid%2C%20U_field%2C%20V_norm_field%2C%20Y_grid%2C%20plt)%3A%0A%20%20%20%20_fig1%2C%20_ax1%20%3D%20plt.subplots(figsize%3D(9%2C%205))%0A%20%20%20%20_ax1.quiver(%0A%20%20%20%20%20%20%20%20T_grid%2C%20Y_grid%2C%20U_field%2C%20V_norm_field%2C%0A%20%20%20%20%20%20%20%20color%3D%22grey%22%2C%20alpha%3D0.45%2C%20scale%3D25%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%20%20%20%20_ax1.set_xlabel(r%22%24t%24%20(min)%22%2C%20fontsize%3D13)%0A%20%20%20%20_ax1.set_ylabel(r%22%24T%24%20(%C2%B0C)%22%2C%20fontsize%3D13)%0A%20%20%20%20_ax1.set_title(%22Direction%20field%20for%20Newton's%20law%20of%20cooling%22%2C%20fontsize%3D13)%0A%20%20%20%20_ax1.grid(True%2C%20alpha%3D0.2)%0A%20%20%20%20_fig1.tight_layout()%0A%20%20%20%20_fig1%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%20Exercise%202%20%E2%80%94%20overlay%20the%20analytic%20solution%0A%0A%20%20%20%20From%20Week%202%20you%20know%20the%20analytic%20solution%3A%0A%20%20%20%20%5C%5BT(t)%20%3D%20T_%5Cinfty%20%2B%20(T_0%20-%20T_%5Cinfty)%5C%2Ce%5E%7B-kt%7D.%5C%5D%0A%0A%20%20%20%20**Task.**%20Replace%20%60np.zeros_like(t_cool)%60%20with%20this%20formula.%0A%20%20%20%20Use%20%60k_cool%60%2C%20%60T0_cool%60%2C%20%60Tinf%60%2C%20and%20%60np.exp(...)%60.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(np)%3A%0A%20%20%20%20t_cool%20%20%20%20%3D%20np.linspace(0%2C%2040%2C%20300)%0A%20%20%20%20T_analytic%20%3D%20np.zeros_like(t_cool)%20%20%20%20%23%20%E2%86%90%20replace%20with%20the%20formula%0A%20%20%20%20return%20T_analytic%2C%20t_cool%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(T_analytic%2C%20T_grid%2C%20Tinf%2C%20U_field%2C%20V_norm_field%2C%20Y_grid%2C%20plt%2C%20t_cool)%3A%0A%20%20%20%20_fig2%2C%20_ax2%20%3D%20plt.subplots(figsize%3D(9%2C%205))%0A%20%20%20%20_ax2.quiver(%0A%20%20%20%20%20%20%20%20T_grid%2C%20Y_grid%2C%20U_field%2C%20V_norm_field%2C%0A%20%20%20%20%20%20%20%20color%3D%22grey%22%2C%20alpha%3D0.35%2C%20scale%3D25%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%20%20%20%20_ax2.plot(t_cool%2C%20T_analytic%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%20label%3Dr%22analytic%20%24T(t)%24%22)%0A%20%20%20%20_ax2.axhline(Tinf%2C%20color%3D%22tomato%22%2C%20linestyle%3D%22%3A%22%2C%20linewidth%3D1.5%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20label%3Drf%22%24T_%5Cinfty%20%3D%20%7BTinf%7D%24%C2%B0C%22)%0A%20%20%20%20_ax2.set_xlabel(r%22%24t%24%20(min)%22%2C%20fontsize%3D13)%0A%20%20%20%20_ax2.set_ylabel(r%22%24T%24%20(%C2%B0C)%22%2C%20fontsize%3D13)%0A%20%20%20%20_ax2.set_title(%22Direction%20field%20with%20analytic%20solution%20overlaid%22%2C%20fontsize%3D13)%0A%20%20%20%20_ax2.legend(fontsize%3D11)%0A%20%20%20%20_ax2.grid(True%2C%20alpha%3D0.2)%0A%20%20%20%20_fig2.tight_layout()%0A%20%20%20%20_fig2%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%20Exercise%203%20%E2%80%94%20check%20tangency%0A%0A%20%20%20%20The%20solution%20curve%20should%20be%20tangent%20to%20the%20direction%20field%20at%20every%20point.%0A%20%20%20%20At%20%24t%20%3D%2010%24%2C%20the%20slope%20of%20the%20analytic%20solution%20equals%20%24f(10%2C%20T(10))%24.%0A%0A%20%20%20%20**Task.**%20Fill%20in%20the%20two%20cells%20below.%20%20%60slope_analytic%60%20is%20the%20numerical%0A%20%20%20%20derivative%20of%20%60T_analytic%60%20at%20%24t%20%3D%2010%24%3B%20%60slope_field%60%20is%20%24f(t%2C%20T)%24%20evaluated%0A%20%20%20%20at%20that%20point.%20%20They%20should%20agree%20to%20several%20decimal%20places.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(T0_cool%2C%20np%2C%20t_cool)%3A%0A%20%20%20%20%23%20index%20closest%20to%20t%20%3D%2010%0A%20%20%20%20idx10%20%3D%20int(np.argmin(np.abs(t_cool%20-%2010)))%0A%20%20%20%20dt_cool%20%3D%20t_cool%5B1%5D%20-%20t_cool%5B0%5D%0A%0A%20%20%20%20T_at_10%20%3D%20T0_cool%20*%200.0%20%20%20%20%23%20%E2%86%90%20replace%20with%20T_analytic%5Bidx10%5D%20%20(needs%20T_analytic%20in%20scope%20below)%0A%20%20%20%20slope_analytic%20%3D%200.0%20%20%20%20%20%20%20%20%23%20%E2%86%90%20replace%20with%20%20(T_analytic%5Bidx10%2B1%5D%20-%20T_analytic%5Bidx10-1%5D)%20%2F%20(2*dt_cool)%0A%20%20%20%20slope_field%20%20%20%20%3D%200.0%20%20%20%20%20%20%20%20%23%20%E2%86%90%20replace%20with%20%20-k_cool%20*%20(T_at_10%20-%20Tinf)%0A%20%20%20%20slope_analytic%2C%20slope_field%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%20Observation%0A%0A%20%20%20%20Change%20%60k_cool%60%20in%20the%20parameters%20cell.%20%20How%20does%20the%20direction%20field%20change%3F%0A%20%20%20%20What%20about%20the%20time%20it%20takes%20the%20solution%20to%20reach%20near%20%24T_%5Cinfty%24%3F%0A%0A%20%20%20%20Now%20try%20%60T0_cool%20%3D%205.0%60%20(a%20cold%20drink%20warming%20up%20in%20a%20warm%20room).%0A%20%20%20%20Does%20the%20field%20and%20solution%20still%20make%20sense%3F%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%20Extension%20%E2%80%94%20Euler's%20method%20for%20Newton's%20cooling%0A%0A%20%20%20%20The%20Euler%20update%20is%20%24T_%7Bn%2B1%7D%20%3D%20T_n%20%2B%20h%5C%2Cf(t_n%2C%20T_n)%24%2C%20where%0A%20%20%20%20%24f(t%2C%20T)%20%3D%20-k(T%20-%20T_%5Cinfty)%24.%0A%0A%20%20%20%20**Task.**%0A%20%20%20%201.%20Replace%20%600.0%60%20in%20the%20Euler%20loop%20below%20with%20the%20correct%20update.%0A%20%20%20%20%20%20%20Use%20%60T_euler%60%2C%20%60k_cool%60%2C%20%60Tinf%60%2C%20and%20%60h_eu%60.%0A%20%20%20%202.%20Change%20%60h_eu%60%20and%20observe%20how%20the%20staircase%20converges%20to%20the%20analytic%20solution.%0A%20%20%20%203.%20At%20what%20step%20size%20does%20the%20Euler%20solution%20start%20to%20diverge%3F%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(T0_cool%2C%20np)%3A%0A%20%20%20%20h_eu%20%20%20%20%3D%205.0%20%20%20%20%23%20step%20size%20in%20minutes%20%E2%80%94%20try%201.0%2C%205.0%2C%2015.0%0A%0A%20%20%20%20t_euler%20%3D%20np.arange(0%2C%2040%20%2B%20h_eu%20%2F%202%2C%20h_eu)%0A%20%20%20%20T_euler%20%3D%20np.zeros(len(t_euler))%0A%20%20%20%20T_euler%5B0%5D%20%3D%20T0_cool%0A%20%20%20%20for%20_j%20in%20range(1%2C%20len(t_euler))%3A%0A%20%20%20%20%20%20%20%20T_euler%5B_j%5D%20%3D%20T_euler%5B_j%20-%201%5D%20%2B%20h_eu%20*%200.0%20%20%20%20%23%20%E2%86%90%20replace%200.0%0A%20%20%20%20return%20T_euler%2C%20h_eu%2C%20t_euler%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(T_analytic%2C%20T_euler%2C%20Tinf%2C%20h_eu%2C%20plt%2C%20t_cool%2C%20t_euler)%3A%0A%20%20%20%20_fig3%2C%20_ax3%20%3D%20plt.subplots(figsize%3D(9%2C%204))%0A%20%20%20%20_ax3.plot(t_cool%2C%20T_analytic%2C%20%22k--%22%2C%20linewidth%3D1.5%2C%20label%3D%22analytic%22)%0A%20%20%20%20_ax3.step(t_euler%2C%20T_euler%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%20label%3Drf%22Euler%20(%24h%20%3D%20%7Bh_eu%3A.1f%7D%24%20min)%22)%0A%20%20%20%20_ax3.axhline(Tinf%2C%20color%3D%22tomato%22%2C%20linestyle%3D%22%3A%22%2C%20linewidth%3D1.2)%0A%20%20%20%20_ax3.set_xlabel(r%22%24t%24%20(min)%22%2C%20fontsize%3D13)%0A%20%20%20%20_ax3.set_ylabel(r%22%24T%24%20(%C2%B0C)%22%2C%20fontsize%3D13)%0A%20%20%20%20_ax3.set_title(%22Euler's%20method%20vs%20analytic%20solution%22%2C%20fontsize%3D13)%0A%20%20%20%20_ax3.legend(fontsize%3D11)%0A%20%20%20%20_ax3.grid(True%2C%20alpha%3D0.3)%0A%20%20%20%20_fig3.tight_layout()%0A%20%20%20%20_fig3%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