obsidian_backup/补课/多体动力学/18_Simulation_and_Visualization.ipynb

1375 lines
102 KiB
Plaintext
Raw Normal View History

2025-03-28 09:27:28 +08:00
{
"cells": [
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"import sympy as sm\n",
"import sympy.physics.mechanics as me\n",
"me.init_vprinting(use_latex='mathjax')"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\left( \\frac{l u_{1}}{2}\\mathbf{\\hat{a}_y}, \\ l u_{1}\\mathbf{\\hat{a}_y}, \\ u_{1}\\mathbf{\\hat{n}_z}, \\ u_{2}\\mathbf{\\hat{a}_x} + u_{1}\\mathbf{\\hat{n}_z}\\right)$"
],
"text/plain": [
"⎛l⋅u₁ ⎞\n",
"⎜──── a_y, l⋅u₁ a_y, u₁ n_z, u₂ a_x + u₁ n_z⎟\n",
"⎝ 2 ⎠"
]
},
"execution_count": 5,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"m, g, kt, kl, l = sm.symbols('m, g, k_t, k_l, l')\n",
"q1, q2, q3 = me.dynamicsymbols('q1, q2, q3')\n",
"u1, u2, u3 = me.dynamicsymbols('u1, u2, u3')\n",
"\n",
"N = me.ReferenceFrame('N')\n",
"A = me.ReferenceFrame('A')\n",
"B = me.ReferenceFrame('B')\n",
"\n",
"A.orient_axis(N, q1, N.z)\n",
"B.orient_axis(A, q2, A.x)\n",
"\n",
"A.set_ang_vel(N, u1*N.z)\n",
"B.set_ang_vel(A, u2*A.x)\n",
"\n",
"O = me.Point('O')\n",
"Ao = me.Point('A_O')\n",
"Bo = me.Point('B_O')\n",
"\n",
"Ao.set_pos(O, l/2*A.x)\n",
"Bo.set_pos(O, l*A.x)\n",
"\n",
"O.set_vel(N, 0)\n",
"Ao.v2pt_theory(O, N, A)\n",
"Bo.v2pt_theory(O, N, A)\n",
"\n",
"Ao.vel(N), Bo.vel(N), A.ang_vel_in(N), B.ang_vel_in(N)"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle - q_{3} u_{1} \\cos{\\left(q_{2} \\right)}\\mathbf{\\hat{b}_x} + u_{3}\\mathbf{\\hat{b}_y} + q_{3} u_{2}\\mathbf{\\hat{b}_z} + l u_{1}\\mathbf{\\hat{a}_y}$"
],
"text/plain": [
"-q₃⋅u₁⋅cos(q₂) b_x + u₃ b_y + q₃⋅u₂ b_z + l⋅u₁ a_y"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Q = me.Point('Q')\n",
"Q.set_pos(Bo, q3*B.y)\n",
"Q.set_vel(B, u3*B.y)\n",
"Q.v1pt_theory(Bo, N, B)\n",
"\n",
"Q.vel(N)"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\left( - \\frac{l u_{1}^{2}}{2}\\mathbf{\\hat{a}_x} + \\frac{l \\dot{u}_{1}}{2}\\mathbf{\\hat{a}_y}, \\ - l u_{1}^{2}\\mathbf{\\hat{a}_x} + l \\dot{u}_{1}\\mathbf{\\hat{a}_y}, \\ \\dot{u}_{1}\\mathbf{\\hat{n}_z}, \\ \\dot{u}_{2}\\mathbf{\\hat{a}_x} + u_{1} u_{2}\\mathbf{\\hat{a}_y} + \\dot{u}_{1}\\mathbf{\\hat{n}_z}\\right)$"
],
"text/plain": [
"⎛ 2 ↪\n",
"⎜-l⋅u₁ l⋅u₁̇ 2 ↪\n",
"⎜─────── a_x + ──── a_y, -l⋅u₁ a_x + l⋅u₁̇ a_y, u₁̇ n_z, u₂̇ a_x + u₁⋅u₂ a_y + u ↪\n",
"⎝ 2 2 ↪\n",
"\n",
"↪ ⎞\n",
"↪ ⎟\n",
"↪ ₁̇ n_z⎟\n",
"↪ ⎠"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Ao.acc(N), Bo.acc(N), A.ang_acc_in(N), B.ang_acc_in(N)"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle (q_{3} u_{1} u_{2} \\sin{\\left(q_{2} \\right)} + q_{3} u_{1} \\sin{\\left(q_{2} \\right)} \\dot{q}_{2} - q_{3} \\cos{\\left(q_{2} \\right)} \\dot{u}_{1} - u_{1} u_{3} \\cos{\\left(q_{2} \\right)} - u_{1} \\cos{\\left(q_{2} \\right)} \\dot{q}_{3})\\mathbf{\\hat{b}_x} + (- q_{3} u_{1}^{2} \\cos^{2}{\\left(q_{2} \\right)} - q_{3} u_{2}^{2} + \\dot{u}_{3})\\mathbf{\\hat{b}_y} + (q_{3} u_{1}^{2} \\sin{\\left(q_{2} \\right)} \\cos{\\left(q_{2} \\right)} + q_{3} \\dot{u}_{2} + u_{2} u_{3} + u_{2} \\dot{q}_{3})\\mathbf{\\hat{b}_z} - l u_{1}^{2}\\mathbf{\\hat{a}_x} + l \\dot{u}_{1}\\mathbf{\\hat{a}_y}$"
],
"text/plain": [
" ↪\n",
"(q₃⋅u₁⋅u₂⋅sin(q₂) + q₃⋅u₁⋅sin(q₂)⋅q₂̇ - q₃⋅cos(q₂)⋅u₁̇ - u₁⋅u₃⋅cos(q₂) - u₁⋅cos( ↪\n",
"\n",
"↪ ⎛ 2 2 2 ⎞ ⎛ 2 ↪\n",
"↪ q₂)⋅q₃̇) b_x + ⎝- q₃⋅u₁ ⋅cos (q₂) - q₃⋅u₂ + u₃̇⎠ b_y + ⎝q₃⋅u₁ ⋅sin(q₂)⋅cos(q₂ ↪\n",
"\n",
"↪ ⎞ 2 \n",
"↪ ) + q₃⋅u₂̇ + u₂⋅u₃ + u₂⋅q₃̇⎠ b_z + -l⋅u₁ a_x + l⋅u₁̇ a_y"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Q.acc(N)"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle (2 q_{3} u_{1} u_{2} \\sin{\\left(q_{2} \\right)} - q_{3} \\cos{\\left(q_{2} \\right)} \\dot{u}_{1} - 2 u_{1} u_{3} \\cos{\\left(q_{2} \\right)})\\mathbf{\\hat{b}_x} + (- q_{3} u_{1}^{2} \\cos^{2}{\\left(q_{2} \\right)} - q_{3} u_{2}^{2} + \\dot{u}_{3})\\mathbf{\\hat{b}_y} + (q_{3} u_{1}^{2} \\sin{\\left(q_{2} \\right)} \\cos{\\left(q_{2} \\right)} + q_{3} \\dot{u}_{2} + 2 u_{2} u_{3})\\mathbf{\\hat{b}_z} - l u_{1}^{2}\\mathbf{\\hat{a}_x} + l \\dot{u}_{1}\\mathbf{\\hat{a}_y}$"
],
"text/plain": [
" ⎛ 2 2 ↪\n",
"(2⋅q₃⋅u₁⋅u₂⋅sin(q₂) - q₃⋅cos(q₂)⋅u₁̇ - 2⋅u₁⋅u₃⋅cos(q₂)) b_x + ⎝- q₃⋅u₁ ⋅cos (q₂ ↪\n",
"\n",
"↪ 2 ⎞ ⎛ 2 ⎞ ↪\n",
"↪ ) - q₃⋅u₂ + u₃̇⎠ b_y + ⎝q₃⋅u₁ ⋅sin(q₂)⋅cos(q₂) + q₃⋅u₂̇ + 2⋅u₂⋅u₃⎠ b_z + -l⋅u ↪\n",
"\n",
"↪ 2 \n",
"↪ ₁ a_x + l⋅u₁̇ a_y"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"t = me.dynamicsymbols._t\n",
"\n",
"qdot_repl = {q1.diff(t): u1,\n",
" q2.diff(t): u2,\n",
" q3.diff(t): u3}\n",
"\n",
"Q.set_acc(N, Q.acc(N).xreplace(qdot_repl))\n",
"Q.acc(N)"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"R_Ao = m*g*N.x\n",
"R_Bo = m*g*N.x + kl*q3*B.y\n",
"R_Q = m/4*g*N.x - kl*q3*B.y\n",
"T_A = -kt*q1*N.z + kt*q2*A.x\n",
"T_B = -kt*q2*A.x"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"I = m*l**2/12\n",
"I_A_Ao = I*me.outer(A.y, A.y) + I*me.outer(A.z, A.z)\n",
"I_B_Bo = I*me.outer(B.x, B.x) + I*me.outer(B.z, B.z)"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"v_Ao_1 = Ao.vel(N).diff(u1, N)\n",
"v_Bo_1 = Bo.vel(N).diff(u1, N)\n",
"v_Q_1 = Q.vel(N).diff(u1, N)\n",
"\n",
"v_Ao_2 = Ao.vel(N).diff(u2, N)\n",
"v_Bo_2 = Bo.vel(N).diff(u2, N)\n",
"v_Q_2 = Q.vel(N).diff(u2, N)\n",
"\n",
"v_Ao_3 = Ao.vel(N).diff(u3, N)\n",
"v_Bo_3 = Bo.vel(N).diff(u3, N)\n",
"v_Q_3 = Q.vel(N).diff(u3, N)\n",
"\n",
"w_A_1 = A.ang_vel_in(N).diff(u1, N)\n",
"w_B_1 = B.ang_vel_in(N).diff(u1, N)\n",
"\n",
"w_A_2 = A.ang_vel_in(N).diff(u2, N)\n",
"w_B_2 = B.ang_vel_in(N).diff(u2, N)\n",
"\n",
"w_A_3 = A.ang_vel_in(N).diff(u3, N)\n",
"w_B_3 = B.ang_vel_in(N).diff(u3, N)"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [],
"source": [
"F1 = v_Ao_1.dot(R_Ao) + v_Bo_1.dot(R_Bo) + v_Q_1.dot(R_Q) + w_A_1.dot(T_A) + w_B_1.dot(T_B)\n",
"F2 = v_Ao_2.dot(R_Ao) + v_Bo_2.dot(R_Bo) + v_Q_2.dot(R_Q) + w_A_2.dot(T_A) + w_B_2.dot(T_B)\n",
"F3 = v_Ao_3.dot(R_Ao) + v_Bo_3.dot(R_Bo) + v_Q_3.dot(R_Q) + w_A_3.dot(T_A) + w_B_3.dot(T_B)"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\left[\\begin{matrix}- \\frac{7 g l m \\sin{\\left(q_{1} \\right)}}{4} - \\frac{g m q_{3} \\cos{\\left(q_{1} \\right)} \\cos{\\left(q_{2} \\right)}}{4} - k_{t} q_{1}\\\\\\frac{g m q_{3} \\sin{\\left(q_{1} \\right)} \\sin{\\left(q_{2} \\right)}}{4} - k_{t} q_{2}\\\\- \\frac{g m \\sin{\\left(q_{1} \\right)} \\cos{\\left(q_{2} \\right)}}{4} - k_{l} q_{3}\\end{matrix}\\right]$"
],
"text/plain": [
"⎡ 7⋅g⋅l⋅m⋅sin(q₁) g⋅m⋅q₃⋅cos(q₁)⋅cos(q₂) ⎤\n",
"⎢- ─────────────── - ────────────────────── - kₜ⋅q₁⎥\n",
"⎢ 4 4 ⎥\n",
"⎢ ⎥\n",
"⎢ g⋅m⋅q₃⋅sin(q₁)⋅sin(q₂) ⎥\n",
"⎢ ────────────────────── - kₜ⋅q₂ ⎥\n",
"⎢ 4 ⎥\n",
"⎢ ⎥\n",
"⎢ g⋅m⋅sin(q₁)⋅cos(q₂) ⎥\n",
"⎢ - ─────────────────── - kₗ⋅q₃ ⎥\n",
"⎣ 4 ⎦"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Fr = sm.Matrix([F1, F2, F3])\n",
"Fr"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
"TAs = -(A.ang_acc_in(N).dot(I_A_Ao) + me.cross(A.ang_vel_in(N), I_A_Ao).dot(A.ang_vel_in(N)))\n",
"TBs = -(B.ang_acc_in(N).dot(I_B_Bo) + me.cross(B.ang_vel_in(N), I_B_Bo).dot(B.ang_vel_in(N)))\n",
"\n",
"F1s = v_Ao_1.dot(-m*Ao.acc(N)) + v_Bo_1.dot(-m*Bo.acc(N)) + v_Q_1.dot(-m/4*Q.acc(N))\n",
"F1s += w_A_1.dot(TAs) + w_B_1.dot(TBs)\n",
"\n",
"F2s = v_Ao_2.dot(-m*Ao.acc(N)) + v_Bo_2.dot(-m*Bo.acc(N)) + v_Q_2.dot(-m/4*Q.acc(N))\n",
"F2s += w_A_2.dot(TAs) + w_B_2.dot(TBs)\n",
"\n",
"F3s = v_Ao_3.dot(-m*Ao.acc(N)) + v_Bo_3.dot(-m*Bo.acc(N)) + v_Q_3.dot(-m/4*Q.acc(N))\n",
"F3s += w_A_3.dot(TAs) + w_B_3.dot(TBs)"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\left[\\begin{matrix}- \\frac{19 l^{2} m \\dot{u}_{1}}{12} - \\frac{l m q_{3} u_{1}^{2} \\cos{\\left(q_{2} \\right)}}{4} + l \\left(- \\frac{m \\left(- q_{3} u_{1}^{2} \\cos^{2}{\\left(q_{2} \\right)} - q_{3} u_{2}^{2} + \\dot{u}_{3}\\right) \\cos{\\left(q_{2} \\right)}}{4} + \\frac{m \\left(q_{3} u_{1}^{2} \\sin{\\left(q_{2} \\right)} \\cos{\\left(q_{2} \\right)} + q_{3} \\dot{u}_{2} + 2 u_{2} u_{3}\\right) \\sin{\\left(q_{2} \\right)}}{4}\\right) + \\frac{m \\left(2 q_{3} u_{1} u_{2} \\sin{\\left(q_{2} \\right)} - q_{3} \\cos{\\left(q_{2} \\right)} \\dot{u}_{1} - 2 u_{1} u_{3} \\cos{\\left(q_{2} \\right)}\\right) q_{3} \\cos{\\left(q_{2} \\right)}}{4} + \\left(- \\frac{l^{2} m \\left(- u_{1} u_{2} \\sin{\\left(q_{2} \\right)} + \\cos{\\left(q_{2} \\right)} \\dot{u}_{1}\\right)}{12} + \\frac{l^{2} m u_{1} u_{2} \\sin{\\left(q_{2} \\right)}}{12}\\right) \\cos{\\left(q_{2} \\right)}\\\\- \\frac{l^{2} m u_{1}^{2} \\sin{\\left(q_{2} \\right)} \\cos{\\left(q_{2} \\right)}}{12} - \\frac{l^{2} m \\dot{u}_{2}}{12} + \\frac{l m q_{3} \\sin{\\left(q_{2} \\right)} \\dot{u}_{1}}{4} - \\frac{m \\left(q_{3} u_{1}^{2} \\sin{\\left(q_{2} \\right)} \\cos{\\left(q_{2} \\right)} + q_{3} \\dot{u}_{2} + 2 u_{2} u_{3}\\right) q_{3}}{4}\\\\- \\frac{l m \\cos{\\left(q_{2} \\right)} \\dot{u}_{1}}{4} - \\frac{m \\left(- q_{3} u_{1}^{2} \\cos^{2}{\\left(q_{2} \\right)} - q_{3} u_{2}^{2} + \\dot{u}_{3}\\right)}{4}\\end{matrix}\\right]$"
],
"text/plain": [
"⎡ 2 2 ⎛ ⎛ 2 2 2 ⎞ ↪\n",
"⎢ 19⋅l ⋅m⋅u₁̇ l⋅m⋅q₃⋅u₁ ⋅cos(q₂) ⎜ m⋅⎝- q₃⋅u₁ ⋅cos (q₂) - q₃⋅u₂ + u₃̇⎠⋅ ↪\n",
"⎢- ────────── - ────────────────── + l⋅⎜- ──────────────────────────────────── ↪\n",
"⎢ 12 4 ⎝ 4 ↪\n",
"⎢ ↪\n",
"⎢ ↪\n",
"⎢ ↪\n",
"⎢ ↪\n",
"⎢ ↪\n",
"⎢ ↪\n",
"⎢ ↪\n",
"⎢ ↪\n",
"⎢ ↪\n",
"⎣ ↪\n",
"\n",
"↪ ⎛ 2 ⎞ ⎞ ↪\n",
"↪ cos(q₂) m⋅⎝q₃⋅u₁ ⋅sin(q₂)⋅cos(q₂) + q₃⋅u₂̇ + 2⋅u₂⋅u₃⎠⋅sin(q₂)⎟ m⋅(2⋅q₃⋅u₁ ↪\n",
"↪ ─────── + ────────────────────────────────────────────────────⎟ + ────────── ↪\n",
"↪ 4 ⎠ ↪\n",
"↪ ↪\n",
"↪ 2 2 2 ⎛ ↪\n",
"↪ l ⋅m⋅u₁ ⋅sin(q₂)⋅cos(q₂) l ⋅m⋅u₂̇ l⋅m⋅q₃⋅sin(q₂)⋅u₁̇ m⋅⎝q₃⋅u ↪\n",
"↪ - ──────────────────────── - ─────── + ───────────────── - ─────── ↪\n",
"↪ 12 12 4 ↪\n",
"↪ ↪\n",
"↪ ⎛ 2 2 ↪\n",
"↪ l⋅m⋅cos(q₂)⋅u₁̇ m⋅⎝- q₃⋅u₁ ⋅cos (q₂) ↪\n",
"↪ - ────────────── - ───────────────────── ↪\n",
"↪ 4 4 ↪\n",
"\n",
"↪ ⎛ 2 ↪\n",
"↪ ⋅u₂⋅sin(q₂) - q₃⋅cos(q₂)⋅u₁̇ - 2⋅u₁⋅u₃⋅cos(q₂))⋅q₃⋅cos(q₂) ⎜ l ⋅m⋅(-u₁⋅u₂⋅ ↪\n",
"↪ ───────────────────────────────────────────────────────── + ⎜- ───────────── ↪\n",
"↪ 4 ⎝ ↪\n",
"↪ ↪\n",
"↪ 2 ⎞ ↪\n",
"↪ ₁ ⋅sin(q₂)⋅cos(q₂) + q₃⋅u₂̇ + 2⋅u₂⋅u₃⎠⋅q₃ ↪\n",
"↪ ──────────────────────────────────────── ↪\n",
"↪ 4 ↪\n",
"↪ ↪\n",
"↪ 2 ⎞ ↪\n",
"↪ - q₃⋅u₂ + u₃̇⎠ ↪\n",
"↪ ────────────── ↪\n",
"↪ ↪\n",
"\n",
"↪ 2 ⎞ ⎤\n",
"↪ sin(q₂) + cos(q₂)⋅u₁̇) l ⋅m⋅u₁⋅u₂⋅sin(q₂)⎟ ⎥\n",
"↪ ───────────────────── + ──────────────────⎟⋅cos(q₂)⎥\n",
"↪ 12 12 ⎠ ⎥\n",
"↪ ⎥\n",
"↪ ⎥\n",
"↪ ⎥\n",
"↪ ⎥\n",
"↪ ⎥\n",
"↪ ⎥\n",
"↪ ⎥\n",
"↪ ⎥\n",
"↪ ⎥\n",
"↪ ⎦"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Frs = sm.Matrix([F1s, F2s, F3s])\n",
"Frs"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\left[\\begin{matrix}\\dot{u}_{1}\\\\\\dot{u}_{2}\\\\\\dot{u}_{3}\\end{matrix}\\right]$"
],
"text/plain": [
"⎡u₁̇⎤\n",
"⎢ ⎥\n",
"⎢u₂̇⎥\n",
"⎢ ⎥\n",
"⎣u₃̇⎦"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"u = sm.Matrix([u1, u2, u3])\n",
"ud = u.diff(t)\n",
"ud"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\left[\\begin{matrix}- \\frac{l^{2} m \\cos^{2}{\\left(q_{2} \\right)}}{12} - \\frac{19 l^{2} m}{12} - \\frac{m q_{3}^{2} \\cos^{2}{\\left(q_{2} \\right)}}{4} & \\frac{l m q_{3} \\sin{\\left(q_{2} \\right)}}{4} & - \\frac{l m \\cos{\\left(q_{2} \\right)}}{4}\\\\\\frac{l m q_{3} \\sin{\\left(q_{2} \\right)}}{4} & - \\frac{l^{2} m}{12} - \\frac{m q_{3}^{2}}{4} & 0\\\\- \\frac{l m \\cos{\\left(q_{2} \\right)}}{4} & 0 & - \\frac{m}{4}\\end{matrix}\\right]$"
],
"text/plain": [
"⎡ 2 2 2 2 2 ⎤\n",
"⎢ l ⋅m⋅cos (q₂) 19⋅l ⋅m m⋅q₃ ⋅cos (q₂) l⋅m⋅q₃⋅sin(q₂) -l⋅m⋅cos(q₂) ⎥\n",
"⎢- ───────────── - ─────── - ────────────── ────────────── ─────────────⎥\n",
"⎢ 12 12 4 4 4 ⎥\n",
"⎢ ⎥\n",
"⎢ 2 2 ⎥\n",
"⎢ l⋅m⋅q₃⋅sin(q₂) l ⋅m m⋅q₃ ⎥\n",
"⎢ ────────────── - ──── - ───── 0 ⎥\n",
"⎢ 4 12 4 ⎥\n",
"⎢ ⎥\n",
"⎢ -l⋅m⋅cos(q₂) -m ⎥\n",
"⎢ ───────────── 0 ─── ⎥\n",
"⎣ 4 4 ⎦"
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Md = Frs.jacobian(ud)\n",
"Md"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\left\\{ \\dot{u}_{1} : 0, \\ \\dot{u}_{2} : 0, \\ \\dot{u}_{3} : 0\\right\\}$"
],
"text/plain": [
"{u₁̇: 0, u₂̇: 0, u₃̇: 0}"
]
},
"execution_count": 19,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"ud_zerod = {udr: 0 for udr in ud}\n",
"ud_zerod\n"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\left[\\begin{matrix}- \\frac{7 g l m \\sin{\\left(q_{1} \\right)}}{4} - \\frac{g m q_{3} \\cos{\\left(q_{1} \\right)} \\cos{\\left(q_{2} \\right)}}{4} - k_{t} q_{1} + \\frac{l^{2} m u_{1} u_{2} \\sin{\\left(q_{2} \\right)} \\cos{\\left(q_{2} \\right)}}{6} - \\frac{l m q_{3} u_{1}^{2} \\cos{\\left(q_{2} \\right)}}{4} + l \\left(- \\frac{m \\left(- q_{3} u_{1}^{2} \\cos^{2}{\\left(q_{2} \\right)} - q_{3} u_{2}^{2}\\right) \\cos{\\left(q_{2} \\right)}}{4} + \\frac{m \\left(q_{3} u_{1}^{2} \\sin{\\left(q_{2} \\right)} \\cos{\\left(q_{2} \\right)} + 2 u_{2} u_{3}\\right) \\sin{\\left(q_{2} \\right)}}{4}\\right) + \\frac{m \\left(2 q_{3} u_{1} u_{2} \\sin{\\left(q_{2} \\right)} - 2 u_{1} u_{3} \\cos{\\left(q_{2} \\right)}\\right) q_{3} \\cos{\\left(q_{2} \\right)}}{4}\\\\\\frac{g m q_{3} \\sin{\\left(q_{1} \\right)} \\sin{\\left(q_{2} \\right)}}{4} - k_{t} q_{2} - \\frac{l^{2} m u_{1}^{2} \\sin{\\left(q_{2} \\right)} \\cos{\\left(q_{2} \\right)}}{12} - \\frac{m \\left(q_{3} u_{1}^{2} \\sin{\\left(q_{2} \\right)} \\cos{\\left(q_{2} \\right)} + 2 u_{2} u_{3}\\right) q_{3}}{4}\\\\- \\frac{g m \\sin{\\left(q_{1} \\right)} \\cos{\\left(q_{2} \\right)}}{4} - k_{l} q_{3} - \\frac{m \\left(- q_{3} u_{1}^{2} \\cos^{2}{\\left(q_{2} \\right)} - q_{3} u_{2}^{2}\\right)}{4}\\end{matrix}\\right]$"
],
"text/plain": [
"⎡ 2 ↪\n",
"⎢ 7⋅g⋅l⋅m⋅sin(q₁) g⋅m⋅q₃⋅cos(q₁)⋅cos(q₂) l ⋅m⋅u₁⋅u₂⋅sin(q₂)⋅cos(q ↪\n",
"⎢- ─────────────── - ────────────────────── - kₜ⋅q₁ + ──────────────────────── ↪\n",
"⎢ 4 4 6 ↪\n",
"⎢ ↪\n",
"⎢ ↪\n",
"⎢ g⋅m ↪\n",
"⎢ ─── ↪\n",
"⎢ ↪\n",
"⎢ ↪\n",
"⎢ ↪\n",
"⎢ ↪\n",
"⎢ ↪\n",
"⎣ ↪\n",
"\n",
"↪ 2 ⎛ ⎛ 2 2 2⎞ ⎛ ↪\n",
"↪ ₂) l⋅m⋅q₃⋅u₁ ⋅cos(q₂) ⎜ m⋅⎝- q₃⋅u₁ ⋅cos (q₂) - q₃⋅u₂ ⎠⋅cos(q₂) m⋅⎝q ↪\n",
"↪ ── - ────────────────── + l⋅⎜- ────────────────────────────────────── + ──── ↪\n",
"↪ 4 ⎝ 4 ↪\n",
"↪ ↪\n",
"↪ 2 2 ⎛ 2 ↪\n",
"↪ ⋅q₃⋅sin(q₁)⋅sin(q₂) l ⋅m⋅u₁ ⋅sin(q₂)⋅cos(q₂) m⋅⎝q₃⋅u₁ ⋅sin(q₂)⋅c ↪\n",
"↪ ─────────────────── - kₜ⋅q₂ - ──────────────────────── - ─────────────────── ↪\n",
"↪ 4 12 ↪\n",
"↪ ↪\n",
"↪ ⎛ 2 2 ↪\n",
"↪ g⋅m⋅sin(q₁)⋅cos(q₂) m⋅⎝- q₃⋅u₁ ⋅cos (q₂) - q₃⋅u₂ ↪\n",
"↪ - ─────────────────── - kₗ⋅q₃ - ──────────────────────────── ↪\n",
"↪ 4 4 ↪\n",
"\n",
"↪ 2 ⎞ ⎞ ↪\n",
"↪ ₃⋅u₁ ⋅sin(q₂)⋅cos(q₂) + 2⋅u₂⋅u₃⎠⋅sin(q₂)⎟ m⋅(2⋅q₃⋅u₁⋅u₂⋅sin(q₂) - 2⋅u₁⋅u₃⋅ ↪\n",
"↪ ────────────────────────────────────────⎟ + ──────────────────────────────── ↪\n",
"↪ 4 ⎠ 4 ↪\n",
"↪ ↪\n",
"↪ ⎞ ↪\n",
"↪ os(q₂) + 2⋅u₂⋅u₃⎠⋅q₃ ↪\n",
"↪ ──────────────────── ↪\n",
"↪ 4 ↪\n",
"↪ ↪\n",
"↪ 2⎞ ↪\n",
"↪ ⎠ ↪\n",
"↪ ── ↪\n",
"↪ ↪\n",
"\n",
"↪ ⎤\n",
"↪ cos(q₂))⋅q₃⋅cos(q₂)⎥\n",
"↪ ───────────────────⎥\n",
"↪ ⎥\n",
"↪ ⎥\n",
"↪ ⎥\n",
"↪ ⎥\n",
"↪ ⎥\n",
"↪ ⎥\n",
"↪ ⎥\n",
"↪ ⎥\n",
"↪ ⎥\n",
"↪ ⎥\n",
"↪ ⎦"
]
},
"execution_count": 20,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"gd = Frs.xreplace(ud_zerod) + Fr\n",
"gd"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\left( \\left[\\begin{matrix}q_{1}\\\\q_{2}\\\\q_{3}\\end{matrix}\\right], \\ \\left[\\begin{matrix}u_{1}\\\\u_{2}\\\\u_{3}\\end{matrix}\\right], \\ \\left[\\begin{matrix}\\dot{q}_{1}\\\\\\dot{q}_{2}\\\\\\dot{q}_{3}\\end{matrix}\\right], \\ \\left[\\begin{matrix}\\dot{u}_{1}\\\\\\dot{u}_{2}\\\\\\dot{u}_{3}\\end{matrix}\\right]\\right)$"
],
"text/plain": [
"⎛⎡q₁⎤ ⎡u₁⎤ ⎡q₁̇⎤ ⎡u₁̇⎤⎞\n",
"⎜⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥⎟\n",
"⎜⎢q₂⎥, ⎢u₂⎥, ⎢q₂̇⎥, ⎢u₂̇⎥⎟\n",
"⎜⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥⎟\n",
"⎝⎣q₃⎦ ⎣u₃⎦ ⎣q₃̇⎦ ⎣u₃̇⎦⎠"
]
},
"execution_count": 22,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"q = sm.Matrix([q1, q2, q3])\n",
"qd = q.diff(t)\n",
"q, u, qd, ud"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\left( \\left[\\begin{matrix}-1 & 0 & 0\\\\0 & -1 & 0\\\\0 & 0 & -1\\end{matrix}\\right], \\ \\left[\\begin{matrix}u_{1}\\\\u_{2}\\\\u_{3}\\end{matrix}\\right]\\right)$"
],
"text/plain": [
"⎛⎡-1 0 0 ⎤ ⎡u₁⎤⎞\n",
"⎜⎢ ⎥ ⎢ ⎥⎟\n",
"⎜⎢0 -1 0 ⎥, ⎢u₂⎥⎟\n",
"⎜⎢ ⎥ ⎢ ⎥⎟\n",
"⎝⎣0 0 -1⎦ ⎣u₃⎦⎠"
]
},
"execution_count": 24,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Mk = -sm.eye(3)\n",
"gk = u\n",
"Mk, gk"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\left( \\left[\\begin{matrix}- \\frac{l^{2} m \\cos^{2}{\\left(q_{2} \\right)}}{12} - \\frac{19 l^{2} m}{12} - \\frac{m q_{3}^{2} \\cos^{2}{\\left(q_{2} \\right)}}{4} & \\frac{l m q_{3} \\sin{\\left(q_{2} \\right)}}{4} & - \\frac{l m \\cos{\\left(q_{2} \\right)}}{4}\\\\\\frac{l m q_{3} \\sin{\\left(q_{2} \\right)}}{4} & - \\frac{l^{2} m}{12} - \\frac{m q_{3}^{2}}{4} & 0\\\\- \\frac{l m \\cos{\\left(q_{2} \\right)}}{4} & 0 & - \\frac{m}{4}\\end{matrix}\\right], \\ \\left[\\begin{matrix}- \\frac{7 g l m \\sin{\\left(q_{1} \\right)}}{4} - \\frac{g m q_{3} \\cos{\\left(q_{1} \\right)} \\cos{\\left(q_{2} \\right)}}{4} - k_{t} q_{1} + \\frac{l^{2} m u_{1} u_{2} \\sin{\\left(q_{2} \\right)} \\cos{\\left(q_{2} \\right)}}{6} - \\frac{l m q_{3} u_{1}^{2} \\cos{\\left(q_{2} \\right)}}{4} + l \\left(- \\frac{m \\left(- q_{3} u_{1}^{2} \\cos^{2}{\\left(q_{2} \\right)} - q_{3} u_{2}^{2}\\right) \\cos{\\left(q_{2} \\right)}}{4} + \\frac{m \\left(q_{3} u_{1}^{2} \\sin{\\left(q_{2} \\right)} \\cos{\\left(q_{2} \\right)} + 2 u_{2} u_{3}\\right) \\sin{\\left(q_{2} \\right)}}{4}\\right) + \\frac{m \\left(2 q_{3} u_{1} u_{2} \\sin{\\left(q_{2} \\right)} - 2 u_{1} u_{3} \\cos{\\left(q_{2} \\right)}\\right) q_{3} \\cos{\\left(q_{2} \\right)}}{4}\\\\\\frac{g m q_{3} \\sin{\\left(q_{1} \\right)} \\sin{\\left(q_{2} \\right)}}{4} - k_{t} q_{2} - \\frac{l^{2} m u_{1}^{2} \\sin{\\left(q_{2} \\right)} \\cos{\\left(q_{2} \\right)}}{12} - \\frac{m \\left(q_{3} u_{1}^{2} \\sin{\\left(q_{2} \\right)} \\cos{\\left(q_{2} \\right)} + 2 u_{2} u_{3}\\right) q_{3}}{4}\\\\- \\frac{g m \\sin{\\left(q_{1} \\right)} \\cos{\\left(q_{2} \\right)}}{4} - k_{l} q_{3} - \\frac{m \\left(- q_{3} u_{1}^{2} \\cos^{2}{\\left(q_{2} \\right)} - q_{3} u_{2}^{2}\\right)}{4}\\end{matrix}\\right]\\right)$"
],
"text/plain": [
"⎛ ↪\n",
"⎜⎡ 2 2 2 2 2 ⎤ ↪\n",
"⎜⎢ l ⋅m⋅cos (q₂) 19⋅l ⋅m m⋅q₃ ⋅cos (q₂) l⋅m⋅q₃⋅sin(q₂) -l⋅m⋅cos(q₂) ⎥ ↪\n",
"⎜⎢- ───────────── - ─────── - ────────────── ────────────── ─────────────⎥ ↪\n",
"⎜⎢ 12 12 4 4 4 ⎥ ↪\n",
"⎜⎢ ⎥ ↪\n",
"⎜⎢ 2 2 ⎥ ↪\n",
"⎜⎢ l⋅m⋅q₃⋅sin(q₂) l ⋅m m⋅q₃ ⎥, ↪\n",
"⎜⎢ ────────────── - ──── - ───── 0 ⎥ ↪\n",
"⎜⎢ 4 12 4 ⎥ ↪\n",
"⎜⎢ ⎥ ↪\n",
"⎜⎢ -l⋅m⋅cos(q₂) -m ⎥ ↪\n",
"⎜⎢ ───────────── 0 ─── ⎥ ↪\n",
"⎝⎣ 4 4 ⎦ ↪\n",
"\n",
"↪ ⎡ 2 ↪\n",
"↪ ⎢ 7⋅g⋅l⋅m⋅sin(q₁) g⋅m⋅q₃⋅cos(q₁)⋅cos(q₂) l ⋅m⋅u₁⋅u₂⋅sin(q₂)⋅cos ↪\n",
"↪ ⎢- ─────────────── - ────────────────────── - kₜ⋅q₁ + ────────────────────── ↪\n",
"↪ ⎢ 4 4 6 ↪\n",
"↪ ⎢ ↪\n",
"↪ ⎢ ↪\n",
"↪ ⎢ g ↪\n",
"↪ ⎢ ─ ↪\n",
"↪ ⎢ ↪\n",
"↪ ⎢ ↪\n",
"↪ ⎢ ↪\n",
"↪ ⎢ ↪\n",
"↪ ⎢ ↪\n",
"↪ ⎣ ↪\n",
"\n",
"↪ 2 ⎛ ⎛ 2 2 2⎞ ↪\n",
"↪ (q₂) l⋅m⋅q₃⋅u₁ ⋅cos(q₂) ⎜ m⋅⎝- q₃⋅u₁ ⋅cos (q₂) - q₃⋅u₂ ⎠⋅cos(q₂) m⋅ ↪\n",
"↪ ──── - ────────────────── + l⋅⎜- ────────────────────────────────────── + ── ↪\n",
"↪ 4 ⎝ 4 ↪\n",
"↪ ↪\n",
"↪ 2 2 ⎛ 2 ↪\n",
"↪ ⋅m⋅q₃⋅sin(q₁)⋅sin(q₂) l ⋅m⋅u₁ ⋅sin(q₂)⋅cos(q₂) m⋅⎝q₃⋅u₁ ⋅sin(q₂) ↪\n",
"↪ ───────────────────── - kₜ⋅q₂ - ──────────────────────── - ───────────────── ↪\n",
"↪ 4 12 ↪\n",
"↪ ↪\n",
"↪ ⎛ 2 2 ↪\n",
"↪ g⋅m⋅sin(q₁)⋅cos(q₂) m⋅⎝- q₃⋅u₁ ⋅cos (q₂) - q₃⋅ ↪\n",
"↪ - ─────────────────── - kₗ⋅q₃ - ────────────────────────── ↪\n",
"↪ 4 4 ↪\n",
"\n",
"↪ ⎛ 2 ⎞ ⎞ ↪\n",
"↪ ⎝q₃⋅u₁ ⋅sin(q₂)⋅cos(q₂) + 2⋅u₂⋅u₃⎠⋅sin(q₂)⎟ m⋅(2⋅q₃⋅u₁⋅u₂⋅sin(q₂) - 2⋅u₁⋅u ↪\n",
"↪ ──────────────────────────────────────────⎟ + ────────────────────────────── ↪\n",
"↪ 4 ⎠ 4 ↪\n",
"↪ ↪\n",
"↪ ⎞ ↪\n",
"↪ ⋅cos(q₂) + 2⋅u₂⋅u₃⎠⋅q₃ ↪\n",
"↪ ────────────────────── ↪\n",
"↪ 4 ↪\n",
"↪ ↪\n",
"↪ 2⎞ ↪\n",
"↪ u₂ ⎠ ↪\n",
"↪ ──── ↪\n",
"↪ ↪\n",
"\n",
"↪ ⎤⎞\n",
"↪ ₃⋅cos(q₂))⋅q₃⋅cos(q₂)⎥⎟\n",
"↪ ─────────────────────⎥⎟\n",
"↪ ⎥⎟\n",
"↪ ⎥⎟\n",
"↪ ⎥⎟\n",
"↪ ⎥⎟\n",
"↪ ⎥⎟\n",
"↪ ⎥⎟\n",
"↪ ⎥⎟\n",
"↪ ⎥⎟\n",
"↪ ⎥⎟\n",
"↪ ⎥⎟\n",
"↪ ⎦⎠"
]
},
"execution_count": 25,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Md, gd"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\left\\{g, k_{l}, k_{t}, l, m, t\\right\\}$"
],
"text/plain": [
"{g, kₗ, kₜ, l, m, t}"
]
},
"execution_count": 26,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Mk.free_symbols | gk.free_symbols | Md.free_symbols | gd.free_symbols"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle \\left[\\begin{matrix}g\\\\k_{l}\\\\k_{t}\\\\l\\\\m\\end{matrix}\\right]$"
],
"text/plain": [
"⎡g ⎤\n",
"⎢ ⎥\n",
"⎢kₗ⎥\n",
"⎢ ⎥\n",
"⎢kₜ⎥\n",
"⎢ ⎥\n",
"⎢l ⎥\n",
"⎢ ⎥\n",
"⎣m ⎦"
]
},
"execution_count": 27,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"p = sm.Matrix([g, kl, kt, l, m])\n",
"p"
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {},
"outputs": [],
"source": [
"eval_eom = sm.lambdify((q, u, p), [Mk, gk, Md, gd])"
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(array([0.43633231, 0.08726646, 0.1 ]), numpy.ndarray, (3,))"
]
},
"execution_count": 29,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"q_vals = np.array([\n",
" np.deg2rad(25.0), # q1, rad\n",
" np.deg2rad(5.0), # q2, rad\n",
" 0.1, # q3, m\n",
"])\n",
"q_vals, type(q_vals), q_vals.shape"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(array([0.1, 2.2, 0.3]), numpy.ndarray, (3,))"
]
},
"execution_count": 30,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"u_vals = np.array([\n",
" 0.1, # u1, rad/s\n",
" 2.2, # u2, rad/s\n",
" 0.3, # u3, m/s\n",
"])\n",
"u_vals, type(u_vals), u_vals.shape"
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(array([9.81, 2. , 0.01, 0.6 , 1. ]), numpy.ndarray, (5,))"
]
},
"execution_count": 31,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"p_vals = np.array([\n",
" 9.81, # g, m/s**2\n",
" 2.0, # kl, N/m\n",
" 0.01, # kt, Nm/rad\n",
" 0.6, # l, m\n",
" 1.0, # m, kg\n",
"])\n",
"p_vals, type(p_vals), p_vals.shape"
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(array([[-1, 0, 0],\n",
" [ 0, -1, 0],\n",
" [ 0, 0, -1]]),\n",
" array([[0.1],\n",
" [2.2],\n",
" [0.3]]),\n",
" array([[-0.60225313, 0.00130734, -0.1494292 ],\n",
" [ 0.00130734, -0.0325 , 0. ],\n",
" [-0.1494292 , 0. , -0.25 ]]),\n",
" array([[-4.48963535],\n",
" [-0.02486744],\n",
" [-1.1112791 ]]))"
]
},
"execution_count": 32,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"Mk_vals, gk_vals, Md_vals, gd_vals = eval_eom(q_vals, u_vals, p_vals)\n",
"Mk_vals, gk_vals, Md_vals, gd_vals"
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([0.1, 2.2, 0.3])"
]
},
"execution_count": 33,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"qd_vals = np.linalg.solve(-Mk_vals, np.squeeze(gk_vals))\n",
"qd_vals"
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([-7.46056427, -1.06525862, 0.01418834])"
]
},
"execution_count": 34,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"ud_vals = np.linalg.solve(-Md_vals, np.squeeze(gd_vals))\n",
"ud_vals"
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {},
"outputs": [],
"source": [
"def eval_rhs(t, x, p):\n",
" \"\"\"Return the right hand side of the explicit ordinary differential\n",
" equations which evaluates the time derivative of the state ``x`` at time\n",
" ``t``.\n",
"\n",
" Parameters\n",
" ==========\n",
" t : float\n",
" Time in seconds.\n",
" x : array_like, shape(6,)\n",
" State at time t: [q1, q2, q3, u1, u2, u3]\n",
" p : array_like, shape(5,)\n",
" Constant parameters: [g, kl, kt, l, m]\n",
"\n",
" Returns\n",
" =======\n",
" xd : ndarray, shape(6,)\n",
" Derivative of the state with respect to time at time ``t``.\n",
"\n",
" \"\"\"\n",
"\n",
" # unpack the q and u vectors from x\n",
" q = x[:3]\n",
" u = x[3:]\n",
"\n",
" # evaluate the equations of motion matrices with the values of q, u, p\n",
" Mk, gk, Md, gd = eval_eom(q, u, p)\n",
"\n",
" # solve for q' and u'\n",
" qd = np.linalg.solve(-Mk, np.squeeze(gk))\n",
" ud = np.linalg.solve(-Md, np.squeeze(gd))\n",
"\n",
" # pack dq/dt and du/dt into a new state time derivative vector dx/dt\n",
" xd = np.empty_like(x)\n",
" xd[:3] = qd\n",
" xd[3:] = ud\n",
"\n",
" return xd"
]
},
{
"cell_type": "code",
"execution_count": 36,
"metadata": {},
"outputs": [],
"source": [
"def euler_integrate(rhs_func, tspan, x0_vals, p_vals, delt=0.03):\n",
" \"\"\"Returns state trajectory and corresponding values of time resulting\n",
" from integrating the ordinary differential equations with Euler's\n",
" Method.\n",
"\n",
" Parameters\n",
" ==========\n",
" rhs_func : function\n",
" Python function that evaluates the derivative of the state and takes\n",
" this form ``dxdt = f(t, x, p)``.\n",
" tspan : 2-tuple of floats\n",
" The initial time and final time values: (t0, tf).\n",
" x0_vals : array_like, shape(2*n,)\n",
" Values of the state x at t0.\n",
" p_vals : array_like, shape(o,)\n",
" Values of constant parameters.\n",
" delt : float\n",
" Integration time step in seconds/step.\n",
"\n",
" Returns\n",
" =======\n",
" ts : ndarray(m, )\n",
" Monotonically equally spaced increasing values of time spaced apart\n",
" by ``delt``.\n",
" xs : ndarray(m, 2*n)\n",
" State values at each time in ts.\n",
"\n",
" \"\"\"\n",
" # generate monotonically increasing values of time.\n",
" duration = tspan[1] - tspan[0]\n",
" num_samples = round(duration/delt) + 1\n",
" ts = np.arange(tspan[0], tspan[0] + delt*num_samples, delt)\n",
"\n",
" # create an empty array to hold the state values.\n",
" x = np.empty((len(ts), len(x0_vals)))\n",
"\n",
" # set the initial conditions to the first element.\n",
" x[0, :] = x0_vals\n",
"\n",
" # use a for loop to sequentially calculate each new x.\n",
" for i, ti in enumerate(ts[:-1]):\n",
" x[i + 1, :] = x[i, :] + delt*rhs_func(ti, x[i, :], p_vals)\n",
"\n",
" return ts, x"
]
},
{
"cell_type": "code",
"execution_count": 37,
"metadata": {},
"outputs": [],
"source": [
"x0 = np.empty(6)\n",
"x0[:3] = q_vals\n",
"x0[3:] = u_vals\n",
"\n",
"t0 = 0.1"
]
},
{
"cell_type": "code",
"execution_count": 38,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([ 0.1 , 2.2 , 0.3 , -7.46056427, -1.06525862,\n",
" 0.01418834])"
]
},
"execution_count": 38,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"eval_rhs(t0, x0, p_vals)"
]
},
{
"cell_type": "code",
"execution_count": 39,
"metadata": {},
"outputs": [],
"source": [
"tf = 2.0\n",
"ts, xs = euler_integrate(eval_rhs, (t0, tf), x0, p_vals)"
]
},
{
"cell_type": "code",
"execution_count": 40,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([0.1 , 0.13, 0.16, 0.19, 0.22, 0.25, 0.28, 0.31, 0.34, 0.37, 0.4 ,\n",
" 0.43, 0.46, 0.49, 0.52, 0.55, 0.58, 0.61, 0.64, 0.67, 0.7 , 0.73,\n",
" 0.76, 0.79, 0.82, 0.85, 0.88, 0.91, 0.94, 0.97, 1. , 1.03, 1.06,\n",
" 1.09, 1.12, 1.15, 1.18, 1.21, 1.24, 1.27, 1.3 , 1.33, 1.36, 1.39,\n",
" 1.42, 1.45, 1.48, 1.51, 1.54, 1.57, 1.6 , 1.63, 1.66, 1.69, 1.72,\n",
" 1.75, 1.78, 1.81, 1.84, 1.87, 1.9 , 1.93, 1.96, 1.99])"
]
},
"execution_count": 40,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"ts"
]
},
{
"cell_type": "code",
"execution_count": 41,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([[ 4.36332313e-01, 8.72664626e-02, 1.00000000e-01,\n",
" 1.00000000e-01, 2.20000000e+00, 3.00000000e-01],\n",
" [ 4.39332313e-01, 1.53266463e-01, 1.09000000e-01,\n",
" -1.23816928e-01, 2.16804224e+00, 3.00425650e-01],\n",
" [ 4.35617805e-01, 2.18307730e-01, 1.18012770e-01,\n",
" -3.48867554e-01, 2.13303503e+00, 2.99423633e-01],\n",
" [ 4.25151779e-01, 2.82298781e-01, 1.26995479e-01,\n",
" -5.72475739e-01, 2.09463526e+00, 2.97361542e-01],\n",
" [ 4.07977506e-01, 3.45137839e-01, 1.35916325e-01,\n",
" -7.91929289e-01, 2.05197871e+00, 2.94628416e-01],\n",
" [ 3.84219628e-01, 4.06697200e-01, 1.44755177e-01,\n",
" -1.00444789e+00, 2.00382451e+00, 2.91554565e-01],\n",
" [ 3.54086191e-01, 4.66811935e-01, 1.53501814e-01,\n",
" -1.20715867e+00, 1.94877129e+00, 2.88335792e-01],\n",
" [ 3.17871431e-01, 5.25275074e-01, 1.62151888e-01,\n",
" -1.39708735e+00, 1.88551903e+00, 2.84974403e-01],\n",
" [ 2.75958810e-01, 5.81840645e-01, 1.70701120e-01,\n",
" -1.57117194e+00, 1.81314156e+00, 2.81247424e-01],\n",
" [ 2.28823652e-01, 6.36234892e-01, 1.79138543e-01,\n",
" -1.72630399e+00, 1.73132870e+00, 2.76708292e-01],\n",
" [ 1.77034532e-01, 6.88174753e-01, 1.87439792e-01,\n",
" -1.85939863e+00, 1.64055711e+00, 2.70722509e-01],\n",
" [ 1.21252574e-01, 7.37391466e-01, 1.95561467e-01,\n",
" -1.96749077e+00, 1.54215738e+00, 2.62531540e-01],\n",
" [ 6.22278505e-02, 7.83656188e-01, 2.03437413e-01,\n",
" -2.04785035e+00, 1.43826241e+00, 2.51334050e-01],\n",
" [ 7.92340069e-04, 8.26804060e-01, 2.10977435e-01,\n",
" -2.09810628e+00, 1.33164459e+00, 2.36370890e-01],\n",
" [-6.21508484e-02, 8.66753398e-01, 2.18068561e-01,\n",
" -2.11636717e+00, 1.22547077e+00, 2.17000797e-01],\n",
" [-1.25641864e-01, 9.03517521e-01, 2.24578585e-01,\n",
" -2.10132625e+00, 1.12301845e+00, 1.92757262e-01],\n",
" [-1.88681651e-01, 9.37208074e-01, 2.30361303e-01,\n",
" -2.05234005e+00, 1.02739894e+00, 1.63382255e-01],\n",
" [-2.50251853e-01, 9.68030042e-01, 2.35262771e-01,\n",
" -1.96947297e+00, 9.41324864e-01, 1.28837671e-01],\n",
" [-3.09336042e-01, 9.96269788e-01, 2.39127901e-01,\n",
" -1.85350427e+00, 8.66944026e-01, 8.92990821e-02],\n",
" [-3.64941170e-01, 1.02227811e+00, 2.41806873e-01,\n",
" -1.70589822e+00, 8.05745012e-01, 4.51378231e-02],\n",
" [-4.16118116e-01, 1.04645046e+00, 2.43161008e-01,\n",
" -1.52874283e+00, 7.58527111e-01, -3.10326579e-03],\n",
" [-4.61980401e-01, 1.06920627e+00, 2.43067910e-01,\n",
" -1.32466575e+00, 7.25419993e-01, -5.47370444e-02],\n",
" [-5.01720374e-01, 1.09096887e+00, 2.41425799e-01,\n",
" -1.09673819e+00, 7.05937072e-01, -1.08961098e-01],\n",
" [-5.34622520e-01, 1.11214698e+00, 2.38156966e-01,\n",
" -8.48377776e-01, 6.99048783e-01, -1.64889059e-01],\n",
" [-5.60073853e-01, 1.13311845e+00, 2.33210294e-01,\n",
" -5.83259685e-01, 7.03266085e-01, -2.21585204e-01],\n",
" [-5.77571643e-01, 1.15421643e+00, 2.26562738e-01,\n",
" -3.05242994e-01, 7.16728855e-01, -2.78102687e-01],\n",
" [-5.86728933e-01, 1.17571830e+00, 2.18219657e-01,\n",
" -1.83153063e-02, 7.37297421e-01, -3.33524242e-01],\n",
" [-5.87278392e-01, 1.19783722e+00, 2.08213930e-01,\n",
" 2.73444426e-01, 7.62647870e-01, -3.87002450e-01],\n",
" [-5.79075060e-01, 1.22071666e+00, 1.96603856e-01,\n",
" 5.65888142e-01, 7.90372680e-01, -4.37795193e-01],\n",
" [-5.62098415e-01, 1.24442784e+00, 1.83470001e-01,\n",
" 8.54811324e-01, 8.18087698e-01, -4.85291412e-01],\n",
" [-5.36454076e-01, 1.26897047e+00, 1.68911258e-01,\n",
" 1.13596271e+00, 8.43544507e-01, -5.29023016e-01],\n",
" [-5.02375194e-01, 1.29427680e+00, 1.53040568e-01,\n",
" 1.40505807e+00, 8.64744480e-01, -5.68660838e-01],\n",
" [-4.60223452e-01, 1.32021914e+00, 1.35980743e-01,\n",
" 1.65780471e+00, 8.80047808e-01, -6.03995480e-01],\n",
" [-4.10489311e-01, 1.34662057e+00, 1.17860878e-01,\n",
" 1.88994245e+00, 8.88268730e-01, -6.34907053e-01],\n",
" [-3.53791037e-01, 1.37326863e+00, 9.88136667e-02,\n",
" 2.09730514e+00, 8.88747420e-01, -6.61330213e-01],\n",
" [-2.90871883e-01, 1.39993105e+00, 7.89737603e-02,\n",
" 2.27590461e+00, 8.81389870e-01, -6.83221797e-01],\n",
" [-2.22594745e-01, 1.42637275e+00, 5.84771064e-02,\n",
" 2.42203585e+00, 8.66668771e-01, -7.00537408e-01],\n",
" [-1.49933669e-01, 1.45237281e+00, 3.74609842e-02,\n",
" 2.53239896e+00, 8.45580400e-01, -7.13220652e-01],\n",
" [-7.39617005e-02, 1.47774023e+00, 1.60643646e-02,\n",
" 2.60422940e+00, 8.19554682e-01, -7.21205288e-01],\n",
" [ 4.16518154e-03, 1.50232687e+00, -5.57179403e-03,\n",
" 2.63542555e+00, 7.90319326e-01, -7.24427319e-01],\n",
" [ 8.32279481e-02, 1.52603645e+00, -2.73046136e-02,\n",
" 2.62466048e+00, 7.59725678e-01, -7.22841890e-01],\n",
" [ 1.61967762e-01, 1.54882822e+00, -4.89898703e-02,\n",
" 2.57146487e+00, 7.29553925e-01, -7.16439450e-01],\n",
" [ 2.39111708e-01, 1.57071483e+00, -7.04830538e-02,\n",
" 2.47627011e+00, 7.01325433e-01, -7.05256667e-01],\n",
" [ 3.13399812e-01, 1.59175460e+00, -9.16407538e-02,\n",
" 2.34040397e+00, 6.76154998e-01, -6.89379799e-01],\n",
" [ 3.83611931e-01, 1.61203925e+00, -1.12322148e-01,\n",
" 2.16603688e+00, 6.54671168e-01, -6.68940608e-01],\n",
" [ 4.48593037e-01, 1.63167938e+00, -1.32390366e-01,\n",
" 1.95608296e+00, 6.37018712e-01, -6.44106831e-01],\n",
" [ 5.07275526e-01, 1.65078994e+00, -1.51713571e-01,\n",
" 1.71406620e+00, 6.22938813e-01, -6.15070314e-01],\n",
" [ 5.58697512e-01, 1.66947811e+00, -1.70165680e-01,\n",
" 1.44396681e+00, 6.11906749e-01, -5.82036029e-01],\n",
" [ 6.02016516e-01, 1.68783531e+00, -1.87626761e-01,\n",
" 1.15006549e+00, 6.03298605e-01, -5.45214730e-01],\n",
" [ 6.36518481e-01, 1.70593427e+00, -2.03983203e-01,\n",
" 8.36802201e-01, 5.96558484e-01, -5.04821104e-01],\n",
" [ 6.61622547e-01, 1.72383102e+00, -2.19127836e-01,\n",
" 5.08663225e-01, 5.91343485e-01, -4.61078217e-01],\n",
" [ 6.76882444e-01, 1.74157133e+00, -2.32960183e-01,\n",
" 1.70104973e-01, 5.87631757e-01, -4.14227908e-01],\n",
" [ 6.81985593e-01, 1.75920028e+00, -2.45387020e-01,\n",
" -1.74482332e-01, 5.85786780e-01, -3.64545604e-01],\n",
" [ 6.76751123e-01, 1.77677388e+00, -2.56323388e-01,\n",
" -5.20773373e-01, 5.86577261e-01, -3.12357072e-01],\n",
" [ 6.61127922e-01, 1.79437120e+00, -2.65694100e-01,\n",
" -8.64475606e-01, 5.91156371e-01, -2.58054046e-01],\n",
" [ 6.35193654e-01, 1.81210589e+00, -2.73435722e-01,\n",
" -1.20127673e+00, 6.01006553e-01, -2.02105739e-01],\n",
" [ 5.99155352e-01, 1.83013609e+00, -2.79498894e-01,\n",
" -1.52678495e+00, 6.17857094e-01, -1.45063860e-01],\n",
" [ 5.53351803e-01, 1.84867180e+00, -2.83850810e-01,\n",
" -1.83647709e+00, 6.43581213e-01, -8.75598764e-02],\n",
" [ 4.98257491e-01, 1.86797924e+00, -2.86477606e-01,\n",
" -2.12567032e+00, 6.80078039e-01, -3.02944186e-02],\n",
" [ 4.34487381e-01, 1.88838158e+00, -2.87386438e-01,\n",
" -2.38953235e+00, 7.29142906e-01, 2.59803906e-02],\n",
" [ 3.62801410e-01, 1.91025587e+00, -2.86607027e-01,\n",
" -2.62314224e+00, 7.92327527e-01, 8.04846656e-02],\n",
" [ 2.84107143e-01, 1.93402569e+00, -2.84192487e-01,\n",
" -2.82160885e+00, 8.70790639e-01, 1.32439931e-01],\n",
" [ 1.99458878e-01, 1.96014941e+00, -2.80219289e-01,\n",
" -2.98024641e+00, 9.65140541e-01, 1.81104004e-01],\n",
" [ 1.10051485e-01, 1.98910363e+00, -2.74786169e-01,\n",
" -3.09479777e+00, 1.07527459e+00, 2.25812638e-01]])"
]
},
"execution_count": 41,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"xs"
]
},
{
"cell_type": "code",
"execution_count": 42,
"metadata": {},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "code",
"execution_count": 43,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"[<matplotlib.lines.Line2D at 0x1e335aa5870>,\n",
" <matplotlib.lines.Line2D at 0x1e335aa5960>,\n",
" <matplotlib.lines.Line2D at 0x1e335aa5a50>,\n",
" <matplotlib.lines.Line2D at 0x1e335aa5b40>,\n",
" <matplotlib.lines.Line2D at 0x1e335aa5c30>,\n",
" <matplotlib.lines.Line2D at 0x1e335aa5d20>]"
]
},
"execution_count": 43,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXIAAAD4CAYAAADxeG0DAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAAsTAAALEwEAmpwYAABaA0lEQVR4nO3dd3hb1fnA8e/VsGzZ8t57ZC+yCCQkBEgCYY+wwoaWWVZpKVAo8CulQKEUKKVllx02LSEQMlgNZJG9vfdeki1rn98f13EcyLalK9nn8zz3kSxd6b6+vn51dKYihECSJEkKXTqtA5AkSZL6RiZySZKkECcTuSRJUoiTiVySJCnEyUQuSZIU4gxaHDQxMVHk5uZqcWhJkqSQ9eOPPzYJIZJ++rgmiTw3N5e1a9dqcWhJkqSQpShK+b4el1UrkiRJIU4mckmSpBAnE7kkSVKIk4lckiQpxMlELkmSFOJkIpckSQpxMpFLkiSFOE36kUuSdGDetja6Nm7EsWMniikMQ1wc+rg49LGxGFJSMKakaB2iFERkIpekICDcbqyff07nDyvp2rABV2npAfePOOooYs47j+jTTkVvsQQoSilYKVosLDF58mQhR3ZKEgiPh/b/fkrTc8/hrqpCHxtLxPjxPVv4mDHg9eBta8Pb2oqntRVXcTFtn3yCq6gYJTwcy8lziLvoIsyTJmn960h+pijKj0KIyT97XCZySQo84fNhXfQ5Tc8+i6usjPBRo0i89RaiZs5EUZSDv14IHJs30/bRR1g/W4TPZiP2gvNJvusu9FFRAfgNJC3IRC5JQcJdXU3V7b/GsXkzpqFDSbz1FiyzZx9SAt8Xn8NB0z+eo/nllzGmppL25z8Teewx/Ry1FAz2l8hlrxVJCqCOFSsonXc+rtJS0v/yGHn/+YToOXOOOIkD6MLDSf7NHeS89SaK0UjFVVdR9/Cf8XV19WPkUjCTiVySAkAIQdMLL1J57XUYkhLJ++B9Ys46C0XXf/+C5gkTyPv4I+Iuu4zWN96g7KKL8TQ19dv7S8FLJnJJ8jNvRyfVt95G45NPEj33FHIXLCDMT/Px68xmUu+7l6wXX8RVWUn55Vfgrq/3y7Gk4CETuST5kddmo+LKK7EtX07yXXeR/te/oouM9Ptxo2ZMJ/ulF/E0NKjJvKbG78eUtCMTuST5ia+zk8rrrsexaxeZz/6dhKuv6lNd+OEyT5pE9isv421tpfyyy3FVVgbs2FJgyUQuSX7gczio/NXNdG3cSMYTT2A58URN4og46iiyX30VX2cn5ZddjvMgA42k0CQTuST1M+FyUXXbbdhXrSL90UeIPuVkTeOJGDOa7NdfQ7jdVF53Pd62Nk3jkfqfTOSS1I+Ex0P1b++k85tvSX3wQWLOOkvrkAAIHz6crOf+gaeujuo7foPweLQOSepHMpFLUj+q//OfsX35JSn33E3cRRdqHc5eIsaPJ+X+P9D5/fc0/O1vWocj9SM5aZYk9ZO2Tz6h9e13iL/6auKvvFLrcPYp7oILcG7fTsvLrxA+ahQxp5+udUhSPwipEnn51mY2Lq+kodyK1+vTOhxJ6uHYto26Bx7EPGUKyb+5Q+twDijl7ruJmDSJ2nvvw7F9u9bhSP0gpOZa+frtnWz9thoAQ5iOlNxoUvNjSBsaS1pBDGHh8guGFHjetjZKz78A4XaT99GHGBIStA7poDxNTZTOOx9Fryf3ww8wxMVpHZJ0CAbMpFkdrQ5qi9upK26nrqSdxsoOhE+g6BSSsi1kDIslY1gcaUNkYpf8T3i9VN5wI50rV5L7xutEjB+vdUiHrGvzZsovvYzIGTPIfPbvAe3jLh2Z/SXykMt0UXHhDJ0cztDJ6gopbqeXupJ2qne1UrOrjY3LKln/ZQU6vUJqfgzZo+PJGhlPUpYFRScvVKl/Nf3jH3R+9x2pDz4QUkkcIGLsWJJuu5WGx5/AuvAzYs48Q+uQpCMUciXyg3E7vdQVt1O5o4XK7S00VXYAEB5lJGd0ArnjEskeHS9L61KfdXz3PyqvvZaYc88l7c8Ph2SJVni9lF9yKa6yMvIXfoohKUnrkKQD8FvViqIoWcDrQAoggBeEEE8f6DWBnI+8s91J1Y5WKrY2U761GWenB51BIXNYHLnjEimYmIw5OiwgsUgDh9dqpeTMs9BFRZH34QfowsO1DumIOUtKKD3nXCKPn0Hm32UVSzDzZyJPA9KEEOsURbEAPwLnCCG27e81Wi0s4fP6qCtpp3RjE6Ubm2hv7EJRIGN4HEMnp5A/IYnwSGPA45JCT83d99D+6afkLniHiLFjtQ6nz5pffoWGxx8n/fHHZRVLEAtYY6eiKP8BnhVCLNnfPsGwQpAQgpbaTorWNlC4pp72xi50eoXsUfGMPC6dnLEJ6PUh1TtTChDb8uVU3fQrEm64nuTbb9c6nH4hq1gCwN0FtZsgZRSYjmzB7IAkckVRcoFvgTFCCOtPnrsOuA4gOzt7Unl5eb8dt6+EEDRVdlC4pp5dq+vobHdhjg5jxNQ0Rh6XRmyyWesQpSDhaW2l5MyzMCQkkPf+eyhhA6darqeKRfZi6TuPExq2Q816qFkH1euhYRsIL1zyPgw7svl3/J7IFUWJAr4BHhZCfHSgfYOhRL4/Pq+P8q0tbPtfDeVbmhE+QeaIOMbPziZ7dLy8uAe56jt+g/XLL8l7/z3CR47UOpx+t7uKJeOpvxE9d67W4QQ/IaCjARq3Q90WqNusbk07wdc9n014LGRMhPSJkD4BcqaBOf6IDufXRK4oihFYCCwWQjx5sP2DOZH31tnmZPv3tWz5tprONidxaZGMn53F8Cmp6I2y2mWwsX6xmOrbbyfx1ltIuukmrcPxC+H1UjrvfLzWdgoWLQrpRtx+JQTYaqFpFzQVQuMOtcTdsB26WvbsF5UKqWO7tzFq4o7Lg34qAPqzsVMBXgNahBC3H8prQiWR7+b1+Cj6sYH1SyporuogIjqM8bOyGHtCJkaTXuvwpADwtLZSctrpGNPTyV3wDopx4DaKd65eTcUVVw7oD6z96mqDlmJoLum+LVK3pkJwdezZzxQNSSMgeaS6JY2AlNEQlezX8PyZyKcD3wGbgd0ToPxeCLFof68JtUS+mxCCqp2trP+ygsptLURYjEyam8vo49MxGGVCH8hq7ruP9k/+Q95HHxI+bJjW4fhd1e2/puPrryn4fBHGtDStw+k/Ph901EFrObSWQktp922Jer936RoFYrMgvgASh0Hi0O7bYWBJ7bdS9uEYMEP0g0VtcTur/ltC9c5WouJMTD4tlxHT0mRPlwGoa8MGyi6eT/w115Dyuzu1DicgXFXVlJx+OpY5c8h44nGtwzk8XW3QVq4m671uy9T7XueefRUdxGSq1R/x+RCfpybuhCEQlwvG4KpakoncT6p2tLDqvyXUlViJSzUz/YKhZI8O/kmTpEMjvF5KL7gAb1Mz+YsWoY/y/8LJwaLxmWdoeu6f5Lz9FuaJE7UOZw+HFdoq1OTcVqFurb3uO9v33j88BmJz1MQct/s2F2JzITYbDKHT80gmcj8SQlC2uZkVHxTS3tBF7rhEjjt/iOy2OAC0vPUW9Q/9iYy/PUn0qadqHU5A+ex2ik89DUNiIrnvv4eiC9C3TadtT1Lu2cr3JGxH2977G81qoo7NVhN1TNaehB2bAxGxgYk7AGQiDwCv28fG5ZWsXVSG1+tj/KxsJp+WKxtEQ5SnuZniuacSMXYMWS+/PCi7nrZ/upCaO+8k7eE/ETtvXv+8qavz50m6d4l6r3pqwBChJundW1xOr59z1a58g+RvIxN5AHW2O/nh42J2rqzDkhDOiZeNIGvkkfUblbRTc8/vaV+4kPz/fIIpP1/rcDQhhFBHfFZUULD4C/RRUQd/kccF7ZVqnfReSbr7vr1p7/31pn0k6Jw9t5GJgyZRH4xM5BqoKWzjqzd30FZvZ+S0NKbNGyLncgkR9nXrKL/kUhKuvTboV/zxt65Nmyi78KI93RGFgI76PY2HPQm7TN2sNajz53XTGdXeH70T9O5qj9hsiEyCQFXbhDiZyDXicXtZ81kZ67+sICLKyPH
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.plot(ts, xs)"
]
},
{
"cell_type": "code",
"execution_count": 44,
"metadata": {},
"outputs": [],
"source": [
"from scipy.integrate import solve_ivp"
]
},
{
"cell_type": "code",
"execution_count": 45,
"metadata": {},
"outputs": [],
"source": [
"result = solve_ivp(eval_rhs, (t0, tf), x0, args=(p_vals,))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [
{
"ename": "ValueError",
"evalue": "x and y must have same first dimension, but have shapes (10,) and (6, 10)",
"output_type": "error",
"traceback": [
"\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[1;31mValueError\u001b[0m Traceback (most recent call last)",
"\u001b[1;32m~\\AppData\\Local\\Temp\\ipykernel_38236\\1859627624.py\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0mplt\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mplot\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mresult\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mt\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mresult\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0my\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
"\u001b[1;32mc:\\Users\\13063\\.conda\\envs\\MinerU\\lib\\site-packages\\matplotlib\\pyplot.py\u001b[0m in \u001b[0;36mplot\u001b[1;34m(scalex, scaley, data, *args, **kwargs)\u001b[0m\n\u001b[0;32m 3706\u001b[0m \u001b[1;33m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 3707\u001b[0m ) -> list[Line2D]:\n\u001b[1;32m-> 3708\u001b[1;33m return gca().plot(\n\u001b[0m\u001b[0;32m 3709\u001b[0m \u001b[1;33m*\u001b[0m\u001b[0margs\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 3710\u001b[0m \u001b[0mscalex\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mscalex\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
"\u001b[1;32mc:\\Users\\13063\\.conda\\envs\\MinerU\\lib\\site-packages\\matplotlib\\axes\\_axes.py\u001b[0m in \u001b[0;36mplot\u001b[1;34m(self, scalex, scaley, data, *args, **kwargs)\u001b[0m\n\u001b[0;32m 1777\u001b[0m \"\"\"\n\u001b[0;32m 1778\u001b[0m \u001b[0mkwargs\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mcbook\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mnormalize_kwargs\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mmlines\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mLine2D\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m-> 1779\u001b[1;33m \u001b[0mlines\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m[\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_get_lines\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m*\u001b[0m\u001b[0margs\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mdata\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mdata\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 1780\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[0mline\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mlines\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 1781\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0madd_line\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mline\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
"\u001b[1;32mc:\\Users\\13063\\.conda\\envs\\MinerU\\lib\\site-packages\\matplotlib\\axes\\_base.py\u001b[0m in \u001b[0;36m__call__\u001b[1;34m(self, axes, data, *args, **kwargs)\u001b[0m\n\u001b[0;32m 294\u001b[0m \u001b[0mthis\u001b[0m \u001b[1;33m+=\u001b[0m \u001b[0margs\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 295\u001b[0m \u001b[0margs\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0margs\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 296\u001b[1;33m yield from self._plot_args(\n\u001b[0m\u001b[0;32m 297\u001b[0m axes, this, kwargs, ambiguous_fmt_datakey=ambiguous_fmt_datakey)\n\u001b[0;32m 298\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
"\u001b[1;32mc:\\Users\\13063\\.conda\\envs\\MinerU\\lib\\site-packages\\matplotlib\\axes\\_base.py\u001b[0m in \u001b[0;36m_plot_args\u001b[1;34m(self, axes, tup, kwargs, return_kwargs, ambiguous_fmt_datakey)\u001b[0m\n\u001b[0;32m 484\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 485\u001b[0m \u001b[1;32mif\u001b[0m \u001b[0mx\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m!=\u001b[0m \u001b[0my\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 486\u001b[1;33m raise ValueError(f\"x and y must have same first dimension, but \"\n\u001b[0m\u001b[0;32m 487\u001b[0m f\"have shapes {x.shape} and {y.shape}\")\n\u001b[0;32m 488\u001b[0m \u001b[1;32mif\u001b[0m \u001b[0mx\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mndim\u001b[0m \u001b[1;33m>\u001b[0m \u001b[1;36m2\u001b[0m \u001b[1;32mor\u001b[0m \u001b[0my\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mndim\u001b[0m \u001b[1;33m>\u001b[0m \u001b[1;36m2\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
"\u001b[1;31mValueError\u001b[0m: x and y must have same first dimension, but have shapes (10,) and (6, 10)"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy80BEi2AAAACXBIWXMAAAsTAAALEwEAmpwYAAANT0lEQVR4nO3cYYjkd33H8ffHO1NpjKb0VpC706T00njYQtIlTRFqirZc8uDugUXuIFgleGAbKVWEFEuU+MiGWhCu1ZOKVdAYfSALntwDjQTEC7chNXgXItvTeheFrDHNk6Ax7bcPZtKdrneZf3Zndy/7fb/gYP7/+e3Mlx97752d2ZlUFZKk7e8VWz2AJGlzGHxJasLgS1ITBl+SmjD4ktSEwZekJqYGP8lnkzyZ5PuXuD5JPplkKcmjSW6c/ZiSpPUa8gj/c8CBF7n+VmDf+N9R4F/WP5YkadamBr+qHgR+/iJLDgGfr5FTwNVJXj+rASVJs7FzBrexGzg/cXxhfO6nqxcmOcrotwCuvPLKP7z++utncPeS1MfDDz/8s6qaW8vXziL4g1XVceA4wPz8fC0uLm7m3UvSy16S/1zr187ir3SeAPZOHO8Zn5MkXUZmEfwF4F3jv9a5GXimqn7t6RxJ0taa+pROki8BtwC7klwAPgK8EqCqPgWcAG4DloBngfds1LCSpLWbGvyqOjLl+gL+emYTSZI2hO+0laQmDL4kNWHwJakJgy9JTRh8SWrC4EtSEwZfkpow+JLUhMGXpCYMviQ1YfAlqQmDL0lNGHxJasLgS1ITBl+SmjD4ktSEwZekJgy+JDVh8CWpCYMvSU0YfElqwuBLUhMGX5KaMPiS1ITBl6QmDL4kNWHwJakJgy9JTRh8SWrC4EtSEwZfkpow+JLUhMGXpCYMviQ1YfAlqYlBwU9yIMnjSZaS3HWR69+Q5IEkjyR5NMltsx9VkrQeU4OfZAdwDLgV2A8cSbJ/1bK/B+6vqhuAw8A/z3pQSdL6DHmEfxOwVFXnquo54D7g0Ko1BbxmfPm1wE9mN6IkaRaGBH83cH7i+ML43KSPArcnuQCcAN5/sRtKcjTJYpLF5eXlNYwrSVqrWb1oewT4XFXtAW4DvpDk1267qo5X1XxVzc/Nzc3oriVJQwwJ/hPA3onjPeNzk+4A7geoqu8CrwJ2zWJASdJsDAn+aWBfkmuTXMHoRdmFVWt+DLwNIMmbGAXf52wk6TIyNfhV9TxwJ3ASeIzRX+OcSXJPkoPjZR8E3pvke8CXgHdXVW3U0JKkl27nkEVVdYLRi7GT5+6euHwWeMtsR5MkzZLvtJWkJgy+JDVh8CWpCYMvSU0YfElqwuBLUhMGX5KaMPiS1ITBl6QmDL4kNWHwJakJgy9JTRh8SWrC4EtSEwZfkpow+JLUhMGXpCYMviQ1YfAlqQmDL0lNGHxJasLgS1ITBl+SmjD4ktSEwZekJgy+JDVh8CWpCYMvSU0YfElqwuBLUhMGX5KaMPiS1ITBl6QmDL4kNTEo+EkOJHk8yVKSuy6x5p1JziY5k+SLsx1TkrReO6ctSLIDOAb8GXABOJ1koarOTqzZB/wd8JaqejrJ6zZqYEnS2gx5hH8TsFRV56rqOeA+4NCqNe8FjlXV0wBV9eRsx5QkrdeQ4O8Gzk8cXxifm3QdcF2S7yQ5leTAxW4oydEki0kWl5eX1zaxJGlNZvWi7U5gH3ALcAT4TJKrVy+qquNVNV9V83NzczO6a0nSEEOC/wSwd+J4z/jcpAvAQlX9qqp+CPyA0Q8ASdJlYkjwTwP7klyb5ArgMLCwas3XGD26J8kuRk/xnJvdmJKk9Zoa/Kp6HrgTOAk8BtxfVWeS3JPk4HjZSeCpJGeBB4APVdVTGzW0JOmlS1VtyR3Pz8/X4uLilty3JL1cJXm4qubX8rW+01aSmjD4ktSEwZekJgy+JDVh8CWpCYMvSU0YfElqwuBLUhMGX5KaMPiS1ITBl6QmDL4kNWHwJakJgy9JTRh8SWrC4EtSEwZfkpow+JLUhMGXpCYMviQ1YfAlqQmDL0lNGHxJasLgS1ITBl+SmjD4ktSEwZekJgy+JDVh8CWpCYMvSU0YfElqwuBLUhMGX5KaMPiS1ITBl6QmBgU/yYEkjydZSnLXi6x7R5JKMj+7ESVJszA1+El2AMeAW4H9wJEk+y+y7irgb4CHZj2kJGn9hjzCvwlYqqpzVfUccB9w6CLrPgZ8HPjFDOeTJM3IkODvBs5PHF8Yn/s/SW4E9lbV11/shpIcTbKYZHF5efklDytJWrt1v2ib5BXAJ4APTltbVcerar6q5ufm5tZ715Kkl2BI8J8A9k4c7xmfe8FVwJuBbyf5EXAzsOALt5J0eRkS/NPAviTXJrkCOAwsvHBlVT1TVbuq6pqqugY4BRysqsUNmViStCZTg19VzwN3AieBx4D7q+pMknuSHNzoASVJs7FzyKKqOgGcWHXu7kusvWX9Y0mSZs132kpSEwZfkpow+JLUhMGXpCYMviQ1YfAlqQmDL0lNGHxJasLgS1ITBl+SmjD4ktSEwZekJgy+JDVh8CWpCYMvSU0YfElqwuBLUhMGX5KaMPiS1ITBl6QmDL4kNWHwJakJgy9JTRh8SWrC4EtSEwZfkpow+JLUhMGXpCYMviQ1YfAlqQmDL0lNGHxJasLgS1ITBl+SmhgU/CQHkjyeZCnJXRe5/gNJziZ5NMk3k7xx9qNKktZjavCT7ACOAbcC+4EjSfavWvYIMF9VfwB8FfiHWQ8qSVqfIY/wbwKWqupcVT0H3AccmlxQVQ9U1bPjw1PAntmOKUlaryHB3w2cnzi+MD53KXcA37jYFUmOJllMsri8vDx8SknSus30RdsktwPzwL0Xu76qjlfVfFXNz83NzfKuJUlT7Byw5glg78TxnvG5/yfJ24EPA2+tql/OZjxJ0qwMeYR/GtiX5NokVwCHgYXJBUluAD4NHKyqJ2c/piRpvaYGv6qeB+4ETgKPAfdX1Zkk9yQ5OF52L/Bq4CtJ/j3JwiVuTpK0RYY8pUNVnQBOrDp398Tlt894LknSjPlOW0lqwuBLUhMGX5KaMPiS1ITBl6QmDL4kNWHwJakJgy9JTRh8SWrC4EtSEwZfkpow+JLUhMGXpCYMviQ1YfAlqQmDL0lNGHxJasLgS1ITBl+SmjD4ktSEwZekJgy+JDVh8CWpCYMvSU0YfElqwuBLUhMGX5KaMPiS1ITBl6QmDL4kNWHwJakJgy9JTRh8SWrC4EtSEwZfkpoYFPwkB5I8nmQpyV0Xuf43knx5fP1DSa6Z+aSSpHWZGvwkO4BjwK3AfuBIkv2rlt0BPF1Vvwv8E/DxWQ8qSVqfIY/wbwKWqupcVT0H3AccWrXmEPBv48tfBd6WJLMbU5K0XjsHrNkNnJ84vgD80aXWVNXzSZ4Bfhv42eSiJEeBo+PDXyb5/lqG3oZ2sWqvGnMvVrgXK9yLFb+31i8cEvyZqarjwHGAJItVNb+Z93+5ci9WuBcr3IsV7sWKJItr/dohT+k8AeydON4zPnfRNUl2Aq8FnlrrUJKk2RsS/NPAviTXJrkCOAwsrFqzAPzl+PJfAN+qqprdmJKk9Zr6lM74Ofk7gZPADuCzVXUmyT3AYlUtAP8KfCHJEvBzRj8Upjm+jrm3G/dihXuxwr1Y4V6sWPNexAfiktSD77SVpCYMviQ1seHB92MZVgzYiw8kOZvk0STfTPLGrZhzM0zbi4l170hSSbbtn+QN2Ysk7xx/b5xJ8sXNnnGzDPg/8oYkDyR5ZPz/5LatmHOjJflskicv9V6ljHxyvE+PJrlx0A1X1Yb9Y/Qi738AvwNcAXwP2L9qzV8BnxpfPgx8eSNn2qp/A/fiT4HfHF9+X+e9GK+7CngQOAXMb/XcW/h9sQ94BPit8fHrtnruLdyL48D7xpf3Az/a6rk3aC/+BLgR+P4lrr8N+AYQ4GbgoSG3u9GP8P1YhhVT96KqHqiqZ8eHpxi952E7GvJ9AfAxRp/L9IvNHG6TDdmL9wLHquppgKp6cpNn3CxD9qKA14wvvxb4ySbOt2mq6kFGf/F4KYeAz9fIKeD
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.plot(result.t, np.transpose(result.y))"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "MinerU",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.10.16"
}
},
"nbformat": 4,
"nbformat_minor": 2
}