Week 4b Lab — The nonlinear pendulum

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%204b%20Lab%20%E2%80%94%20The%20nonlinear%20pendulum%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%204b%20lab%20%E2%80%94%20The%20nonlinear%20pendulum%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%20You%20will%20measure%20the%20period%20of%20the%20nonlinear%20pendulum%20numerically%2C%20compare%20it%0A%20%20%20%20to%20the%20linearised%20formula%2C%20and%20find%20the%20angle%20at%20which%20the%20approximation%20breaks%20down.%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%20from%20scipy.integrate%20import%20solve_ivp%0A%20%20%20%20return%20np%2C%20plt%2C%20solve_ivp%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20Parameters%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%20g_lab%20%3D%209.81%20%20%20%20%23%20gravitational%20acceleration%20(m%2Fs%C2%B2)%0A%20%20%20%20l_lab%20%3D%201.0%20%20%20%20%20%23%20pendulum%20length%20(m)%0A%20%20%20%20return%20g_lab%2C%20l_lab%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%20solver%0A%0A%20%20%20%20The%20nonlinear%20pendulum%20%24%5Ctheta''%20%2B%20(g%2Fl)%5Csin%5Ctheta%20%3D%200%24%20is%20rewritten%20as%0A%20%20%20%20%24u%20%3D%20%5B%5Ctheta%2C%20%5Cdot%5Ctheta%5D%24%2C%20%24u'%20%3D%20%5B%5Cdot%5Ctheta%2C%5C%2C%20-(g%2Fl)%5Csin%5Ctheta%5D%24.%0A%20%20%20%20The%20cell%20below%20solves%20it%20and%20is%20already%20complete%20%E2%80%94%20you%20will%20use%20%60sol_single%60%0A%20%20%20%20in%20the%20exercises%20that%20follow.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(g_lab%2C%20l_lab%2C%20np%2C%20solve_ivp)%3A%0A%20%20%20%20theta0_single%20%3D%20np.deg2rad(30.0)%20%20%20%20%23%20initial%20amplitude%20in%20radians%0A%0A%20%20%20%20def%20rhs_pend_lab(t%2C%20u)%3A%0A%20%20%20%20%20%20%20%20return%20%5Bu%5B1%5D%2C%20-(g_lab%20%2F%20l_lab)%20*%20np.sin(u%5B0%5D)%5D%0A%0A%20%20%20%20sol_single%20%3D%20solve_ivp(%0A%20%20%20%20%20%20%20%20rhs_pend_lab%2C%20%5B0%2C%2020%5D%2C%20%5Btheta0_single%2C%200.0%5D%2C%0A%20%20%20%20%20%20%20%20max_step%3D0.01%2C%20dense_output%3DTrue%2C%0A%20%20%20%20)%0A%20%20%20%20return%20rhs_pend_lab%2C%20sol_single%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%20measure%20the%20period%20numerically%0A%0A%20%20%20%20The%20period%20is%20the%20time%20between%20two%20successive%20peaks%20of%20%24%5Ctheta(t)%24.%0A%0A%20%20%20%20**Task.**%20Replace%20%600.0%60%20with%20the%20formula%20for%20the%20**linearised**%20period%3A%0A%20%20%20%20%5C%5BT_%5Ctext%7Blin%7D%20%3D%20%5Cfrac%7B2%5Cpi%7D%7B%5Comega%7D%2C%20%5Cqquad%20%5Comega%20%3D%20%5Csqrt%7Bg%2Fl%7D.%5C%5D%0A%20%20%20%20Use%20%60g_lab%60%2C%20%60l_lab%60%2C%20and%20%60np.sqrt(...)%60.%0A%0A%20%20%20%20The%20cell%20below%20finds%20the%20first%20two%20peaks%20of%20%60sol_single%60%20automatically%20and%0A%20%20%20%20reports%20the%20numerical%20period%20%E2%80%94%20compare%20the%20two.%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%20T_lin%20%3D%200.0%20%20%20%20%23%20%E2%86%90%20replace%20with%20%202%20*%20np.pi%20%2F%20np.sqrt(g_lab%20%2F%20l_lab)%0A%20%20%20%20T_lin%0A%20%20%20%20return%20(T_lin%2C)%0A%0A%0A%40app.cell%0Adef%20_(np%2C%20sol_single)%3A%0A%20%20%20%20_theta%20%3D%20sol_single.y%5B0%5D%0A%20%20%20%20_t%20%20%20%20%20%3D%20sol_single.t%0A%20%20%20%20_peaks%20%3D%20np.where((_theta%5B1%3A-1%5D%20%3E%20_theta%5B%3A-2%5D)%20%26%20(_theta%5B1%3A-1%5D%20%3E%20_theta%5B2%3A%5D))%5B0%5D%20%2B%201%0A%20%20%20%20T_nonlin_single%20%3D%20float(_t%5B_peaks%5B1%5D%5D%20-%20_t%5B_peaks%5B0%5D%5D)%20if%20len(_peaks)%20%3E%3D%202%20else%20float(%22nan%22)%0A%20%20%20%20T_nonlin_single%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%20period%20vs%20amplitude%0A%0A%20%20%20%20Now%20measure%20the%20period%20at%20many%20amplitudes%20from%205%C2%B0%20to%2080%C2%B0.%0A%0A%20%20%20%20**Task.**%20The%20loop%20below%20is%20almost%20complete.%20%20Replace%20the%20%600.0%60%20placeholder%3A%0A%20%20%20%20%60_theta0_r%60%20should%20be%20%60_theta0_deg%60%20converted%20to%20radians%20using%20%60np.deg2rad%60.%0A%0A%20%20%20%20When%20you%20run%20the%20cell%20the%20plot%20will%20show%20both%20the%20nonlinear%20period%20(computed%0A%20%20%20%20numerically)%20and%20the%20constant%20linearised%20period.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(np%2C%20rhs_pend_lab%2C%20solve_ivp)%3A%0A%20%20%20%20theta0_deg_arr%20%3D%20np.linspace(5%2C%2080%2C%2030)%0A%20%20%20%20T_nonlin_arr%20%20%20%3D%20np.zeros(len(theta0_deg_arr))%0A%0A%20%20%20%20for%20_ii%2C%20_theta0_deg%20in%20enumerate(theta0_deg_arr)%3A%0A%20%20%20%20%20%20%20%20_theta0_r%20%3D%200.0%20%20%20%20%23%20%E2%86%90%20replace%20with%20%20np.deg2rad(_theta0_deg)%0A%20%20%20%20%20%20%20%20_sol_ii%20%3D%20solve_ivp(%0A%20%20%20%20%20%20%20%20%20%20%20%20rhs_pend_lab%2C%20%5B0%2C%2030%5D%2C%20%5B_theta0_r%2C%200.0%5D%2C%20max_step%3D0.01%2C%0A%20%20%20%20%20%20%20%20)%0A%20%20%20%20%20%20%20%20_th_ii%20%20%3D%20_sol_ii.y%5B0%5D%0A%20%20%20%20%20%20%20%20_pk_ii%20%20%3D%20np.where(%0A%20%20%20%20%20%20%20%20%20%20%20%20(_th_ii%5B1%3A-1%5D%20%3E%20_th_ii%5B%3A-2%5D)%20%26%20(_th_ii%5B1%3A-1%5D%20%3E%20_th_ii%5B2%3A%5D)%0A%20%20%20%20%20%20%20%20)%5B0%5D%20%2B%201%0A%20%20%20%20%20%20%20%20T_nonlin_arr%5B_ii%5D%20%3D%20(%0A%20%20%20%20%20%20%20%20%20%20%20%20float(_sol_ii.t%5B_pk_ii%5B1%5D%5D%20-%20_sol_ii.t%5B_pk_ii%5B0%5D%5D)%0A%20%20%20%20%20%20%20%20%20%20%20%20if%20len(_pk_ii)%20%3E%3D%202%20else%20float(%22nan%22)%0A%20%20%20%20%20%20%20%20)%0A%20%20%20%20return%20T_nonlin_arr%2C%20theta0_deg_arr%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(T_lin%2C%20T_nonlin_arr%2C%20plt%2C%20theta0_deg_arr)%3A%0A%20%20%20%20_fig%2C%20_ax%20%3D%20plt.subplots(figsize%3D(8%2C%204))%0A%20%20%20%20_ax.plot(theta0_deg_arr%2C%20T_nonlin_arr%2C%20color%3D%22%23922b21%22%2C%20linewidth%3D2.5%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20label%3D%22nonlinear%20period%20(numerical)%22)%0A%20%20%20%20_ax.axhline(T_lin%2C%20color%3D%22%231a5276%22%2C%20linestyle%3D%22--%22%2C%20linewidth%3D1.8%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20label%3Drf%22linearised%20period%20%24T%20%3D%20%7BT_lin%3A.3f%7D%24%20s%22)%0A%20%20%20%20_ax.set_xlabel(r%22Initial%20amplitude%20%24%5Ctheta_0%24%20(degrees)%22%2C%20fontsize%3D13)%0A%20%20%20%20_ax.set_ylabel(r%22Period%20%24T%24%20(s)%22%2C%20fontsize%3D13)%0A%20%20%20%20_ax.set_title(%22Period%20of%20the%20pendulum%20vs%20initial%20amplitude%22%2C%20fontsize%3D13)%0A%20%20%20%20_ax.legend(fontsize%3D11)%0A%20%20%20%20_ax.grid(True%2C%20alpha%3D0.3)%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%20Exercise%203%20%E2%80%94%20where%20does%20the%20approximation%20break%20down%3F%0A%0A%20%20%20%20**Task.**%20Use%20the%20arrays%20%60theta0_deg_arr%60%20and%20%60T_nonlin_arr%60%20(and%20%60T_lin%60)%0A%20%20%20%20to%20find%20the%20smallest%20amplitude%20at%20which%20the%20nonlinear%20period%20differs%20from%0A%20%20%20%20%24T_%5Ctext%7Blin%7D%24%20by%20more%20than%205%25.%0A%0A%20%20%20%20Fill%20in%20the%20%60threshold%60%20and%20the%20condition%20inside%20%60np.where%60.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(T_lin%2C%20T_nonlin_arr%2C%20np%2C%20theta0_deg_arr)%3A%0A%20%20%20%20threshold%20%3D%200.0%20%20%20%20%23%20%E2%86%90%20replace%20with%20%200.05%20*%20T_lin%0A%20%20%20%20_idx_5pct%20%3D%20np.where(np.abs(T_nonlin_arr%20-%20T_lin)%20%3E%20threshold)%5B0%5D%0A%20%20%20%20angle_5pct%20%3D%20float(theta0_deg_arr%5B_idx_5pct%5B0%5D%5D)%20if%20len(_idx_5pct)%20%3E%200%20else%20float(%22nan%22)%0A%20%20%20%20angle_5pct%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%20energy%20conservation%0A%0A%20%20%20%20For%20the%20nonlinear%20pendulum%20the%20total%20mechanical%20energy%0A%20%20%20%20%5C%5BE%20%3D%20%5Ctfrac%7B1%7D%7B2%7Dl%5E2%5C%2C%5Cdot%5Ctheta%5E2%20%2B%20gl%5C%2C(1%20-%20%5Ccos%5Ctheta)%5C%5D%0A%20%20%20%20is%20conserved.%20%20If%20the%20solver%20is%20accurate%2C%20%24E%24%20should%20stay%20constant.%0A%0A%20%20%20%20**Task.**%0A%20%20%20%201.%20Replace%20%60np.zeros_like(sol_single.t)%60%20with%20the%20formula%20for%20%24E%24.%0A%20%20%20%20%20%20%20Use%20%60g_lab%60%2C%20%60l_lab%60%2C%20%60sol_single.y%5B0%5D%60%20(%24%5Ctheta%24)%2C%20and%20%60sol_single.y%5B1%5D%60%20(%24%5Cdot%5Ctheta%24).%0A%20%20%20%202.%20Run%20the%20cell.%20%20How%20constant%20is%20%24E%24%3F%20%20What%20is%20the%20maximum%20variation%3F%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(np%2C%20sol_single)%3A%0A%20%20%20%20E_total%20%3D%20np.zeros_like(sol_single.t)%20%20%20%20%23%20%E2%86%90%20replace%20with%20the%20formula%0A%20%20%20%20E_total%0A%20%20%20%20return%20(E_total%2C)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(E_total%2C%20plt%2C%20sol_single)%3A%0A%20%20%20%20_fig3%2C%20_ax3%20%3D%20plt.subplots(figsize%3D(8%2C%203))%0A%20%20%20%20_ax3.plot(sol_single.t%2C%20E_total%2C%20color%3D%22%231e8449%22%2C%20linewidth%3D2)%0A%20%20%20%20_ax3.set_xlabel(r%22%24t%24%20(s)%22%2C%20fontsize%3D13)%0A%20%20%20%20_ax3.set_ylabel(r%22%24E%24%20(J%2Fkg)%22%2C%20fontsize%3D13)%0A%20%20%20%20_ax3.set_title(%22Total%20mechanical%20energy%20along%20the%20nonlinear%20solution%22%2C%20fontsize%3D13)%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