{
"cells": [
{
"cell_type": "markdown",
"id": "b8b9cbb6-7d46-4378-9df9-9a6287861d89",
"metadata": {},
"source": [
"# MT94/P26/TD - Equations différentielles, systèmes dynamiques"
]
},
{
"cell_type": "markdown",
"id": "0fb02a9b-e519-49de-8664-987ab9e2f2e4",
"metadata": {},
"source": [
"## Entrée en matière\n",
"### Schémas à un pas"
]
},
{
"cell_type": "markdown",
"id": "edd2ba70-a179-4a0f-b6a4-ea1ea6dff465",
"metadata": {},
"source": [
"On considère la fonction $y:[0,T]\\rightarrow\\mathbb{R}$ solution d'un problème de Cauchy, c'est à dire une équation différentielle munie d'une condition initiale en $t=0$\n",
"\\begin{align}\n",
"y'&=f(t,y),\\,t\\in[0,T]\\\\\n",
"y(0)&=y_0.\n",
"\\end{align}\n",
"On suppose que $f$ vérifie les conditions du théorème de Cauchy-Lipschitz, ce qui garantit l'existence et l'unicité de la solution. Pour $N\\in\\mathbb{N}^*$ on considère une discrétisation uniforme de l'intervalle $[0,T]$ sous la forme d'une suite $(t_n)_{n=0\\dots N}$ définie par $t_n=nh$ où $h=T/N$. \n",
"\n",
"Un schéma à un pas a pour but de construire une suite $(y_n)_{n=0\\dots N}$ approchant les valeurs de la solution $y$ aux points $t_n$, c'est à dire $(y(t_n))_{n=0\\dots N}$, en considérant l'équation de récurrence\n",
"$$\n",
"y_{n+1}=y_n+h\\phi(t_n,y_n;h),\\,n\\in\\{0\\dots N-1\\}.\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "ac3a58c2-8178-4b39-96cd-4a07d83fe1db",
"metadata": {},
"source": [
"En cours nous avons vu les schémas suivants :\n",
"1. Schéma d'Euler\n",
"$$\n",
"y_{n+1}=y_n+hf(t_n,y_n).\n",
"$$\n",
"2. Schéma d'Euler-Cauchy\n",
"\\begin{align*}\n",
"k_1&=f(t_n,y_n),\\\\\n",
"k_2&=f(t_n+h,y_n+hk_1),\\\\\n",
"y_{n+1}&=y_n+\\frac{h}{2}(k_1+k_2).\n",
"\\end{align*}\n",
"\n",
"3. Schéma du point milieu\n",
"\\begin{align*}\n",
"k_1&=f(t_n,y_n),\\\\\n",
"k_2&=f(t_n+\\frac{h}{2},y_n+\\frac{h}{2} k_1),\\\\\n",
"y_{n+1}&=y_n+hk_2.\n",
"\\end{align*}\n",
"4. Schéma de Runge et Kutta d'ordre 4\n",
"\\begin{align*}\n",
"k_1&=f(t_n,y_n),\\\\\n",
"k_2&=f(t_n+\\frac{h}{2},y_n+\\frac{h}{2} k_1),\\\\\n",
"k_3&=f(t_n+\\frac{h}{2},y_n+\\frac{h}{2} k_2),\\\\\n",
"k_4&=f(t_n+h,y_n+hk_3),\\\\\n",
"y_{n+1}&=y_n+\\frac{h}{6}(k_1+2k_2+2k_3+k_4).\n",
"\\end{align*}\n"
]
},
{
"cell_type": "markdown",
"id": "cae4309c-bff4-424c-abd7-a4f910ed9aa5",
"metadata": {},
"source": [
"## Evaluation de l'ordre des schémas\n",
"Nous allons évaluer expérimentalement l'ordre de ces schémas, c'est à dire la valeur de l'entier $p$ tel que\n",
"\\begin{equation}\n",
"e(h)=Ch^p+h^p\\epsilon(h),\n",
"\\end{equation}\n",
"où\n",
"$$\n",
"e(h)=\\max_{n\\in\\{0\\dots N\\}} |y(t_n)-y_n|,\n",
"$$\n",
"et où $t\\rightarrow y(t)$ est la solution de l'équation différentielle :\n",
"\\begin{align}\n",
"y'&=-ty+t,\\,t\\in[0,4],\\tag{2}\\\\\n",
"y(0)&=0.\n",
"\\end{align}\n",
"Voici comment organiser votre travail :"
]
},
{
"cell_type": "markdown",
"id": "2fb0430d-683b-4deb-98f0-2ad3284c6dc3",
"metadata": {},
"source": [
"1. Calculez la solution exacte de (2)."
]
},
{
"cell_type": "code",
"execution_count": 62,
"id": "e2284bf9-1677-4435-a0d6-57e649e726e3",
"metadata": {},
"outputs": [],
"source": [
"ytrue=#(t)->(1-exp(-0.5*t.^2));"
]
},
{
"cell_type": "markdown",
"id": "55a7ee82-35f5-41fd-a516-06f400048abc",
"metadata": {},
"source": [
"2. Ecrivez une fonction scilab calculant le second-membre de l'équation différentielle (E) :"
]
},
{
"cell_type": "code",
"execution_count": 63,
"id": "ad419712-5526-4047-97c2-6ec31a6021e7",
"metadata": {},
"outputs": [],
"source": [
"function dydt=f1(t,y)\n",
" dydt=-t*y+t;\n",
"endfunction"
]
},
{
"cell_type": "markdown",
"id": "409b605e-c60b-4ce0-9d71-62ab7684e4a0",
"metadata": {},
"source": [
"3. Ecrivez une fonction Scilab acceptant comme argument $y_0$, un vecteur contenant $(t_n)_{n=0\\dots N}$ et l'identificateur de la fonction Scilab calculant $f$, et renvoyant dans un vecteur les approximations $y_n$ pour $n=0\\dots N$ données par le Schéma d'Euler :"
]
},
{
"cell_type": "code",
"execution_count": 64,
"id": "0d77bc15-f4ef-4632-9282-dc7087e7daa4",
"metadata": {},
"outputs": [],
"source": [
"function y=euler(f,t,y0)\n",
" n=length(t);\n",
" y=zeros(1,n);\n",
" h=t(2)-t(1);\n",
" y(1)=y0;\n",
" for k=1:n-1\n",
" y(k+1)=y(k)+h*f(t(k),y(k));\n",
" end\n",
"endfunction"
]
},
{
"cell_type": "markdown",
"id": "4ab6db49-8d13-4d31-9bd3-bc726128b884",
"metadata": {},
"source": [
"4. Testez votre fonction en représentant graphiquement les approximations $y_n$ superposées à la solution exacte $y(t)$, ceci pour des valeurs croissantes de $N$."
]
},
{
"cell_type": "code",
"execution_count": 65,
"id": "662c2ce3-7001-4caa-9fb6-6185ca8f6f47",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
""
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"for N = [10]\n",
" t = linspace(0,4,N+1);\n",
" u = linspace(0,4,500);\n",
" y = euler(f1,t,0);\n",
" scf(); // cree une nouvelle fenetre graphique\n",
" plot(t,y,\"x\",u,ytrue(u))\n",
" legend(\"Euler\",\"$1-\\exp(-t^2/2)$\",4)\n",
" title(\"N=\"+string(N))\n",
"end\n"
]
},
{
"cell_type": "markdown",
"id": "fa1ac3ad-1222-4531-ab65-36188dd9d8c6",
"metadata": {},
"source": [
"5. Calculez et représentez graphiquement l'erreur de convergence $e(h)$ pour\n",
"$$\n",
"N\\in\\{1,3,5,10,22,47,100,216,465,1000\\}\n",
"$$\n",
"en fonction de $h$ en utilisant une échelle logarithmique pour les abscisses et les ordonnées. Vérifiez que l'on obtient une droite de pente égale à l'ordre du schéma (on pourra utiliser la fonction `reglin`)."
]
},
{
"cell_type": "code",
"execution_count": 66,
"id": "6a7463be-a872-4231-867a-1232205cdf99",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
""
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"e = [];\n",
"valN = [5,10,22,47,100,216,465,1000];\n",
"h = 4./valN;\n",
"for N = valN\n",
" t = linspace(0,4,N+1);\n",
" y = euler(f1,t,0);\n",
" e = [e, max(abs(y-ytrue(t)))];\n",
"end\n",
"\n",
"a = reglin(log(h),log(e));\n",
"plot(\"ll\",h,e,'o')\n",
"xlabel(\"h\")\n",
"ylabel(\"e(h)\")\n",
"title(sprintf(\"ordre experimental = %4.2f\",a))\n",
"e_euler = e;"
]
},
{
"cell_type": "markdown",
"id": "2c5c8822-346f-4b87-8553-902970a2d227",
"metadata": {},
"source": [
"6. Recommencez les deux points précédents pour les trois autres schémas."
]
},
{
"cell_type": "markdown",
"id": "da2e97f3-f263-4693-a214-f736769a2be9",
"metadata": {},
"source": [
"### Schéma d'Euler-Cauchy"
]
},
{
"cell_type": "code",
"execution_count": 67,
"id": "57fe3ec3-44ea-4c61-8743-2b6facefc969",
"metadata": {},
"outputs": [],
"source": [
"function y=euler_cauchy(f,t,y0)\n",
" n=length(t);\n",
" y=zeros(1,n);\n",
" h=t(2)-t(1);\n",
" y(1)=y0;\n",
" for k=1:n-1\n",
" k1 = f(t(k),y(k));\n",
" k2 = f(t(k)+h,y(k)+h*k1);\n",
" y(k+1) = y(k) + h/2*(k1+k2);\n",
" end\n",
"endfunction"
]
},
{
"cell_type": "code",
"execution_count": 68,
"id": "0bf27f5a-24c8-47b9-9c85-23a23db1f44f",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
""
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"for N = [10]\n",
" t = linspace(0,4,N+1);\n",
" u = linspace(0,4,500);\n",
" y = euler_cauchy(f1,t,0);\n",
" scf(); // cree une nouvelle fenetre graphique\n",
" plot(t,y,\"x\",u,ytrue(u))\n",
" legend(\"Euler\",\"$1-\\exp(-t^2/2)$\",4)\n",
" title(\"N=\"+string(N))\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 69,
"id": "14d9eed2-6633-4313-8719-82c957c42305",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
""
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"e = [];\n",
"for N = valN\n",
" t = linspace(0,4,N+1);\n",
" y = euler_cauchy(f1,t,0);\n",
" e = [e, max(abs(y-ytrue(t)))];\n",
"end\n",
"\n",
"a = reglin(log(h),log(e));\n",
"plot(\"ll\",h,e,'o')\n",
"xlabel(\"h\")\n",
"ylabel(\"e(h)\")\n",
"title(sprintf(\"ordre experimental = %4.2f\",a))\n",
"e_euler_cauchy = e;"
]
},
{
"cell_type": "markdown",
"id": "c5bba872-16d6-49e8-a816-2292c5cd0867",
"metadata": {},
"source": [
"### Point milieu"
]
},
{
"cell_type": "code",
"execution_count": 70,
"id": "ee6476d5-b8d1-4d0f-b683-96337618a9a9",
"metadata": {},
"outputs": [],
"source": [
"function y=point_milieu(f,t,y0)\n",
" n=length(t);\n",
" y=zeros(1,n);\n",
" h=t(2)-t(1);\n",
" y(1)=y0;\n",
" for k=1:n-1\n",
" k1 = f(t(k),y(k));\n",
" k2 = f(t(k)+h/2,y(k)+h/2*k1);\n",
" y(k+1) = y(k) + h*k2;\n",
" end\n",
"endfunction"
]
},
{
"cell_type": "code",
"execution_count": 71,
"id": "e38c5b0f-b654-46c2-af1d-885dd3329033",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
""
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"for N = [10]\n",
" t = linspace(0,4,N+1);\n",
" u = linspace(0,4,500);\n",
" y = point_milieu(f1,t,0);\n",
" scf(); // cree une nouvelle fenetre graphique\n",
" plot(t,y,\"x\",u,ytrue(u))\n",
" legend(\"Point milieu\",\"$1-\\exp(-t^2/2)$\",4)\n",
" title(\"N=\"+string(N))\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 72,
"id": "e0b6a61f-b70a-4f4b-89c0-7bfb13d106bc",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
""
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"e = [];\n",
"for N = valN\n",
" t = linspace(0,4,N+1);\n",
" y = point_milieu(f1,t,0);\n",
" e = [e, max(abs(y-ytrue(t)))];\n",
"end\n",
"\n",
"a = reglin(log(h),log(e));\n",
"plot(\"ll\",h,e,'o')\n",
"xlabel(\"h\")\n",
"ylabel(\"e(h)\")\n",
"title(sprintf(\"ordre experimental = %4.2f\",a))\n",
"e_point_milieu = e;"
]
},
{
"cell_type": "markdown",
"id": "25b50f74-7111-4763-94e3-3e7eee2146bb",
"metadata": {},
"source": [
"### Runge-Kutta d'ordre 4"
]
},
{
"cell_type": "code",
"execution_count": 73,
"id": "13f87429-0eb7-4454-9f88-88752047bffc",
"metadata": {},
"outputs": [],
"source": [
"function y=rk4(f,t,y0)\n",
" n=length(t);\n",
" y=zeros(1,n);\n",
" h=t(2)-t(1);\n",
" y(1)=y0;\n",
" for k=1:n-1\n",
" k1 = f(t(k),y(k));\n",
" k2 = f(t(k)+h/2,y(k)+h/2*k1);\n",
" k3 = f(t(k)+h/2,y(k)+h/2*k2);\n",
" k4 = f(t(k)+h,y(k)+h*k3);\n",
" y(k+1) = y(k)+h/6*(k1+2*k2+2*k3+k4);\n",
" end\n",
"endfunction"
]
},
{
"cell_type": "code",
"execution_count": 74,
"id": "667409dc-63df-400b-98ca-5cfcd9edca5d",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
""
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"for N = [10]\n",
" t = linspace(0,4,N+1);\n",
" u = linspace(0,4,500);\n",
" y = rk4(f1,t,0);\n",
" scf(); // cree une nouvelle fenetre graphique\n",
" plot(t,y,\"x\",u,ytrue(u))\n",
" legend(\"RK4\",\"$1-\\exp(-t^2/2)$\",4)\n",
" title(\"N=\"+string(N))\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 75,
"id": "038d8fb7-66e3-4b94-91aa-560fa8ea414a",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
""
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"e = [];\n",
"for N = valN\n",
" t = linspace(0,4,N+1);\n",
" y = rk4(f1,t,0);\n",
" e = [e, max(abs(y-ytrue(t)))];\n",
"end\n",
"\n",
"a = reglin(log(h),log(e));\n",
"plot(\"ll\",h,e,'o')\n",
"xlabel(\"h\")\n",
"ylabel(\"e(h)\")\n",
"title(sprintf(\"ordre experimental = %4.2f\",a))\n",
"e_rk4 = e;"
]
},
{
"cell_type": "markdown",
"id": "00d816b6-ff52-43db-942e-a95daf9f6096",
"metadata": {},
"source": [
"7. Pour terminer, superposer sur un même graphique les graphes des erreurs $h\\rightarrow e(h)$ (pour les 4 schémas) ainsi que l'approximation de $y(t)$ renvoyée par la macro `cvode` de Scilab. Cherchez dans la page d'aide de `cvode` comment il est possible de contrôler cette précision."
]
},
{
"cell_type": "code",
"execution_count": 108,
"id": "6f08a816-062c-49a0-966d-b3da35cde044",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
""
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"e = [];\n",
"for N = valN\n",
" t = linspace(0,4,N+1);\n",
" [_,y] = cvode(f1,t,0,atol=1e-14,rtol=1e-14);\n",
" e = [e, max(abs(y-ytrue(t)))];\n",
"end\n",
"e_cvode = e;\n",
"plot(\"ll\",h,[e_euler;e_euler_cauchy;e_point_milieu;e_rk4;e_cvode],\"-o\");\n",
"legend(\"Euler\",\"Euler Cauchy\",\"Point milieu\",\"RK4\",\"cvode\",-1);\n",
"xlabel(\"h\")\n",
"ylabel(\"e(h)\")"
]
},
{
"cell_type": "markdown",
"id": "dfa81591-b394-424a-a83b-49f0d002c82c",
"metadata": {},
"source": [
"## Application : le pendule\n",
"\n",
"
\n",
"$\\theta(t)$\n",
"
\n",
"Le pendule que nous allons considérer est représenté sur la figure ci-dessus. Nous supposerons que la tige reliant un poids de masse $M$ à l'axe de rotation est de masse négligeable devant $M$. Nous choisissons de mesurer la déviation du pendule de la position verticale d'équilibre stable par l'angle $\\theta(t)$ mesuré positivement comme défini sur la figure."
]
},
{
"cell_type": "markdown",
"id": "0779f58a-5dd7-4d38-a8c9-a93f556556f9",
"metadata": {},
"source": [
"Si l'on applique les relations de la dynamique pour les solides en rotation autour d'un axe, on obtient l'équation différentielle suivante :\n",
"\\begin{equation}\n",
"\\left\\{\\begin{array}{rcl}\n",
"\\theta''(t) &=& - \\displaystyle\\frac{g}{L}\\sin\\theta(t),\\\\\n",
"\\theta(0)&=&\\theta_0,\\\\\n",
"\\theta'(0)&=&0.\n",
"\\end{array}\\right.\n",
"\\tag{3}\n",
"\\end{equation}\n",
"La valeur de $\\theta_0$ donne la déviation initiale du pendule, et on a considéré que la vitesse angulaire initiale est nulle. \n"
]
},
{
"cell_type": "markdown",
"id": "40979f0b-0d35-4501-a720-72c13f6b37b4",
"metadata": {},
"source": [
"1. Montrer que l'on peut approcher (sous certaines conditions que l'on précisera) la solution de cette équation différentielle par la fonction\n",
"$$\n",
"\\phi(t) = \\theta_0\\cos\\left(\\sqrt{\\frac{g}{L}}t\\right)\n",
"$$\n",
"que l'on appellera dans la suite \"solution linéarisée\"."
]
},
{
"cell_type": "markdown",
"id": "9eeb83fd-ee24-4ff5-bd50-a01236cf28e5",
"metadata": {},
"source": [
"Si on suppose que $\\theta_0$ est petit, on peut approcher $\\sin(\\theta)$ par son développement à l'ordre 1, ce qui donne l'équation différentielle\n",
"\\begin{equation}\n",
"\\left\\{\\begin{array}{rcl}\n",
"\\phi''(t) &=& - \\displaystyle\\frac{g}{L}\\phi(t),\\\\\n",
"\\phi(0)&=&\\theta_0,\\\\\n",
"\\phi'(0)&=&0,\n",
"\\end{array}\\right.\n",
"\\tag{4}\n",
"\\end{equation}\n",
"dont on cherche les solutions sous la forme $\\phi(t)=a\\cos(\\omega_0t)+b\\sin(\\omega_0t)$, où $\\omega_0=\\sqrt{g/L}$. Les deux conditions initiales permettent de trouver $b=0$ et $a=\\theta_0$."
]
},
{
"cell_type": "markdown",
"id": "553bf165-d622-40e4-9755-d8f37807d160",
"metadata": {},
"source": [
"2. On pose $y=\\left(\\theta,\\theta'\\right)$, \n",
"écrire (3) sous la forme d'un système de deux équations du premier ordre en temps $y'=f(t,y)$ et écrire la fonction Scilab calculant $f$ :"
]
},
{
"cell_type": "code",
"execution_count": 128,
"id": "78b8f050-269b-4fa9-8fab-8285ed48ec39",
"metadata": {},
"outputs": [],
"source": [
"function dydt=fpendule(t,y)\n",
" g=9.81;\n",
" L=1;\n",
" theta=y(1);\n",
" thetaprime=y(2);\n",
" dydt = [thetaprime; -g/L*sin(theta)];\n",
"endfunction"
]
},
{
"cell_type": "markdown",
"id": "40bee13a-ab5f-4b33-a9f1-77623f76edf6",
"metadata": {},
"source": [
"3. Superposer sur un même graphique la solution linéarisée et l'approximation de $\\theta(t)$ donnée par Scilab avec la fonction `cvode` pour les valeurs de $\\theta_0$ dans l'ensemble\n",
"$$\n",
"\\left\\{\\frac{\\pi}{32},~\\frac{\\pi}{16},~\\frac{\\pi}{8},~\\frac{\\pi}{2},~\\frac{3\\pi}{4},~\\frac{63\\pi}{64}\\right\\}\\cdot\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 124,
"id": "fcba8f8e-35b9-40aa-bf19-05f5f3c87774",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
""
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"t = linspace(0,10,500);\n",
"theta = [];\n",
"for theta0 = %pi./[32,16,8,2,4/3,64/63]\n",
" [_,y] = cvode(fpendule,t,[theta0;0],rtol=1e-6);\n",
" theta = [theta;y(1,:)];\n",
"end\n",
"plot(t,theta)\n",
"xlabel(\"t\")\n",
"ylabel(\"$\\theta(t)$\")"
]
},
{
"cell_type": "markdown",
"id": "31fa4a35-aa33-4436-8082-32fae9989bde",
"metadata": {},
"source": [
"On peut aussi représenter les courbes paramétriques $(\\theta(t),\\theta'(t))$ dans le plan de phase :"
]
},
{
"cell_type": "code",
"execution_count": 127,
"id": "1a1ae580-3d42-45f8-9f16-7b662f239006",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
""
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"t = linspace(0,10,500);\n",
"theta = [];\n",
"for theta0 = %pi./[32,16,8,2,4/3,64/63]\n",
" [_,y] = cvode(fpendule,t,[theta0;0],rtol=1e-6);\n",
" plot(y(1,:),y(2,:))\n",
"end\n",
"xlabel(\"$\\theta$\")\n",
"ylabel(\"$\\theta$\")\n",
"title(\"Portrait de phase\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "731d9be3-f022-402d-aeb5-2c0e5b135a0f",
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Scilab",
"language": "scilab",
"name": "scilab"
},
"language_info": {
"file_extension": ".sci",
"help_links": [
{
"text": "MetaKernel Magics",
"url": "https://metakernel.readthedocs.io/en/latest/source/README.html"
}
],
"mimetype": "text/x-scilab",
"name": "scilab",
"version": "0.10.2"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
![](\"thependule.svg\")
$\\theta(t)$\n",
"