{"id":193,"date":"2019-01-10T00:48:26","date_gmt":"2019-01-10T08:48:26","guid":{"rendered":"http:\/\/gantovnik.com\/bio-tips\/?p=193"},"modified":"2024-04-26T17:50:03","modified_gmt":"2024-04-27T00:50:03","slug":"double-pendulum","status":"publish","type":"post","link":"https:\/\/gantovnik.com\/bio-tips\/2019\/01\/double-pendulum\/","title":{"rendered":"#39 Double pendulum using python"},"content":{"rendered":"

import os
\nimport numpy as np
\nimport matplotlib.pyplot as plt
\nfrom scipy import integrate
\nimport sympy
\nos.chdir(r’D:\\projects\\wordpress\\ex39′)
\nos.getcwd()
\nt, g, m1, l1, m2, l2 = sympy.symbols(“t, g, m_1, l_1, m_2, l_2”)
\ntheta1, theta2 = sympy.symbols(“theta_1, theta_2”, cls=sympy.Function)
\node1 = sympy.Eq((m1+m2)*l1 * theta1(t).diff(t,t) +
\nm2*l2 * theta2(t).diff(t,t) +
\nm2*l2 * theta2(t).diff(t)**2 * sympy.sin(theta1(t)-theta2(t)) +
\ng*(m1+m2) * sympy.sin(theta1(t)))
\node1
\node2 = sympy.Eq(m2*l2 * theta2(t).diff(t,t) +
\nm2*l1 * theta1(t).diff(t,t) * sympy.cos(theta1(t)-theta2(t)) –
\nm2*l1 * theta1(t).diff(t)**2 * sympy.sin(theta1(t) – theta2(t)) +
\nm2*g * sympy.sin(theta2(t)))
\node2
\n# sympy cannot solve these ODEs
\ntry:
\n sympy.dsolve(ode1, ode2)
\nexcept Exception as e:
\n print(e)<\/p>\n

y1, y2, y3, y4 = sympy.symbols(“y_1, y_2, y_3, y_4”, cls=sympy.Function)
\nvarchange = {theta1(t).diff(t, t): y2(t).diff(t),
\ntheta1(t): y1(t),
\ntheta2(t).diff(t, t): y4(t).diff(t),
\ntheta2(t): y3(t)}
\node1_vc = ode1.subs(varchange)
\node2_vc = ode2.subs(varchange)
\node3 = y1(t).diff(t) – y2(t)
\node4 = y3(t).diff(t) – y4(t)
\ny = sympy.Matrix([y1(t), y2(t), y3(t), y4(t)])
\nvcsol = sympy.solve((ode1_vc, ode2_vc, ode3, ode4), y.diff(t), dict=True)
\nf = y.diff(t).subs(vcsol[0])
\nsympy.Eq(y.diff(t), f)
\nparams = {m1: 5.0, l1: 2.0,
\nm2: 1.0, l2: 1.0, g: 10.0}
\nf_np = sympy.lambdify((t, y), f.subs(params), ‘numpy’)
\njac = sympy.Matrix([[fj.diff(yi) for yi in y] for fj in f])
\njac_np = sympy.lambdify((t, y), jac.subs(params), ‘numpy’)
\ny0 = [2.0, 0, 0.0, 0]\nt = np.linspace(0, 20, 1000)
\njac_np(0, y0)
\nr = integrate.ode(f_np, jac_np).set_initial_value(y0, t[0]);
\ndt = t[1] – t[0]\ny = np.zeros((len(t), len(y0)))
\nidx = 0
\nwhile r.successful() and r.t < t[-1]:
\n y[idx, :] = r.y
\n r.integrate(r.t + dt)
\n idx += 1<\/p>\n

fig = plt.figure(figsize=(10, 4))
\nax1 = plt.subplot2grid((2, 5), (0, 0), colspan=3)
\nax2 = plt.subplot2grid((2, 5), (1, 0), colspan=3)
\nax3 = plt.subplot2grid((2, 5), (0, 3), colspan=2, rowspan=2)
\nax1.plot(t, y[:, 0], ‘r’)
\nax1.set_ylabel(r’$\\theta_1$’, fontsize=18)
\nax2.plot(t, y[:, 2], ‘b’)
\nax2.set_xlabel(‘$t$’, fontsize=18)
\nax2.set_ylabel(r’$\\theta_2$’, fontsize=18)
\nax3.plot(y[:, 0], y[:, 2], ‘k’)
\nax3.set_xlabel(r’$\\theta_1$’, fontsize=18)
\nax3.set_ylabel(r’$\\theta_2$’, fontsize=18)
\nfig.tight_layout()
\nplt.savefig(“example39.png”, dpi=100)<\/p>\n

\"example39\"<\/p>\n

theta1_np, theta2_np = y[:, 0], y[:, 2]\nx1 = params[l1] * np.sin(theta1_np)
\ny1 = -params[l1] * np.cos(theta1_np)
\nx2 = x1 + params[l2] * np.sin(theta2_np)
\ny2 = y1 – params[l2] * np.cos(theta2_np)
\nfig = plt.figure(figsize=(10, 4))
\nax1 = plt.subplot2grid((2, 5), (0, 0), colspan=3)
\nax2 = plt.subplot2grid((2, 5), (1, 0), colspan=3)
\nax3 = plt.subplot2grid((2, 5), (0, 3), colspan=2, rowspan=2)
\nax1.plot(t, x1, ‘r’)
\nax1.plot(t, y1, ‘b’)
\nax1.set_ylabel(‘$x_1, y_1$’, fontsize=18)
\nax1.set_yticks([-3, 0, 3])
\nax2.plot(t, x2, ‘r’)
\nax2.plot(t, y2, ‘b’)
\nax2.set_xlabel(‘$t$’, fontsize=18)
\nax2.set_ylabel(‘$x_2, y_2$’, fontsize=18)
\nax2.set_yticks([-3, 0, 3])
\nax3.plot(x1, y1, ‘r’)
\nax3.plot(x2, y2, ‘b’, lw=0.5)
\nax3.set_xlabel(‘$x$’, fontsize=18)
\nax3.set_ylabel(‘$y$’, fontsize=18)
\nax3.set_xticks([-3, 0, 3])
\nax3.set_yticks([-3, 0, 3])
\nfig.tight_layout()
\nplt.savefig(“example39_2.png”, dpi=100)
\nplt.show()
\nplt.close()<\/p>\n

\"example39_2\"<\/pre>\n","protected":false},"excerpt":{"rendered":"

import os import numpy as np import matplotlib.pyplot as plt from scipy import integrate import sympy os.chdir(r’D:\\projects\\wordpress\\ex39′) os.getcwd() t, g, m1, l1, m2, l2 = sympy.symbols(“t, g, m_1, l_1, m_2, l_2”) theta1, theta2 = sympy.symbols(“theta_1, theta_2”, cls=sympy.Function) ode1 = sympy.Eq((m1+m2)*l1 * theta1(t).diff(t,t) + m2*l2 * theta2(t).diff(t,t) + m2*l2 * theta2(t).diff(t)**2 * sympy.sin(theta1(t)-theta2(t)) + g*(m1+m2) * […]<\/p>\n","protected":false},"author":1,"featured_media":0,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"nf_dc_page":"","_et_pb_use_builder":"","_et_pb_old_content":"","_et_gb_content_width":"","_lmt_disableupdate":"yes","_lmt_disable":"","jetpack_post_was_ever_published":false,"_jetpack_newsletter_access":"","_jetpack_dont_email_post_to_subs":false,"_jetpack_newsletter_tier_id":0,"_jetpack_memberships_contains_paywalled_content":false,"_jetpack_memberships_contains_paid_content":false,"footnotes":""},"categories":[2],"tags":[],"class_list":["post-193","post","type-post","status-publish","format-standard","hentry","category-python"],"modified_by":"gantovnik","jetpack_featured_media_url":"","jetpack_sharing_enabled":true,"jetpack_shortlink":"https:\/\/wp.me\/p8bH0k-37","jetpack_likes_enabled":true,"jetpack-related-posts":[{"id":148,"url":"https:\/\/gantovnik.com\/bio-tips\/2019\/01\/damped-harmonic-oscillator\/","url_meta":{"origin":193,"position":0},"title":"Damped harmonic oscillator","author":"gantovnik","date":"2019-01-06","format":false,"excerpt":"import os import numpy as np import matplotlib.pyplot as plt import matplotlib as mpl mpl.rcParams['text.usetex'] = True import sympy from IPython.display import display sympy.init_printing() def apply_ics(sol, ics, x, known_params): free_params = sol.free_symbols - set(known_params) eqs = [(sol.lhs.diff(x, n) - sol.rhs.diff(x, n)).subs(x, 0).subs(ics) for n in range(len(ics))] sol_params = sympy.solve(eqs, free_params)\u2026","rel":"","context":"In "python"","block_context":{"text":"python","link":"https:\/\/gantovnik.com\/bio-tips\/category\/python\/"},"img":{"alt_text":"example31","src":"https:\/\/i0.wp.com\/gantovnik.com\/bio-tips\/wp-content\/uploads\/2019\/01\/example31.png?resize=350%2C200","width":350,"height":200,"srcset":"https:\/\/i0.wp.com\/gantovnik.com\/bio-tips\/wp-content\/uploads\/2019\/01\/example31.png?resize=350%2C200 1x, https:\/\/i0.wp.com\/gantovnik.com\/bio-tips\/wp-content\/uploads\/2019\/01\/example31.png?resize=525%2C300 1.5x"},"classes":[]},{"id":1204,"url":"https:\/\/gantovnik.com\/bio-tips\/2021\/11\/210-parametric-curve-in-3d-2-2\/","url_meta":{"origin":193,"position":1},"title":"#212 The double pendulum simulation using python","author":"gantovnik","date":"2021-11-27","format":false,"excerpt":"","rel":"","context":"In "animation"","block_context":{"text":"animation","link":"https:\/\/gantovnik.com\/bio-tips\/category\/animation\/"},"img":{"alt_text":"","src":"https:\/\/i0.wp.com\/gantovnik.com\/bio-tips\/wp-content\/uploads\/2021\/11\/ex212.png?resize=350%2C200&ssl=1","width":350,"height":200},"classes":[]},{"id":143,"url":"https:\/\/gantovnik.com\/bio-tips\/2019\/01\/143\/","url_meta":{"origin":193,"position":2},"title":"#30 Newton’s Law of Cooling","author":"gantovnik","date":"2019-01-06","format":false,"excerpt":"import os import numpy as np import matplotlib.pyplot as plt import matplotlib as mpl mpl.rcParams['text.usetex'] = True import sympy from IPython.display import display sympy.init_printing() #%matplotlib inline #%config InlineBackend.figure_format='retina' os.chdir(r'D:\\projects\\wordpress\\ex30') os.getcwd() #Symbolic ODE solving with SymPy #Newton's law of cooling t, k, T0, Ta = sympy.symbols(\"t, k, T_0, T_a\") T =\u2026","rel":"","context":"In "python"","block_context":{"text":"python","link":"https:\/\/gantovnik.com\/bio-tips\/category\/python\/"},"img":{"alt_text":"example30","src":"https:\/\/i0.wp.com\/gantovnik.com\/bio-tips\/wp-content\/uploads\/2019\/01\/example30.png?resize=350%2C200","width":350,"height":200,"srcset":"https:\/\/i0.wp.com\/gantovnik.com\/bio-tips\/wp-content\/uploads\/2019\/01\/example30.png?resize=350%2C200 1x, https:\/\/i0.wp.com\/gantovnik.com\/bio-tips\/wp-content\/uploads\/2019\/01\/example30.png?resize=525%2C300 1.5x"},"classes":[]},{"id":154,"url":"https:\/\/gantovnik.com\/bio-tips\/2019\/01\/inexact-solutions-to-odes\/","url_meta":{"origin":193,"position":3},"title":"Inexact solutions to ODEs","author":"gantovnik","date":"2019-01-08","format":false,"excerpt":"\u00a0 import os import numpy as np import matplotlib.pyplot as plt import matplotlib as mpl import sympy from IPython.display import display sympy.init_printing() mpl.rcParams['text.usetex'] = True import sympy os.chdir(r'D:\\projects\\wordpress\\ex33') os.getcwd() def plot_direction_field(x, y_x, f_xy, x_lim=(-5, 5), y_lim=(-5, 5), ax=None): f_np = sympy.lambdify((x, y_x), f_xy, 'numpy') x_vec = np.linspace(x_lim[0], x_lim[1], 20) y_vec\u2026","rel":"","context":"In "python"","block_context":{"text":"python","link":"https:\/\/gantovnik.com\/bio-tips\/category\/python\/"},"img":{"alt_text":"example33","src":"https:\/\/i0.wp.com\/gantovnik.com\/bio-tips\/wp-content\/uploads\/2019\/01\/example33.png?resize=350%2C200","width":350,"height":200,"srcset":"https:\/\/i0.wp.com\/gantovnik.com\/bio-tips\/wp-content\/uploads\/2019\/01\/example33.png?resize=350%2C200 1x, https:\/\/i0.wp.com\/gantovnik.com\/bio-tips\/wp-content\/uploads\/2019\/01\/example33.png?resize=525%2C300 1.5x"},"classes":[]},{"id":190,"url":"https:\/\/gantovnik.com\/bio-tips\/2019\/01\/coupled-damped-springs-jacobian-is-available\/","url_meta":{"origin":193,"position":4},"title":"Coupled damped springs (Jacobian is available)","author":"gantovnik","date":"2019-01-10","format":false,"excerpt":"\u00a0 import os import numpy as np import matplotlib.pyplot as plt from scipy import integrate os.chdir(r'D:\\projects\\wordpress\\ex38') os.getcwd() def f(t, y, args): m1, k1, g1, m2, k2, g2 = args return [y[1], - k1\/m1 * y[0] + k2\/m1 * (y[2] - y[0]) - g1\/m1 * y[1], y[3], - k2\/m2 * (y[2]\u2026","rel":"","context":"In "python"","block_context":{"text":"python","link":"https:\/\/gantovnik.com\/bio-tips\/category\/python\/"},"img":{"alt_text":"example38","src":"https:\/\/i0.wp.com\/gantovnik.com\/bio-tips\/wp-content\/uploads\/2019\/01\/example38.png?resize=350%2C200","width":350,"height":200,"srcset":"https:\/\/i0.wp.com\/gantovnik.com\/bio-tips\/wp-content\/uploads\/2019\/01\/example38.png?resize=350%2C200 1x, https:\/\/i0.wp.com\/gantovnik.com\/bio-tips\/wp-content\/uploads\/2019\/01\/example38.png?resize=525%2C300 1.5x"},"classes":[]},{"id":151,"url":"https:\/\/gantovnik.com\/bio-tips\/2019\/01\/direction-fields\/","url_meta":{"origin":193,"position":5},"title":"Direction fields","author":"gantovnik","date":"2019-01-07","format":false,"excerpt":"","rel":"","context":"In "python"","block_context":{"text":"python","link":"https:\/\/gantovnik.com\/bio-tips\/category\/python\/"},"img":{"alt_text":"example32","src":"https:\/\/i0.wp.com\/gantovnik.com\/bio-tips\/wp-content\/uploads\/2019\/01\/example32.png?resize=350%2C200","width":350,"height":200,"srcset":"https:\/\/i0.wp.com\/gantovnik.com\/bio-tips\/wp-content\/uploads\/2019\/01\/example32.png?resize=350%2C200 1x, https:\/\/i0.wp.com\/gantovnik.com\/bio-tips\/wp-content\/uploads\/2019\/01\/example32.png?resize=525%2C300 1.5x"},"classes":[]}],"_links":{"self":[{"href":"https:\/\/gantovnik.com\/bio-tips\/wp-json\/wp\/v2\/posts\/193","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/gantovnik.com\/bio-tips\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/gantovnik.com\/bio-tips\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/gantovnik.com\/bio-tips\/wp-json\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/gantovnik.com\/bio-tips\/wp-json\/wp\/v2\/comments?post=193"}],"version-history":[{"count":2,"href":"https:\/\/gantovnik.com\/bio-tips\/wp-json\/wp\/v2\/posts\/193\/revisions"}],"predecessor-version":[{"id":2151,"href":"https:\/\/gantovnik.com\/bio-tips\/wp-json\/wp\/v2\/posts\/193\/revisions\/2151"}],"wp:attachment":[{"href":"https:\/\/gantovnik.com\/bio-tips\/wp-json\/wp\/v2\/media?parent=193"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/gantovnik.com\/bio-tips\/wp-json\/wp\/v2\/categories?post=193"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/gantovnik.com\/bio-tips\/wp-json\/wp\/v2\/tags?post=193"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}