Week 4b — ODEs: 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%20%E2%80%94%20ODEs%3A%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%3A%20Ordinary%20differential%20equations%20%E2%80%94%20the%20nonlinear%20pendulum%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%20derived%20the%20small-angle%20pendulum%20%24%5Ctheta''%20%2B%20%5Cfrac%7Bg%7D%7Bl%7D%5C%2C%5Ctheta%20%3D%200%24%0A%20%20%20%20by%20replacing%20%24%5Csin%5Ctheta%24%20with%20%24%5Ctheta%24.%20%20That%20gives%20a%20clean%20analytic%20solution.%0A%20%20%20%20Here%20we%20solve%20the%20*exact*%20equation%20%24%5Ctheta''%20%2B%20%5Cfrac%7Bg%7D%7Bl%7D%5Csin%5Ctheta%20%3D%200%24%0A%20%20%20%20numerically%20and%20ask%3A%20when%20does%20the%20approximation%20break%20down%3F%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%20The%20two%20pendulums%0A%0A%20%20%20%20Both%20start%20from%20rest%20at%20angle%20%24%5Ctheta_0%24.%20%20The%20linearised%20model%20has%20analytic%0A%20%20%20%20solution%20%24%5Ctheta(t)%20%3D%20%5Ctheta_0%5Ccos(%5Comega%20t)%24%20with%20%24%5Comega%20%3D%20%5Csqrt%7Bg%2Fl%7D%24%2C%0A%20%20%20%20independent%20of%20amplitude.%20%20The%20nonlinear%20model%20is%20solved%20with%20%60solve_ivp%60%20%E2%80%94%0A%20%20%20%20rewritten%20as%20a%20two-component%20system%20%24u%20%3D%20%5B%5Ctheta%2C%20%5Cdot%5Ctheta%5D%24%3A%0A%20%20%20%20%5C%5Bu'%20%3D%20%5Cbigl%5B%5Cdot%5Ctheta%2C%5C%3B%20-(g%2Fl)%5Csin%5Ctheta%5Cbigr%5D.%5C%5D%0A%0A%20%20%20%20Drag%20the%20amplitude%20slider.%20%20At%20small%20angles%20the%20curves%20agree%3B%20at%20large%20angles%0A%20%20%20%20they%20diverge.%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%20theta0_deg_slider%20%3D%20mo.ui.slider(5%2C%2080%2C%20step%3D5%2C%20value%3D15%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%20%20%20label%3Dr%22Initial%20amplitude%20%24%5Ctheta_0%24%20(degrees)%22)%0A%20%20%20%20theta0_deg_slider%0A%20%20%20%20return%20(theta0_deg_slider%2C)%0A%0A%0A%40app.cell%0Adef%20_(np%2C%20solve_ivp%2C%20theta0_deg_slider)%3A%0A%20%20%20%20_g_p%20%3D%209.81%0A%20%20%20%20_l_p%20%3D%201.0%0A%20%20%20%20omega_pend%20%20%20%20%20%20%3D%20np.sqrt(_g_p%20%2F%20_l_p)%0A%20%20%20%20theta0_rad_pend%20%3D%20np.deg2rad(theta0_deg_slider.value)%0A%20%20%20%20T_max_pend%20%20%20%20%20%20%3D%203%20*%202%20*%20np.pi%20%2F%20omega_pend%20%20%20%20%23%20three%20full%20linearised%20periods%0A%0A%20%20%20%20def%20_rhs_pend(t%2C%20u)%3A%0A%20%20%20%20%20%20%20%20return%20%5Bu%5B1%5D%2C%20-(_g_p%20%2F%20_l_p)%20*%20np.sin(u%5B0%5D)%5D%0A%0A%20%20%20%20sol_pend%20%3D%20solve_ivp(%0A%20%20%20%20%20%20%20%20_rhs_pend%2C%20%5B0%2C%20T_max_pend%5D%2C%20%5Btheta0_rad_pend%2C%200.0%5D%2C%0A%20%20%20%20%20%20%20%20dense_output%3DTrue%2C%20max_step%3D0.01%2C%0A%20%20%20%20)%0A%20%20%20%20g_pend%20%3D%20_g_p%0A%20%20%20%20l_pend%20%3D%20_l_p%0A%20%20%20%20return%20T_max_pend%2C%20l_pend%2C%20omega_pend%2C%20sol_pend%2C%20theta0_rad_pend%0A%0A%0A%40app.cell%0Adef%20_(T_max_pend%2C%20np%2C%20omega_pend%2C%20plt%2C%20sol_pend%2C%20theta0_rad_pend)%3A%0A%20%20%20%20_t_plot%20%3D%20np.linspace(0%2C%20T_max_pend%2C%20800)%0A%20%20%20%20_theta_lin%20%3D%20theta0_rad_pend%20*%20np.cos(omega_pend%20*%20_t_plot)%0A%20%20%20%20_theta_nonlin%20%3D%20sol_pend.sol(_t_plot)%5B0%5D%0A%0A%20%20%20%20_fig1%2C%20_ax1%20%3D%20plt.subplots(figsize%3D(10%2C%204))%0A%20%20%20%20_ax1.plot(_t_plot%2C%20np.rad2deg(_theta_lin)%2C%20color%3D%22%231a5276%22%2C%20linewidth%3D2%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20label%3D%22linearised%20%20%24%5C%5Ctheta_0%5C%5Ccos(%5C%5Comega%20t)%24%22)%0A%20%20%20%20_ax1.plot(_t_plot%2C%20np.rad2deg(_theta_nonlin)%2C%20color%3D%22%23922b21%22%2C%20linewidth%3D2%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20linestyle%3D%22--%22%2C%20label%3D%22nonlinear%20(solve_ivp)%22)%0A%20%20%20%20_ax1.axhline(0%2C%20color%3D%22grey%22%2C%20linewidth%3D0.8%2C%20alpha%3D0.5)%0A%20%20%20%20_ax1.set_xlabel(r%22%24t%24%20(s)%22%2C%20fontsize%3D13)%0A%20%20%20%20_ax1.set_ylabel(r%22%24%5Ctheta%24%20(degrees)%22%2C%20fontsize%3D13)%0A%20%20%20%20_ax1.set_title(%0A%20%20%20%20%20%20%20%20rf%22Linearised%20vs%20nonlinear%20pendulum%20%20%20(%24%5Ctheta_0%20%3D%20%7Bnp.rad2deg(theta0_rad_pend)%3A.0f%7D%C2%B0%24)%22%2C%0A%20%20%20%20%20%20%20%20fontsize%3D13%2C%0A%20%20%20%20)%0A%20%20%20%20_ax1.legend(fontsize%3D11)%0A%20%20%20%20_ax1.grid(True%2C%20alpha%3D0.3)%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%20Animation%0A%0A%20%20%20%20The%20time%20slider%20below%20drives%20a%20side-by-side%20drawing%20of%20both%20pendulums.%0A%20%20%20%20Watch%20the%20nonlinear%20pendulum%20(red)%20fall%20behind%20the%20linearised%20one%20(blue)%0A%20%20%20%20as%20the%20amplitude%20increases%20%E2%80%94%20the%20nonlinear%20period%20is%20longer.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(T_max_pend%2C%20mo)%3A%0A%20%20%20%20t_anim_slider%20%3D%20mo.ui.slider(0.0%2C%20T_max_pend%2C%20step%3D0.05%2C%20value%3D0.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%20label%3Dr%22%24t%24%20(s)%22)%0A%20%20%20%20t_anim_slider%0A%20%20%20%20return%20(t_anim_slider%2C)%0A%0A%0A%40app.cell%0Adef%20_(l_pend%2C%20np%2C%20omega_pend%2C%20plt%2C%20sol_pend%2C%20t_anim_slider%2C%20theta0_rad_pend)%3A%0A%20%20%20%20_t_now%20%3D%20t_anim_slider.value%0A%20%20%20%20_th_lin%20%3D%20theta0_rad_pend%20*%20np.cos(omega_pend%20*%20_t_now)%0A%20%20%20%20_th_non%20%3D%20sol_pend.sol(_t_now)%5B0%5D%0A%0A%20%20%20%20_fig2%2C%20_axes2%20%3D%20plt.subplots(1%2C%202%2C%20figsize%3D(8%2C%205))%0A%0A%20%20%20%20for%20_ax2%2C%20_theta%2C%20_title%2C%20_color%20in%20%5B%0A%20%20%20%20%20%20%20%20(_axes2%5B0%5D%2C%20_th_lin%2C%20%20%22Linearised%22%2C%20%22%231a5276%22)%2C%0A%20%20%20%20%20%20%20%20(_axes2%5B1%5D%2C%20_th_non%2C%20%20%22Nonlinear%22%2C%20%20%22%23922b21%22)%2C%0A%20%20%20%20%5D%3A%0A%20%20%20%20%20%20%20%20_xb%20%3D%20l_pend%20*%20np.sin(_theta)%0A%20%20%20%20%20%20%20%20_yb%20%3D%20-l_pend%20*%20np.cos(_theta)%0A%20%20%20%20%20%20%20%20%23%20pivot%20bar%0A%20%20%20%20%20%20%20%20_ax2.plot(%5B-0.25%2C%200.25%5D%2C%20%5B0.04%2C%200.04%5D%2C%20color%3D%22%23333%22%2C%20linewidth%3D5%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20solid_capstyle%3D%22round%22)%0A%20%20%20%20%20%20%20%20%23%20rod%0A%20%20%20%20%20%20%20%20_ax2.plot(%5B0%2C%20_xb%5D%2C%20%5B0%2C%20_yb%5D%2C%20color%3D_color%2C%20linewidth%3D3.5%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20solid_capstyle%3D%22round%22)%0A%20%20%20%20%20%20%20%20%23%20bob%0A%20%20%20%20%20%20%20%20_circle%20%3D%20plt.Circle((_xb%2C%20_yb)%2C%200.09%2C%20color%3D_color%2C%20zorder%3D5)%0A%20%20%20%20%20%20%20%20_ax2.add_patch(_circle)%0A%20%20%20%20%20%20%20%20%23%20pivot%20pin%0A%20%20%20%20%20%20%20%20_ax2.plot(0%2C%200%2C%20%22o%22%2C%20color%3D%22%23333%22%2C%20markersize%3D7%2C%20zorder%3D6)%0A%0A%20%20%20%20%20%20%20%20_ax2.set_xlim(-1.3%2C%201.3)%0A%20%20%20%20%20%20%20%20_ax2.set_ylim(-1.3%2C%200.2)%0A%20%20%20%20%20%20%20%20_ax2.set_aspect(%22equal%22)%0A%20%20%20%20%20%20%20%20_ax2.set_title(_title%2C%20fontsize%3D13%2C%20color%3D_color%2C%20fontweight%3D%22bold%22)%0A%20%20%20%20%20%20%20%20_ax2.axis(%22off%22)%0A%0A%20%20%20%20_fig2.suptitle(rf%22%24t%20%3D%20%7B_t_now%3A.2f%7D%24%20s%22%2C%20fontsize%3D12%2C%20color%3D%22%23555%22)%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%20The%20period%20is%20not%20constant%0A%0A%20%20%20%20For%20a%20linear%20oscillator%20the%20period%20%24T%20%3D%202%5Cpi%2F%5Comega%24%20is%20the%20same%20regardless%20of%0A%20%20%20%20amplitude.%20%20For%20the%20nonlinear%20pendulum%20the%20period%20grows%20with%20amplitude%20%E2%80%94%20the%0A%20%20%20%20larger%20the%20swing%2C%20the%20longer%20it%20takes%20to%20complete%20one%20oscillation.%0A%0A%20%20%20%20This%20is%20impossible%20to%20see%20from%20the%20linearised%20formula.%20%20It%20is%20immediately%20visible%0A%20%20%20%20here%3A%20the%20red%20pendulum%20lags%20further%20and%20further%20behind%20the%20blue%20one%20as%20%24%5Ctheta_0%24%0A%20%20%20%20increases.%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-%20The%20exact%20pendulum%20ODE%20is%20nonlinear.%20%20The%20linearised%20ODE%20(%24%5Csin%5Ctheta%20%5Capprox%20%5Ctheta%24)%0A%20%20%20%20%20%20has%20an%20analytic%20solution%3B%20the%20exact%20one%20requires%20a%20numerical%20solver.%0A%20%20%20%20-%20Rewrite%20the%20second-order%20ODE%20as%20a%20two-component%20system%20before%20calling%20%60solve_ivp%60.%0A%20%20%20%20-%20The%20linearised%20solution%20is%20a%20good%20approximation%20for%20small%20angles.%20%20At%20large%0A%20%20%20%20%20%20angles%20the%20period%20differs%20and%20the%20curves%20drift%20apart.%0A%20%20%20%20-%20The%20small-angle%20approximation%20is%20not%20a%20curiosity%20%E2%80%94%20it%20is%20the%20default%20for%20most%0A%20%20%20%20%20%20pen-and-paper%20pendulum%20problems.%20%20Computation%20lets%20us%20see%20exactly%20where%20it%20fails.%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