diff --git a/examples/Demonstration-Physics.ipynb b/examples/Demonstration-Physics.ipynb index d47105e..faebdbf 100644 --- a/examples/Demonstration-Physics.ipynb +++ b/examples/Demonstration-Physics.ipynb @@ -1,654 +1,645 @@ { "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", " Example extracted from the Mécanique Générale repository (MecaDRIL) by Cécile Hébert, more on: https://github.com/c-hebert/MecaDRIL.\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# EXPERIENCE DE LA BILLE SUR GLISSIÈRE" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import urllib\n", "import os\n", - "from notebook import notebookapp\n", + "#from notebook import notebookapp\n", + "from jupyter_server import serverapp\n", "\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "import pandas as pd\n", "import math\n", "\n", "import io\n", "import base64\n", "import ipywidgets as widgets\n", "from ipywidgets import interact, interactive, fixed, interact_manual, Label, VBox\n", "from IPython.display import HTML, Video, IFrame\n", "\n", "from bokeh.plotting import figure, curdoc\n", "from bokeh.models.widgets import Slider, CheckboxButtonGroup, PreText\n", "from bokeh.layouts import row, column, widgetbox\n", "from bokeh.models import ColumnDataSource, Slider, Button, TextInput, Arrow, OpenHead, NormalHead, VeeHead\n", "from bokeh.plotting import figure, show, ColumnDataSource\n", "\n", "from bokeh.io import output_notebook, show, export_png\n", "\n", "output_notebook()\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Le problème :" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Une bille est lâchée dans une glissière qui forme un looping. Si elle n'est pas lâchée d'une hauteur suffisante, elle ne complète pas le looping. De quelle hauteur faut-il la lâcher pour qu'elle complète le looping ?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Détail des calculs et de la modélisation" ] }, { - "cell_type": "code", - "execution_count": null, + "cell_type": "markdown", "metadata": {}, - "outputs": [], "source": [ - "from IPython.display import Image\n", - "Image(\"figs/Glissiere.png\")" + "\"expérience" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "On lâche la bille d'une hauteur h et on observe son comportement au point P en supposant qu'elle suit jusqu'en P la trajectoire imposée par la glissière. On prend $P$ dans le quart supérieur droit de la glissière ($\\theta >0$)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 1. On cherche la norme de $\\vec v$ au point P: $v_p$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "La conservation de l'énergie donne: \n", "$$mgh = mgy_{p} + \\frac{1}{2}mv_{p}^{2}$$\n", "\n", "$$y_p = R ( 1 + \\sin\\theta)$$\n", "\n", "$$v_p = \\sqrt{2g[h - R ( 1 + \\sin\\theta)]}$$\n", "\n", "Un calcul plus juste tient compte de l'énergie de rotation de la bille. Dans ce cas la conservation de l'énergie devient: \n", "$$mgh = mgy_p + \\frac{1}{2}mv_{p}^{2} + \\frac{1}{2}Iw_{p}^{2}$$\n", "\n", "$$I = \\frac{2}{5}mr^2$$\n", "\n", "$$w_p = \\frac{v_p}{r}$$ avec r le rayon de la bille \n", "\n", "$$mgh = mgy_p + \\frac{1}{2}mv_{p}^{2} + \\frac{1}{2}\\frac{2}{5}mr^2\\frac{v_p^2}{r^2}$$\n", "\n", "$$gh = gR(1 + \\sin\\theta) + \\frac{7}{10}v_p^2$$\n", "\n", "$$v_p = \\sqrt{\\frac{10}{7}g[h - R(1 + \\sin\\theta)]}$$\n", "\n", "de manière générale, on peut donc écrire\n", "\n", "$$v_p = \\sqrt{A g[h - R(1 + \\sin\\theta)]}$$\n", "\n", "Avec $A$ un facteur qui vaut 2 pour un objet qui glisse et $\\frac{10}{7}$ pour une sphère homogène qui roule sans glisser.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 2. On cherche les composantes du vecteur $\\vec v_p$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "La direction du vecteur vitesse en P s'obtient par le fait qu'il est tangent à la glissière\n", "\n", "$$\\vec v=\\begin{pmatrix}-v_p\\sin\\theta\\\\v_p\\cos\\theta \\end{pmatrix}$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 3. On calcule la trajectoire parabolique qu'aurait une bille libre dans le champ de pesanteur si elle était lancée en P avec $\\vec v_p$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$\\vec r_P=\\begin{pmatrix}-v_p\\sin\\theta \\, t + R \\cos \\theta \\\\-\\frac{1}{2} g t^2 + v_p\\cos\\theta \\, t \n", "+ R (1+\\sin \\theta)\n", "\\end{pmatrix}$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 4. Analyse des forces exercées sur la bille en P." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Les forces sont le poids $\\vec P$ et la réaction de la glissière $\\vec N$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$$ \\vec P = m \\vec g = - mg \\vec e_y $$ \n", "\n", "$$ \\vec N = N \\vec n$$ \n", "avec $\\vec n$ vecteur normal au support pointant vers le centre du cercle. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "$\\vec N$ est forcément dirigée vers le centre du cercle donc $N \\geq 0$. \n", "\n", "$N <0$ impliquerait que la glissière attire la bille. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Au point P, \n", "$$ \\sum \\vec F = m \\vec a = m \\frac{v_p^2}{R}\\vec n + m \\left ( \\frac{\\mathrm{d} v}{\\mathrm{d} t} \\right )\\vec \\tau = -mg\\vec e_y + N\\vec n$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Projeté sur $\\vec n $, on obtient, \n", "\n", "$$m g \\sin \\theta + N = m \\frac{v_P^2}{R}$$\n", "\n", "$$N= m(\\frac{v_P^2}{R} - g \\sin \\theta)$$\n", "\n", "$$N= m g \\Big(A \\frac{h}{R} - A (1 + \\sin\\theta) - \\sin \\theta \\Big)$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "On remarque donc que si $h$ n'est pas assez grand, $N$ va devenir négatif pour un $\\theta < \\frac{\\pi}{2}$. Cela signifie que au point P, la bille devrait être attiré par la glissière pour rester dessus... elle décolle. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Tracé de la figure " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Paramètres pour la figure" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Définition de l' image de fond\n", "import base64\n", "\n", "def file2src(filename, data_type='image/png'):\n", " with open(filename, \"rb\") as image_file:\n", " data = base64.b64encode(image_file.read())\n", " image_file.close()\n", " return 'data:{};base64,{}'.format(data_type, data.decode())\n", "\n", "# Définition de l'image de fond\n", "img_path = file2src(\"resources/background.png\")\n", "\n", "# Import des données du video\n", "df = pd.read_csv('resources/data-05-04-19.csv', index_col=0)\n", "df = df.rename({'VideoAnalysis: X (cm)': 'x', 'VideoAnalysis: Y (cm)':'y',\n", " 'VideoAnalysis: X Velocity (cm/s)':'v_x', 'VideoAnalysis: Y Velocity (cm/s)':'v_y'}, axis='columns')\n", "\n", "# Calculer les limites \n", "y_min = df['y'].min()\n", "x_min = df['x'].min()\n", "x_max = df['x'].max()\n", "x_lowest = df.loc[df['y'] == y_min]['x'].iloc[0]\n", "\n", "# Calculer le rayon de la boucle\n", "df_h = df.loc[df['x'] < 20]\n", "df_h = df_h.loc[df_h['x'] > -20]\n", "R = df_h['y'].max() - df_h['y'].min()\n", "idx = df_h.loc[df_h['y'] == df_h['y'].min()]\n", "R = R/2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Équations et paramètres initiaux" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Choissir A: 2 ou 10/7\n", "A = 10/7\n", "\n", "# Définir parametres initiaux\n", "g = 9.81 # pesanteur\n", "h = 2*R # diamètre\n", "theta = 20 # angle par rapport à l'horizontale\n", "to_rad = lambda x: np.pi*x/180.0\n", "thetarad = to_rad(theta)\n", "\n", "t = np.linspace(0, 2, 100) # temps après que la bille est arrivéee en P (tracé parabole)\n", "neg_t = np.linspace(0, -1, 100) # temps avant que la bille soit arrivéee en P (tracé parabole)\n", "\n", "# Eq. parametrique de la parabole\n", "x_p = lambda v, ang, t: -v*np.sin(ang)*t + R*np.cos(ang)\n", "y_p = lambda v, ang, t: -0.5*g*t*t + v*np.cos(ang)*t + R*(1 + np.sin(ang))\n", "\n", "# Définition du vecteur vitesse (tracé depuis P)\n", "vx = lambda v, ang: np.linspace(x_p(v, ang, t)[0], x_p(v, ang, t)[0] - 0.5*v*np.sin(ang), 10) \n", "vy = lambda v, ang: np.linspace(y_p(v, ang, t)[0], y_p(v, ang, t)[0] + 0.5*v*np.cos(ang), 10) \n", "\n", "# Define equations\n", "vitesse_p = lambda h, A, theta_rad: math.sqrt(A*g*R*(h - (1 + np.sin(theta_rad)))) #Vitesse au point A\n", "normale_p = lambda hsurR, A, angle: g*(A*hsurR-A*(1+np.sin(angle))-np.sin(angle)) #Force normale au point A\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Figure" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "## PLOT\n", "\n", "def modify_doc(doc):\n", " \n", " # Définir information à visualiser\n", " x = -df['x']\n", " y = df['y'] - y_min\n", " x0 = x.iloc[0]\n", " y0 = y.iloc[0]\n", "\n", " # Définir glissière\n", " slope2 = (y.iloc[6] - y.iloc[5]) / (x.iloc[6] - x.iloc[5])\n", " get_x = lambda y_ini: x0 - (y0 - y_ini)*(1/slope2)\n", "\n", " line_ini_y = np.linspace(2, y.iloc[0], 20)\n", " line_ini_x = get_x(line_ini_y)\n", "\n", " slope = np.arctan((y.iloc[6] - y.iloc[5]) / (x.iloc[6] - x.iloc[5]))\n", " line_end_x = x.iloc[-1] - np.linspace(0, 40, 20)*np.cos(-slope-0.02)\n", " line_end_y = y.iloc[-1] - np.linspace(0, 40, 20)*np.sin(-slope -0.02)\n", "\n", " angle = np.linspace(0, 2*np.pi, 100)\n", " \n", " # Conditions initiales \n", " g = 9.81\n", " h = 2.0\n", " A = 10/7\n", " theta = 45\n", " thetarad = to_rad(theta)\n", " t = np.linspace(0, 2, 100)\n", " neg_t = np.linspace(0, -1, 100)\n", " vit_P = vitesse_p(h,A, thetarad)\n", " lwidth = 2.5 \n", "\n", "\n", " def initialize(h_init, thetarad):\n", " # Calculer position, vitesse et forces au point P\n", " vit_P = vitesse_p(h_init, A, thetarad) # vitesse en P\n", " n_P = normale_p(h_init, A, thetarad) # force normale en P\n", " \n", " # coordonnées parabole après point P\n", " x_parab = x_p(vit_P, thetarad, t) \n", " y_parab = y_p(vit_P, thetarad, t)\n", " \n", " # coordonnées parabole avant point P\n", " x_parab_neg = x_p(vit_P, thetarad, neg_t)\n", " y_parab_neg = y_p(vit_P, thetarad, neg_t)\n", " \n", " # vecteur force normale\n", " norm_end_x = x_parab[0] - n_P*np.cos(thetarad)\n", " norm_end_y = y_parab[0] - n_P*np.sin(thetarad)\n", " \n", " # vecteur vitesse\n", " v_x, v_y = vx(vit_P, thetarad), vy(vit_P, thetarad)\n", " \n", " P = {\n", " 'x_init': [get_x(h_init*R)], 'y_init':[h_init*R],\n", " 'x_point_P': [x_parab[0]], 'y_point_P': [y_parab[0]]\n", " }\n", " \n", " CURV = {\n", " 'x_point_P_parab_pos': x_parab, 'y_point_P_parab_pos': y_parab,\n", " 'x_point_P_parab_neg': x_parab_neg, 'y_point_P_parab_neg': y_parab_neg\n", " }\n", " \n", " ANG = {\n", " 'l_hor_x':np.linspace(0,R, 100), 'l_hor_y':[R]*100, \n", " 'l_ver_x': np.linspace(0, x_parab[0], 100), 'l_ver_y': np.linspace(R,y_parab[0],100)\n", " \n", " }\n", " \n", " FORCES = {\n", " 'grav_x': [x_parab[0]]*10, 'grav_y': np.linspace(y_parab[0], -9.8+y_parab[0], 10),\n", " 'arrow_up_grav_x': np.linspace(x_parab[0], 0.5 + x_parab[0], 10), \n", " 'arrow_up_grav_y': np.linspace(-9.8 + y_parab[0], 0.6 -9.8 + y_parab[0], 10),\n", " 'arrow_down_grav_x': np.linspace(x_parab[0], -0.5 + x_parab[0], 10), \n", " 'arrow_down_grav_y': np.linspace(-9.8 + y_parab[0], 0.6 -9.8 + y_parab[0], 10),\n", " 'norm_x': np.linspace(x_parab[0], norm_end_x, 10), 'norm_y': np.linspace(y_parab[0], norm_end_y, 10),\n", " 'arrow_up_norm_x': norm_end_x + np.linspace(0, 0.1*n_P*np.cos(thetarad - np.pi*20/180), 10), \n", " 'arrow_up_norm_y': norm_end_y + np.linspace(0, 0.1*n_P*np.sin(thetarad - np.pi*20/180), 10), \n", " 'arrow_down_norm_x': norm_end_x + np.linspace(0, 0.1*n_P*np.cos(thetarad + np.pi*20/180), 10), \n", " 'arrow_down_norm_y': norm_end_y + np.linspace(0, 0.1*n_P*np.sin(thetarad + np.pi*20/180), 10)\n", " }\n", " \n", " VEL = {\n", " 'v_x': v_x,\n", " 'v_y': v_y, \n", " 'arrow_up_vel_x': v_x[-1] + np.linspace(0, 0.1*vit_P*np.cos(thetarad - np.pi*70/180), 10),\n", " 'arrow_up_vel_y': v_y[-1] + np.linspace(0, 0.1*vit_P*np.sin(thetarad - np.pi*70/180), 10), \n", " 'arrow_down_vel_x': v_x[-1] + np.linspace(0,- 0.1*vit_P*np.cos(thetarad + np.pi*70/180), 10),\n", " 'arrow_down_vel_y': v_y[-1] + np.linspace(0,- 0.1*vit_P*np.sin(thetarad + np.pi*70/180), 10) \n", "\n", " }\n", " \n", " return P, CURV, ANG, FORCES, VEL\n", " \n", " \n", " # Calculer les valeurs initiales\n", " P, C, ANG, FORCES, VEL = initialize(h, thetarad)\n", "\n", " # Créer dictionnaire pour acceder aux données\n", " source_P = ColumnDataSource(data = P)\n", " source_C = ColumnDataSource(data = C)\n", " source_ANGLE = ColumnDataSource(data = ANG)\n", " source_FORCES = ColumnDataSource(data = FORCES)\n", " source_VITESSE = ColumnDataSource(data= VEL)\n", " \n", " \n", " # Figure\n", " p = figure(title=\"Trajectory\", plot_height=425, plot_width=950, x_range=(-57,53), y_range=(-5,45),\\\n", " background_fill_color='#ffffff')\n", " \n", " # Définir image de fond\n", " p.image_url(url=[img_path], x=-56.5, y=45, w=110, h=45, alpha=0.2, angle=-np.pi*1/180)\n", "\n", "\n", " # Partie fixe de la figure:\n", " # Glissière\n", - " p.line(line_ini_x, line_ini_y, color='#635c74', alpha=0.7, line_width=2.5, legend='Piste')\n", - " p.line(R*np.sin(angle), R + R*np.cos(angle), color='#635c74', alpha=0.6, line_width=2.5, legend='Piste')\n", + " p.line(line_ini_x, line_ini_y, color='#635c74', alpha=0.7, line_width=2.5, legend_label='Piste')\n", + " p.line(R*np.sin(angle), R + R*np.cos(angle), color='#635c74', alpha=0.6, line_width=2.5, legend_label='Piste')\n", " # Centre de la boucle\n", - " p.circle(0,R,size=3, fill_color='#635c74', line_color='#635c74', legend='Centre')\n", + " p.circle(0,R,size=3, fill_color='#635c74', line_color='#635c74', legend_label='Centre')\n", "\n", "\n", " # Variables:\n", " # Point initial\n", - " p_ini = p.circle('x_init', 'y_init', source=source_P, size=6, fill_color='#e32020', line_color='#e32020', legend='Point Initial')\n", + " p_ini = p.circle('x_init', 'y_init', source=source_P, size=6, fill_color='#e32020', line_color='#e32020', legend_label='Point Initial')\n", "\n", " # Point A\n", " p_a = p.circle('x_point_P','y_point_P', source=source_P, size=6, \\\n", - " legend='A', line_color='#e32020', fill_color='#e32020')\n", + " legend_label='A', line_color='#e32020', fill_color='#e32020')\n", "\n", " # Parabole\n", " parab = p.line('x_point_P_parab_pos', 'y_point_P_parab_pos', source=source_C, color='#e32020', line_width=1.5, alpha=0.8, \\\n", - " legend='Parabole')\n", + " legend_label='Parabole')\n", "\n", " parab_neg = p.line('x_point_P_parab_neg', 'y_point_P_parab_neg', source=source_C, color='#e32020', line_width=1.5, alpha=0.8, \\\n", - " legend='Parabole')\n", + " legend_label='Parabole')\n", " \n", "\n", " # Angle theta\n", " l_h = p.line('l_hor_x', 'l_hor_y', source=source_ANGLE, line_width=lwidth, color='#000000', alpha=0.2)\n", " l_v = p.line('l_ver_x', 'l_ver_y', source=source_ANGLE, line_width=lwidth ,color='#000000', alpha=0.2)\n", "\n", " # Force gravité\n", - " f_g = p.line('grav_x', 'grav_y', source=source_FORCES, line_width=lwidth, color='#178717', alpha=0.8, legend='mg')\n", - " a_u_f = p.line('arrow_up_grav_x', 'arrow_up_grav_y', source=source_FORCES, line_width=lwidth, color='#178717', alpha=0.8, legend='mg')\n", - " a_d_f = p.line('arrow_down_grav_x', 'arrow_down_grav_y', source=source_FORCES, line_width=lwidth, color='#178717', alpha=0.8, legend='mg')\n", + " f_g = p.line('grav_x', 'grav_y', source=source_FORCES, line_width=lwidth, color='#178717', alpha=0.8, legend_label='mg')\n", + " a_u_f = p.line('arrow_up_grav_x', 'arrow_up_grav_y', source=source_FORCES, line_width=lwidth, color='#178717', alpha=0.8, legend_label='mg')\n", + " a_d_f = p.line('arrow_down_grav_x', 'arrow_down_grav_y', source=source_FORCES, line_width=lwidth, color='#178717', alpha=0.8, legend_label='mg')\n", "\n", " # Force normale\n", - " f_n = p.line('norm_x', 'norm_y', source=source_FORCES, line_width=lwidth*1.2, color='#064d06', alpha=0.8, legend='N')\n", - " a_u_n = p.line('arrow_up_norm_x', 'arrow_up_norm_y', source=source_FORCES, line_width=lwidth, color='#064d06', alpha=0.8, legend='N')\n", + " f_n = p.line('norm_x', 'norm_y', source=source_FORCES, line_width=lwidth*1.2, color='#064d06', alpha=0.8, legend_label='N')\n", + " a_u_n = p.line('arrow_up_norm_x', 'arrow_up_norm_y', source=source_FORCES, line_width=lwidth, color='#064d06', alpha=0.8, legend_label='N')\n", "\n", - " a_d_n = p.line('arrow_down_norm_x', 'arrow_down_norm_y', source=source_FORCES, line_width=lwidth, color='#064d06', alpha=0.8, legend='N')\n", + " a_d_n = p.line('arrow_down_norm_x', 'arrow_down_norm_y', source=source_FORCES, line_width=lwidth, color='#064d06', alpha=0.8, legend_label='N')\n", "\n", " # Vitesse\n", - " v = p.line('v_x', 'v_y', source=source_VITESSE, color='#01B9FF', legend='Vitesse', line_width=lwidth) \n", + " v = p.line('v_x', 'v_y', source=source_VITESSE, color='#01B9FF', legend_label='Vitesse', line_width=lwidth) \n", " a_u = p.line('arrow_up_vel_x','arrow_up_vel_y', source=source_VITESSE, line_width=lwidth)\n", " a_d = p.line('arrow_down_vel_x','arrow_down_vel_y', source=source_VITESSE, line_width=lwidth)\n", " \n", " \n", " # Définir sliders \n", " slider_h = Slider(start=2, end=3, step=0.05, value=2, title='Hauteur initiale')\n", " slider_om = Slider(start=0, end=90, step=1, value=45, title='Theta')\n", " \n", " pre = PreText(text='Point initiale.\\tx:{:3.2f} \\ty:{:3.2f} \\nPoint A.\\tx: {:3.2f} \\ty: {:3.2f} \\nVitesse.\\t{:3.2f}[cm/s]\\tx:{:3.2f}[cm/s]\\ty:{:3.2f}[cm/s] \\nForce normale. \\t{:3.2f} [N]'\\\n", " .format(source_P.data['x_init'][0], source_P.data['y_init'][0],source_P.data['x_point_P'][0], source_P.data['y_point_P'][0], \\\n", " vit_P, -vit_P*np.cos(thetarad), vit_P*np.sin(thetarad), normale_p(h, A, thetarad)),\n", " width=500, height=100)\n", " \n", " \n", " def refresh_source(attrname, old, new):\n", " \n", " # Mettre à jour les valeurs à afficher \n", " om = slider_om.value\n", " h = slider_h.value\n", " #h = h*R # hauteur par rapport au rayon\n", " omrad = to_rad(om)\n", "\n", " P, C, ANG, FORCES, VEL = initialize(h, omrad)\n", " \n", " source_P.data = P\n", " source_C.data = C\n", " source_ANGLE.data = ANG\n", "\n", " source_FORCES.data = FORCES\n", " source_VITESSE.data = VEL\n", " \n", " vit_P = vitesse_p(h, A, omrad)\n", " x_p = source_C.data['x_point_P_parab_pos'][0]\n", " y_p = source_C.data['y_point_P_parab_pos'][0]\n", " \n", " pre.text='Point initiale.\\tx:{:3.2f} \\ty:{:3.2f} \\nPoint A.\\tx: {:3.2f} \\ty: {:3.2f} \\nVitesse.\\t{:3.2f}[cm/s]\\tx:{:3.2f}[cm/s]\\ty:{:3.2f}[cm/s] \\nForce normale. \\t{:3.2f} [N]'\\\n", " .format(source_P.data['x_init'][0], source_P.data['y_init'][0],source_P.data['x_point_P'][0], source_P.data['y_point_P'][0], \\\n", " vit_P, -vit_P*np.sin(omrad), vit_P*np.cos(omrad), normale_p(h, A, omrad))\n", " \n", " \n", " slider_h.on_change('value', refresh_source)\n", " slider_om.on_change('value', refresh_source)\n", " \n", " layout = column(\n", " row(p),\n", " row(slider_h, slider_om),\n", " row(pre)\n", " )\n", "\n", " \n", " doc.add_root(layout)\n", " " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def remote_jupyter_proxy_url(port):\n", " \"\"\"\n", " Callable to configure Bokeh's show method when a proxy must be\n", " configured.\n", "\n", " If port is None we're asking about the URL\n", " for the origin header.\n", " \"\"\"\n", " \n", " base_url = os.environ['EXTERNAL_URL']\n", " host = urllib.parse.urlparse(base_url).netloc\n", "\n", " # If port is None we're asking for the URL origin\n", " # so return the public hostname.\n", " if port is None:\n", " return host\n", "\n", " service_url_path = os.environ['JUPYTERHUB_SERVICE_PREFIX']\n", " proxy_url_path = 'proxy/%d' % port\n", "\n", " user_url = urllib.parse.urljoin(base_url, service_url_path)\n", " full_url = urllib.parse.urljoin(user_url, proxy_url_path)\n", " return full_url\n", "\n", "\n", "\n", "def show_document(doc):\n", - " servers = list(notebookapp.list_running_servers())[0]\n", + " servers = list(serverapp.list_running_servers())[0]\n", " if servers['hostname'] == 'localhost':\n", " show(doc) \n", " else:\n", " show(doc, notebook_url=remote_jupyter_proxy_url)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Utilisation de la figure" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "On varie h et $\\theta$ et on observe: \n", "\n", "- l'intensité et la direction de $\\vec N$, calculé par la relation $N= m g \\Big(A \\frac{h}{R} - A (1 + \\sin\\theta) - \\sin \\theta \\Big)$\n", "\n", "- et l'allure de la parabole qu'aurait une bille en chute libre avec $\\vec v_p$ en P. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "***Comment se place la parabole par rapport à la glissière quand N=0 ? Et quand N>0 ? Et N<0 ?***" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ - "show_document(modify_doc) " + "show_document(modify_doc);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Feedback" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "IFrame('https://www.surveymonkey.com/r/NOTOSURVEY?notebook_set=MecaDRIL¬ebook_id=bille_sur_glissiere', 600, 800)" ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] } ], "metadata": { "hide_input": false, "kernelspec": { - "display_name": "Python 3", + "display_name": "Python", "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.6.9" + "version": "3.8.10" } }, "nbformat": 4, "nbformat_minor": 4 } diff --git a/examples/ExercisesA-Structures-ExerciseSheet.ipynb b/examples/ExercisesA-Structures-ExerciseSheet.ipynb index 7ca5a80..dc5bd61 100644 --- a/examples/ExercisesA-Structures-ExerciseSheet.ipynb +++ b/examples/ExercisesA-Structures-ExerciseSheet.ipynb @@ -1,560 +1,560 @@ { "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", " Example extracted from the MNSS repository by Guillaume Anciaux, more on: https://gitlab.epfl.ch/anciaux/mnss-notebooks.\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", " CIVIL-321 \"Modélisation Numérique des Solides et Structures\"

\n", " Comment utiliser ce Jupyter Notebook?\n", "

\n", " Ce Notebook est constitué de cellules de texte et de cellule de code. Les cellules de codes doivent être executées pour voir le résultat du programme. Certaines cellules doivent être remplies par vos soins. Pour exécuter une cellule, cliquez dessus simplement et ensuite cliquez sur le bouton \"play\" (►) dans la barre de menu au dessus du notebook. Vous pouvez aussi taper la combinaison de touches shift + enter. Il est important d'éxécuter les cellules de code en respectant leur ordre d'arrivée dans le notebook.\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "
\n", " Vous avez une question ?

\n", " SpeakUp: https://speakup.epfl.ch/ng/room/05096\n", "
" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from lib.plot import *" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Série d'exercices : Systèmes de ressorts" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Introduction\n", "\n", "Cette séance d'exercice contient des exercices à faire à la main ainsi que par programmation sur Python. Tout peut être fait dans ce *Notebook*. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## À lire avant de commencer\n", "\n", "\n", "- Les index sur Python commence à 0 (et non à 1 **contrairement à Matlab**). \n", "- Les matrices sont crées via le module *Numpy* (https://numpy.org/doc/stable/reference/arrays.ndarray.html)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# exemple, creation d'une matrice \n", "import numpy as np\n", "A = np.array([[1,2],[3, 4]])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- La fonction `plot_matrix(matrix, matrix_name)` permet de visualiser les matrices. \n", " Elle prend en entrée: \n", " - matrix : le nom de la variable \n", " - matrix_name = 'nom_matrix'" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# exemple d'affichage\n", "plot_matrix(A, 'A')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- l'extraction de bloc d'un vecteur ou d'une matrice se fait via les `slice` (https://numpy.org/doc/stable/reference/arrays.indexing.html)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# exraction d'une ligne\n", "plot_matrix(A[0, :], 'A_{{0, x}}')\n", "\n", "# exraction d'une colonne\n", "plot_matrix(A[:, 1], 'A_{{x, 1}}')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exercice 1: Assemblage de matrice (à la main)\n", "\n", "![](figs/ex1.svg)\n", "\n", "*Système de ressorts considéré. $u_i$ sont les déplacements nodaux, $F_i$ sont les forces nodales (appliquées aux noeuds).*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "1. Calculer la matrice de rigidité gloable:\n", " 1. en utilisant les équations d'équilibre pour chaque noeud \n", " 2. en assemblant les matrices de rigidités locales des ressorts." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "\n", " **Place your answer here** \n", "\n", " ---\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " \n", "2. On applique les conditions limites $u_1 = u_3 = 0$ et $u_4 = \\delta > 0$, et on pose $F_2 = 0$ (aucune force appliquée sur 2). Quelle est la valeur de $u_2$ ? Quelle est la force de réaction au noeud 4 ($F_4$) ?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "\n", " **Place your answer here** \n", "\n", " ---\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "3. Sur le système initial, quel déplacement faut-il bloquer pour que les degrés de libertés restants soient indépendants les uns des autres ? Quelle est alors la particularité de la matrice de rigidité réduite ?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "\n", " **Place your answer here** \n", "\n", " ---\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exercice 2: Approximation d'une solution analytique (à la main + Python)\n", "\n", "On considère un système de quatre ressorts en série. La rigidité totale du système est $k$. Il découle donc que $k_l = 4k$.\n", "\n", "![](figs/ex2.svg)\n", "\n", "1. Donner, en fonction de $k_l$, la matrice de rigidité d'un élément." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "\n", " **Place your answer here** \n", "\n", " ---\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "2. Assembler la matrice de rigidité globale, et donner son expression en fonction de $k$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "\n", " **Place your answer here** \n", "\n", " ---\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "On applique les conditions limites $u_1 = u_5 = 0$. On applique les forces:\n", "$$F_2 = \\frac{f}{4},\\quad F_3 = \\frac{f}{4},\\quad F_4 = \\frac{f}{4}$$\n", "\n", "\n", "3. Calculer les déplacements $u_2$, $u_3$ et $u_4$ (en fonction de $f$ et de $k$)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "\n", " **Place your answer here** \n", "\n", " ---\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "4.Comparer ces valeurs avec le champ de déplacement analytique d'une barre bi-encastrée de longueur~$1$, de rigidité axiale $EA/L$ = $k$, soumise à une force axiale répartie: $q(x) = f$. Pour rappel, l'équation différentielle est:\n", " $$k\\frac{\\text{d}^2 u}{\\text{d}x^2} = -q(x),\\quad \\text{avec}\\ u(0) = 0\\ \\text{et}\\ u(1) = 0$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "\n", " **Place your answer here** \n", "\n", " ---\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "5. **Avec Python**, nous allons refaire les questions 1 à 4 pour un système à $n$ ressorts ($n$ sera une variable). \n", "\n", " **Attention**: la rigidité de chaque ressort dépend de n afin de conserver la rigidité totale\n", "k.\n", "\n", " - Puisque chaque ressort a la même matrice de rigidité, créer une variable la contenant." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Place your answer here" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " - Écrire la boucle d'assemblage de la matrice de rigidité globale (boucle sur les éléments). " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Place your answer here" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " - Écrire la boucle qui remplit le vecteur des forces ($F_i = \\frac{f}{n}$) (boucle sur les noeuds). " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Place your answer here" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " - Appliquer les conditions limites puis résoudre le système $[K]\\{u\\} = \\{F\\}$. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Place your answer here" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " - Afficher la solution obtenue. Vous pouvez résoudre un système linéaire avec le fonction `solve` du Numpy (https://numpy.org/doc/stable/reference/generated/numpy.linalg.solve.html)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Place your answer here" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- Afficher la solution analytique sous forme de graphique. Vous pourrez utiliser la fonction `plot` du module `matplotlib.pyplot` (https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.plot.html)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Place your answer here" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "6. Faire varier le nombre de ressorts pour montrer que la solution éléments finis se rapproche de la solution analytique quand le nombre de ressorts est grand." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Place your answer here" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exercice 3: Système de ressorts (à la main)\n", "\n", "Pour le système de ressorts représenté ci dessous:\n", "\n", "![](figs/ex3.svg)\n", "\n", "1. Pour chaque ressort, calculez l'énergie élastique du système ($E_\\text{el}$) en fonction des déplacements.\n", "\n", " Rappel de l'énergie potentielle d'un ressort:\n", " $$E = \\frac{1}{2}k\\Delta l^2$$\n", " Où $\\Delta l$ est la variation de longeur du ressort." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "\n", " **Place your answer here** \n", "\n", " ---\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "2. Calculez le travail des forces extérieures $F_1$ et $F_3$ ($E_\\text{Fext}$)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "\n", " **Place your answer here** \n", "\n", " ---\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "3. Calculez l'énergie potentielle totale du système $\\Pi = E_\\text{el} - E_\\text{Fext}$, en fonction de $u_1$, $u_2$ et $u_3$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "\n", " **Place your answer here** \n", "\n", " ---\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "4. Donnez les équations que doivent satisfaire $u_1$, $u_2$ et $u_3$ pour minimiser l'énergie potentielle $\\Pi$.\n", " Donner l'expression de la matrice de rigidité globale $K$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "\n", " **Place your answer here** \n", "\n", " ---\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "5. Retrouver $K$ par la méthode de rigidité directe." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "\n", " **Place your answer here** \n", "\n", " ---\n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "votre_opinion_compte('Exercice-Ressorts')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "hide_input": false, "kernelspec": { - "display_name": "Python 3", + "display_name": "Python", "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.6.9" + "version": "3.8.10" }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": true, "sideBar": true, "skip_h1_title": false, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": false, "toc_position": { "height": "calc(100% - 180px)", "left": "10px", "top": "150px", "width": "384px" }, "toc_section_display": true, "toc_window_display": true } }, "nbformat": 4, "nbformat_minor": 4 } diff --git a/examples/ExercisesA-Structures-Solution.ipynb b/examples/ExercisesA-Structures-Solution.ipynb index 6130b3c..92585e9 100644 --- a/examples/ExercisesA-Structures-Solution.ipynb +++ b/examples/ExercisesA-Structures-Solution.ipynb @@ -1,855 +1,855 @@ { "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", " Example extracted from the MNSS repository by Guillaume Anciaux, more on: https://gitlab.epfl.ch/anciaux/mnss-notebooks.\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", " CIVIL-321 \"Modélisation Numérique des Solides et Structures\"

\n", " Comment utiliser ce Jupyter Notebook?\n", "

\n", " Ce Notebook est constitué de cellules de texte et de cellule de code. Les cellules de codes doivent être executées pour voir le résultat du programme. Certaines cellules doivent être remplies par vos soins. Pour exécuter une cellule, cliquez dessus simplement et ensuite cliquez sur le bouton \"play\" (►) dans la barre de menu au dessus du notebook. Vous pouvez aussi taper la combinaison de touches shift + enter. Il est important d'éxécuter les cellules de code en respectant leur ordre d'arrivée dans le notebook.\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", " Vous avez une question ?

\n", " SpeakUp: https://speakup.epfl.ch/ng/room/05096\n", "
" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from lib.plot import *" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Série d'exercices : Systèmes de ressorts" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Introduction\n", "\n", "Cette séance d'exercice contient des exercices à faire à la main ainsi que par programmation sur Python. Tout peut être fait dans ce *Notebook*." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## À lire avant de commencer\n", "\n", "\n", "- Les index sur Python commence à 0 (et non à 1 **contrairement à Matlab**). \n", "- Les matrices sont crées via le module *Numpy* (https://numpy.org/doc/stable/reference/arrays.ndarray.html)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# exemple, creation d'une matrice \n", "import numpy as np\n", "A = np.array([[1,2],[3, 4]])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- La fonction `plot_matrix(matrix, matrix_name)` permet de visualiser les matrices. \n", " Elle prend en entrée: \n", " - matrix : le nom de la variable \n", " - matrix_name = 'nom_matrix'" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# exemple d'affichage\n", "plot_matrix(A, 'A')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- l'extraction de bloc d'un vecteur ou d'une matrice se fait via les `slice` (https://numpy.org/doc/stable/reference/arrays.indexing.html)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# exraction d'une ligne\n", "plot_matrix(A[0, :], 'A_{{0, x}}')\n", "\n", "# exraction d'une colonne\n", "plot_matrix(A[:, 1], 'A_{{x, 1}}')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exercice 1: Assemblage de matrice (à la main)\n", "\n", "![](figs/ex1.svg)\n", "\n", "*Système de ressorts considéré. $u_i$ sont les déplacements nodaux, $F_i$ sont les forces nodales (appliquées aux noeuds).*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "1. Calculer la matrice de rigidité gloable:\n", " 1. en utilisant les équations d'équilibre pour chaque noeud \n", " 2. en assemblant les matrices de rigidités locales des ressorts." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "\n", " **Solution:**\n", "\n", " ---\n", "\n", "\n", "\n", "- Isolons le noeud 2:\n", " ![](figs/ex1-1.svg)\n", " \n", " *Figure: Équilibre du noeud 2. La force appliquée est $F_2$. Les forces $f_1$, $f_2$ et $f_3$ sont les forces de rappel des ressorts 1, 2, et 3.*\n", "\n", " L'équilibre donne comme équation: $F_2 - f_1 - f_2 - f_3 = 0$\n", " \n", " On remplace $f_1$, $f_2$, $f_3 \\qquad \\Longrightarrow \\qquad F_2 -k_1(u_2 - u_1) + k_2(u_3 - u_2) + k_3(u_4-u_2) = 0$\n", " \n", " En répétant le processus pour chaque noeud, nous trouvons un système d'équations linéaires dont la matrice est la matrice de rigidité globale:\n", " $$ \\mathbf{K} = \\begin{bmatrix}\n", " k_1 & -k_1 & 0 & 0\\\\\n", " -k_1 & k_1 + k_2 + k_3 & -k_2 & -k_3\\\\\n", " 0 & -k_2 & k_2 & 0\\\\\n", " 0 & -k_3 & 0 & k_3\n", " \\end{bmatrix}$$\n", " \n", "- La matrice de rigidité globale est l'assemblage des matrices de rigidité de chaque ressort:\n", " $$\\mathbf{K} = \\sum_{e=1}^N {\\mathbf{k}^e}$$\n", " \n", " Pour chaque ressort, la matrice de rigidité locale est:\n", " $$\\mathbf{k}^e = \\begin{bmatrix}\n", " k_e & -k_e\\\\\n", " -k_e & k_e\n", " \\end{bmatrix}$$\n", "\n", " Il faut cependant être attentif et se demander pour chaque noeud quels ressorts contribuent à la rigidité. \n", " Par exemple:\n", " - Pour le noeud 1, seul le ressort 1 contribue. La première ligne et la première colonne de $\\mathbf{K}$ vont donc être déterminées par la matrice de rigidité du ressort 1.\n", " - Le noeud 2 est relié à tous les ressorts. La deuxième ligne et la deuxième colonne de $\\mathbf{K}$ seront donc déterminées par les trois matrices de rigidité.\n", "\n", " En applicant ce raisonnement pour chaque noeud, on obtient:\n", " \n", " $$\\mathbf{K} = \\begin{bmatrix}\n", " k_1 & -k_1 & 0 & 0\\\\\n", " -k_1 & k_1 + k_2 + k_3 & -k_2 & -k_3\\\\\n", " 0 & -k_2 & k_2 & 0\\\\\n", " 0 & -k_3 & 0 & k_3\n", " \\end{bmatrix}$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "2. On applique les conditions limites $u_1 = u_3 = 0$ et $u_4 = \\delta > 0$, et on pose $F_2 = 0$ (aucune force appliquée sur 2). Quelle est la valeur de $u_2$ ? Quelle est la force de réaction au noeud 4 ($F_4$) ?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "\n", " **Solution:**\n", "\n", " ---\n", "\n", "\n", "\n", "On applique les conditions limite suivantes:\n", " - $u_1 = 0$\n", " - $u_3 = 0$\n", " - $u_4 = \\delta$\n", "\n", "Les deux premières conditions sont des conditions **homogènes**. On peut donc supprimer les lignes et colones de $\\mathbf{K}$ correspondant à $u_1$ et $u_3$. Le système $\\mathbf{K}u = F$ devient:\n", "\n", "$$\\begin{pmatrix}\n", "k_1 + k_2 + k_3 & -k_3\\\\\n", "-k_3 & k_3\n", "\\end{pmatrix} \\begin{pmatrix}\n", "u_2\\\\\n", "u_4\n", "\\end{pmatrix} = \\begin{pmatrix}\n", "0\\\\\n", "F_4\n", "\\end{pmatrix}$$\n", "\n", "Ce système a deux équations et deux inconnues: $u_2$ et $F_4$, puisque $u_4 = \\delta$. La résolution donne:\n", "$$\\begin{eqnarray*}\n", "u_2 & = & \\frac{k_3 \\delta}{k_1 + k_2 + k_3}\\\\\n", "F_4 & = & k_3 \\delta \\left( 1 - \\frac{k_3}{k_1 + k_2 + k_3} \\right)\n", "\\end{eqnarray*}$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "3. Sur le système initial, quel déplacement faut-il bloquer pour que les degrés de libertés restants soient indépendants les uns des autres ? Quelle est alors la particularité de la matrice de rigidité réduite ?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "\n", " **Solution:**\n", "\n", " ---\n", "\n", "\n", "\n", "- Si le noeud 2 est bloqué, des forces ou des déplacements appliqués aux noeuds restants n'auront pas d'effet sur les autres noeuds. Si l'on applique la condition limite homogène $u_2 = 0$, la matrice de rigidité devient:\n", " $$\\mathbf{K} = \\begin{bmatrix}\n", " k_1 & 0 & 0\\\\\n", " 0 & k_2 & 0\\\\\n", " 0 & 0 & k_3\n", " \\end{bmatrix}$$\n", " \n", " $\\mathbf{K}$ est diagonale, ce qui exprime bien le fait qu'il n'y ai plus de relation entre les degrés de liberté du système." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exercice 2: Approximation d'une solution analytique (à la main + Python)\n", "\n", "On considère un système de quatre ressorts en série. La rigidité totale du système est $k$. Il découle donc que $k_l = 4k$.\n", "\n", "![](figs/ex2.svg)\n", "\n", "1. Donner, en fonction de $k_l$, la matrice de rigidité d'un élément." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "\n", " **Solution:**\n", "\n", " ---\n", "\n", "\n", "\n", "$$\\mathbf{k}^e = \\begin{bmatrix}\n", " k_l & -k_l\\\\\n", " -k_l & k_l\n", " \\end{bmatrix}$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "2. Assembler la matrice de rigidité globale, et donner son expression en fonction de $k$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "\n", " **Solution:**\n", "\n", " ---\n", "\n", "\n", "\n", "- $$\\mathbf{K} = k \\begin{pmatrix}\n", " 4 & -4 & 0 & 0 & 0\\\\\n", " -4 & 8 & -4 & 0 & 0\\\\\n", " 0 & -4 & 8 & -4 & 0\\\\\n", " 0 & 0 & -4 & 8 & -4\\\\\n", " 0 & 0 & 0 & -4 & 4\\\\\n", " \\end{pmatrix}$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "On applique les conditions limites $u_1 = u_5 = 0$. On applique les forces:\n", "$$F_2 = \\frac{f}{4},\\quad F_3 = \\frac{f}{4},\\quad F_4 = \\frac{f}{4}$$\n", "\n", "\n", "3. Calculer les déplacements $u_2$, $u_3$ et $u_4$ (en fonction de $f$ et de $k$)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "\n", " **Solution:**\n", "\n", " ---\n", "\n", "\n", "\n", "$$k\\begin{pmatrix}\n", " 8 & -4 & 0\\\\\n", " -4 & 8 & -4\\\\\n", " 0 & -4 & 8\n", " \\end{pmatrix} \\begin{pmatrix}\n", " u_2\\\\\n", " u_3\\\\\n", " u_4\n", " \\end{pmatrix} = \\frac{f}{4}\\begin{pmatrix}\n", " 1\\\\\n", " 1\\\\\n", " 1\n", " \\end{pmatrix}$$\n", " Ce qui donne:\n", " $$\\begin{pmatrix}\n", " u_2\\\\\n", " u_3\\\\\n", " u_4\n", " \\end{pmatrix} = \\frac{f}{4k}\n", " \\begin{pmatrix}\n", " \\frac{3}{8}\\\\[0.1cm]\n", " \\frac{1}{2}\\\\[0.1cm]\n", " \\frac{3}{8}\n", " \\end{pmatrix}$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "4.Comparer ces valeurs avec le champ de déplacement analytique d'une barre bi-encastrée de longueur~$1$, de rigidité axiale $EA/L$ = $k$, soumise à une force axiale répartie: $q(x) = f$. Pour rappel, l'équation différentielle est:\n", " $$k\\frac{\\text{d}^2 u}{\\text{d}x^2} = -q(x),\\quad \\text{avec}\\ u(0) = 0\\ \\text{et}\\ u(1) = 0$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "\n", " **Solution:**\n", "\n", " ---\n", "\n", "\n", "\n", "- L'équation différentielle se résoud par intégration successive, pour\n", " donner: $u(x) = \\frac{f}{2k}(x-x^2)$\n", " \n", " Les déplacements aux positions de $u_2$, $u_3$ et $u_4$ sont:\n", " $$\\begin{eqnarray*}\n", " u(\\frac{1}{4}) & = & \\frac{3f}{32k} = u_2\\\\\n", " u(\\frac{1}{2}) & = & \\frac{f}{8k} = u_3\\\\\n", " u(\\frac{3}{4}) & = & \\frac{3f}{32k} = u_4\n", " \\end{eqnarray*}$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "5. **Avec Python**, nous allons refaire les questions 1 à 4 pour un système à $n$ ressorts ($n$ sera une variable). \n", "\n", " **Attention**: la rigidité de chaque ressort dépend de n afin de conserver la rigidité totale\n", "k.\n", "\n", " - Puisque chaque ressort a la même matrice de rigidité, créer une variable la contenant." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "###########\n", "# Solution:\n", "##########\n", "\n", "\n", "n = 4\n", "k = 1\n", "kl = n*k\n", "kloc = kl*np.array([[1, -1], [-1, 1]])\n", "plot_matrix(kloc, 'k^{{loc}}')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- Écrire la boucle d'assemblage de la matrice de rigidité globale (boucle sur les éléments)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "###########\n", "# Solution:\n", "##########\n", "\n", "\n", "# On se sert de la connectivité pour assembler les matrices locales\n", "# dans la matrice globale\n", "# La taille de K est (N+1, N+1) car il y a N+1 noeuds\n", "\n", "K = np.zeros((n+1, n+1))\n", "\n", "for e in range(n):\n", " for i in range(2):\n", " for j in range(2):\n", " K[e+i, e+j] = K[e+i, e+j] + kloc[i, j]\n", " # ce qui s'écrit aussi\n", " # K[e+i, e+j] += kloc[i, j]\n", " \n", "plot_matrix(K, 'K')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- Écrire la boucle qui remplit le vecteur des forces ($F_i = \\frac{f}{n}$) (boucle sur les noeuds)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "###########\n", "# Solution:\n", "##########\n", "\n", "\n", "# On remplit le vecteur force, de taille N+1\n", "# les valeurs de F(1) et F(N+1) ne sont pas importantes,\n", "# elles sont ignorées lors de la résolution de Ku=F\n", "\n", "f = 1\n", "F = np.ones(n+1)*f/n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- Appliquer les conditions limites puis résoudre le système $[K]\\{u\\} = \\{F\\}$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "###########\n", "# Solution:\n", "##########\n", "\n", "\n", "# On bloque les noeuds 1 et N+1\n", "# Les noeuds libres sont donc les noeuds de 2 à N\n", "\n", "noeuds_libres = range(1, n)\n", "K_libre = K[noeuds_libres, :][:, noeuds_libres]\n", "F_libre = F[noeuds_libres]\n", "\n", "plot_matrix(K_libre, 'K^{{libre}}')\n", "plot_matrix(F_libre, 'F^{{libre}}')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- Afficher la solution obtenue. Vous pouvez résoudre un système linéaire avec le fonction `solve` du Numpy (https://numpy.org/doc/stable/reference/generated/numpy.linalg.solve.html)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "###########\n", "# Solution:\n", "##########\n", "\n", "\n", "# On construit le vecteur des déplacements complet\n", "u = np.zeros(n+1)\n", "\n", "# Résolution de Ku = F\n", "u[noeuds_libres] = np.linalg.solve(K_libre, F_libre)\n", "plot_matrix(u, 'u')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "- Afficher la solution analytique sous forme de graphique. Vous pourrez utiliser la fonction `plot` du module `matplotlib.pyplot` (https://matplotlib.org/stable/api/_as_gen/matplotlib.pyplot.plot.html)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "###########\n", "# Solution:\n", "##########\n", "\n", "\n", "import matplotlib.pyplot as plt\n", "\n", "x_numeric = np.linspace(0, 1, n+1)\n", "x_analytic = np.linspace(0, 1, 1000)\n", "u_analytic = f / (2 * k) * (x_analytic - x_analytic**2);\n", "plt.plot(x_numeric, u, label='numeric')\n", "plt.plot(x_analytic, u_analytic, '--', label='analytic')\n", "plt.legend()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "6. Faire varier le nombre de ressorts pour montrer que la solution éléments finis se rapproche de la solution analytique quand le nombre de ressorts est grand." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "###########\n", "# Solution:\n", "##########\n", "\n", "\n", "# On écrit une fonction qui résoud notre problème pour un n donné\n", "def solve_problem(n):\n", " kl = n*k\n", " kloc = kl*np.array([[1, -1], [-1, 1]])\n", "\n", " K = np.zeros((n+1, n+1))\n", "\n", " for e in range(n):\n", " for i in range(2):\n", " for j in range(2):\n", " K[e+i, e+j] += kloc[i, j]\n", " \n", " f = 1\n", " F = np.ones(n+1)*f/n\n", " \n", " noeuds_libres = range(1, n)\n", " K_libre = K[noeuds_libres, :][:, noeuds_libres]\n", " F_libre = F[noeuds_libres]\n", "\n", " u = np.zeros(n+1)\n", " u[noeuds_libres] = np.linalg.solve(K_libre, F_libre)\n", " x = np.linspace(0, 1, n+1)\n", " return x, u\n", "\n", "# on fait une boucle pour faire varier n\n", "for n in range(2, 10, 2):\n", " x, u = solve_problem(n)\n", " plt.plot(x, u, label=f'analytic $n = {n}$')\n", " \n", "plt.plot(x_analytic, u_analytic, '--', label='analytic')\n", "plt.legend()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exercice 3: Système de ressorts (à la main)\n", "\n", "Pour le système de ressorts représenté ci dessous:\n", "\n", "![](figs/ex3.svg)\n", "\n", "1. Pour chaque ressort, calculez l'énergie élastique du système ($E_\\text{el}$) en fonction des déplacements.\n", "\n", " Rappel de l'énergie potentielle d'un ressort:\n", " $$E = \\frac{1}{2}k\\Delta l^2$$\n", " Où $\\Delta l$ est la variation de longeur du ressort." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "\n", " **Solution:**\n", "\n", " ---\n", "\n", "\n", " \n", "$$\\begin{equation*}\n", " E_\\text{el} = \\frac{1}{2} k_1 {\\delta_1}^2 + \\frac{1}{2} k_2 {\\delta_2}^2 +\n", " \\frac{1}{2} k_3 {\\delta_3}^2 + \\frac{1}{2} k_4 {\\delta_4}^2\n", " \\end{equation*}$$\n", "\n", "où $\\delta_i$ est l'allongement du $i^\\text{ème}$ ressort. Exprimons ces allongement en fonction des déplacements:\n", "\n", "$$\\begin{align*}\n", " \\delta_{1} &= u_1 - u_2 \\\\\n", " \\delta_{2} &= u_2 \\\\\n", " \\delta_{3} &= u_3 - u_2 \\\\\n", " \\delta_{4} &= -u_3.\n", " \\end{align*}$$\n", "\n", "L'énergie potentielle totale peut donc être réécrite en fonction des déplacements:\n", "\n", "$$\\begin{equation*}\n", " E_\\text{el} = \\frac{1}{2} k_1 (u_1 - u_2)^2 + \\frac{1}{2} k_2 {u_2}^2 +\n", " \\frac{1}{2} k_3 (u_3 - u_2)^2 + \\frac{1}{2} k_4 {u_3}^2\n", " \\end{equation*}$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "2. Calculez le travail des forces extérieures $F_1$ et $F_3$ ($E_\\text{Fext}$)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "\n", " **Solution:**\n", "\n", " ---\n", "\n", "\n", "\n", "Le travail des forces extérieures est:\n", " $$E_\\text{Fext} = F_1 u_1 + F_3 u_3$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "3. Calculez l'énergie potentielle totale du système $\\Pi = E_\\text{el} - E_\\text{Fext}$, en fonction de $u_1$, $u_2$ et $u_3$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "\n", " **Solution:**\n", "\n", " ---\n", "\n", "\n", "\n", "L'énergie potentielle totale du système est:\n", " $$\\Pi(u_1, u_2, u_3) = \\frac{1}{2} k_1 (u_1 - u_2)^2 + \\frac{1}{2} k_2 {u_2}^2 +\n", " \\frac{1}{2} k_3 (u_3 - u_2)^2 + \\frac{1}{2} k_4 {u_3}^2 - F_1\n", " u_1 - F_3 u_3$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "4. Donnez les équations que doivent satisfaire $u_1$, $u_2$ et $u_3$ pour minimiser l'énergie potentielle $\\Pi$.\n", " Donner l'expression de la matrice de rigidité globale $K$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "\n", " **Solution:**\n", "\n", " ---\n", "\n", "\n", " \n", "Minimiser $\\Pi$ par rapport à $u_1$,$u_2$ et $u_3$ nécessite\n", "\n", "$$\\begin{equation*}\n", " \\frac{\\partial \\Pi}{\\partial u_i} = 0 \\quad\\quad i=1,2,3.\n", " \\end{equation*}$$\n", "\n", "Les équations en résultant sont:\n", "\n", "$$\\begin{align*}\n", " \\frac{\\partial \\Pi}{\\partial u_1} &= k_1 \\left(u_1 - u_2\\right) - F_1 = 0 \\\\\n", " \\frac{\\partial \\Pi}{\\partial u_2} &= -k_1 \\left(u_1 - u_2\\right) + k_2 u_2 - k_3 \\left(u_3 - u_2\\right) = 0\\\\\n", " \\frac{\\partial \\Pi}{\\partial u_3} &= k_3 \\left(u_3 - u_2\\right) + k_4 u_3 - F_3 = 0.\n", " \\end{align*}$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "5. Retrouver $K$ par la méthode de rigidité directe." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "\n", " **Solution:**\n", "\n", " ---\n", "\n", "\n", "\n", "Elles peuvent être réécrites sous la forme matricielle\n", "\n", "$$\\begin{equation*}\n", " \\left[\\begin{array}{ccc}\n", " k_1 & -k_1 & 0 \\\\\n", " -k_1 & k_1+k_2+k_3 & -k_3 \\\\\n", " 0 & -k_3 & k_3+k_4\n", " \\end{array}\\right] \\left\\{\\begin{array}{c}\n", " u_1 \\\\ u_2 \\\\ u_3\n", " \\end{array}\\right\\} = \\left\\{\\begin{array}{c}\n", " F_1 \\\\ 0 \\\\ F_3\n", " \\end{array}\\right\\}\n", " \\end{equation*}$$" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "votre_opinion_compte('Exercice-Ressorts')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "hide_input": false, "kernelspec": { - "display_name": "Python 3", + "display_name": "Python", "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.6.9" + "version": "3.8.10" }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": true, "sideBar": true, "skip_h1_title": false, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": false, "toc_position": { "height": "calc(100% - 180px)", "left": "10px", "top": "150px", "width": "384px" }, "toc_section_display": true, "toc_window_display": true } }, "nbformat": 4, "nbformat_minor": 4 } diff --git a/examples/TextbookA-SignalProcessing.ipynb b/examples/TextbookA-SignalProcessing.ipynb index f2063c8..3dc237f 100644 --- a/examples/TextbookA-SignalProcessing.ipynb +++ b/examples/TextbookA-SignalProcessing.ipynb @@ -1,396 +1,396 @@ { "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", " Example extracted from the Signal Processing repository (COM303) by Paolo Prandoni, more on: https://github.com/prandoni/COM303.\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Hearing the phase of a sound\n", "\n", "In this notebook we will investigate the effect of phase on the perceptual quality of a sound. It is often said that the human ear is largely insensitive to phase and that's why most of the equalization in commercial-grade audio equipment takes place in the magnitude domain only.\n", "\n", "But is it really so? Let's find out." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%matplotlib inline\n", "import matplotlib\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "import IPython\n", "from scipy.io import wavfile" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.rcParams[\"figure.figsize\"] = (14,4)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We will be synthesizing audio clips so let's set the sampling rate for the rest of the notebook:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "Fs = 16000 # sampling freqency\n", "TWOPI = 2 * np.pi" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We will be synthesizing and playing audio clips so let's define a convenience function to \"beautify\" the resulting sound: basically, we want a gentle fade-in and fade-out to avoid abrupt \"clicks\" when the waveform begins and ends. \n", "\n", "Also, there is a \"bug\" in the current version of IPython whereby audio data is forcibly normalized prior to playing (see [here](https://github.com/ipython/ipython/issues/8608) for details; this may have been solved in the meantime). On the other hand, we want to avoid normalization in order to keep control over the volume of the sound. A way to do so is to make sure that all audio clips have at least one sample at a pre-defined maximum value, and this value is the same for all clips. To do so we add a slow \"tail\" to the data which will not result in an audible sound but will set a common maximum value in all clips." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def prepare(x, max_value = 3):\n", " N = len(x)\n", " # fade-in and fade-out times max 0.2 seconds\n", " tf = min(int(0.2 * Fs), int(0.1 * N))\n", " for n in range(0, int(tf)):\n", " s = float(n) / float(tf)\n", " x[n] *= s\n", " x[N-n-1] *= s\n", " # let's append an anti-normalization tail; drawback is one second of silence in the end\n", " x = np.concatenate((x, np.linspace(0, max_value, int(Fs/2)), np.linspace(max_value, 0, int(Fs/2))))\n", " return x" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 1) Sustained sounds\n", "\n", "The first experiment will use sustained sounds, i.e. waveform whose global envelope does not change in time. A periodic sustained waveform is simply the sum of harmonically-related sinusoidal components, i.e. a set of sines or cosines with different amplitudes and phases at frequencies multiple of a fundamental. The fundamental frequency is also known as the pitch of the sound.\n", "\n", "A class of instruments that produce good sustained sounds is the woodwinds (clarinet, saxophone, etc). The mechanics of sound generation in woodwinds are beyond the scope of this notebook and plenty of references can be found on the internet. For our purposes we will choose to simulate a clarinet according to the analysis described [here](http://www.phy.mtu.edu/~suits/clarinet.html).\n", "\n", "![title](figs/clarinet.png)\n", "\n", "A clarinet-like sustained sound will contain frequencies at odd multiples of the fundamental. We will just keep the fundamental and five harmonics and we be able to specify the initial phase offset for each component:" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "outputs": [], "source": [ "def clarinet(f, phase = []):\n", " # length in seconds of audio clips\n", " T = 3\n", " \n", " # we will keep 5 harmonics and the fundamental\n", " # amplitude of components: \n", " ha = [0.75, 0.5, 0.14, 0.5, 0.12, 0.17]\n", " \n", " # phase\n", " phase = np.concatenate((phase, np.zeros(len(ha)-len(phase))))\n", "\n", " x = np.zeros((T * Fs))\n", " # clarinet has only odd harmonics\n", " n = np.arange(len(x))\n", " for k, h in enumerate(ha):\n", " x += h * np.sin(phase[k] + TWOPI * (2*k + 1) * (float(f)/Fs) * n)\n", " return x" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# fundamental frequency: D4\n", "D4 = 293.665\n", "x = clarinet(D4)\n", "\n", "# let's look at the waveform, nice odd-harmonics shape: \n", "plt.plot(x[0:300])\n", "plt.show()\n", "\n", "# and of course we can play it (using our preparing function):\n", "IPython.display.Audio(prepare(x), rate=Fs)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Ok, so it's not the best clarinet sound in the universe but it's not bad for just a few lines of code. Now let's see how changing the phase affects the sound. Let's just use random phase offsets for the components: we can see that the waveform doesn't look too nice anymore:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "xrp = clarinet(D4, [3.84, 0.90, 3.98, 4.50, 4.80, 2.96])\n", "\n", "plt.plot(xrp[0:300])\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# but if we play it, it sounds the same! \n", "IPython.display.Audio(prepare(xrp), rate=Fs)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "OK, so it seems that phase is not important after all. To check once again, run the following notebook cell as many times as you want and see if you can tell the difference between the original zero-phase and a random-phase sustained note (the phases will be different every time you run the cell):" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "xrp = clarinet(D4, np.random.rand(6) * TWOPI)\n", "plt.plot(xrp[0:300])\n", "plt.show() \n", "IPython.display.display(IPython.display.Audio(prepare(x), rate=Fs))\n", "IPython.display.display(IPython.display.Audio(prepare(xrp), rate=Fs))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 1) Dynamic sounds\n", "\n", "In the second experiment we will use real-world dynamic sounds, i.e. sounds that display time-varying characteristics. Typically, a physical musical instrument will produce sounds whose envelope displays four subsequent portions:\n", "\n", "* the **attack** time is the time taken for the sound to go from silence to max amplitude\n", "* the **decay** time is the time taken for the sound to decrease to sustain level\n", "* the **sustain** time is the time during the sound is kept at the same amplitude\n", "* the **release** time is the time taken for sound to go to zero after the stimulation is stopped.\n", "\n", "![title](figs/piano.jpg)\n", "\n", "Consider for instance a piano note: the attack time is very quick (the hammer hits the string); the decay is quite rapid as the string settles into harmonic equilibrium but there is no sustain since once the hammer hits, the stimulation ends. So a piano note has a distinct volume envelope that rises very fast and then releases slowly:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from scipy.io import wavfile\n", "Fs, x = wavfile.read(\"resources/piano.wav\")\n", "plt.plot(x)\n", "plt.show() \n", "IPython.display.Audio(x, rate=Fs)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "By now we know that the \"shape\" of a waveform is largely encoded in the phase. It is no surprise, therefore, that if we mess up with the phase of the piano sample above, we will get something that looks very different.\n", "\n", "To see this, let's take the DFT of the audio data, set the phase to zero and take the IDFT:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# first some prep work; let's make sure that\n", "# the length of the signal is even \n", "# (it will be useful later)\n", "if len(x) % 2 != 0:\n", " x = x[:-1]\n", "\n", "# let's also store the maximum value for our \n", "# \"prepare\" function \n", "mv = int(max(abs(x)) * 1.2)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Let's take the Fourier transform\n", "X = np.fft.fft(x)\n", "\n", "# we can plot the DFT and verify we have a nice \n", "# harmonic spectrum\n", "plt.plot(np.abs(X[0:int(len(X)/2)]))\n", "plt.show() " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# now we set the phase to zero; we just need to\n", "# take the magnitude of the DFT\n", "xzp = np.fft.ifft(np.abs(X))\n", "\n", "# in theory, xzp should be real; however, because\n", "# of numerical imprecision, we're left with some imaginary crumbs:\n", "print (max(np.imag(xzp)) / max(np.abs(xzp)))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# the imaginary part is negligible, as expected, \n", "# so let's just get rid of it\n", "xzp = np.real(xzp)\n", "\n", "# and now we can plot:\n", "plt.plot(xzp)\n", "plt.show() " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Gee, what happened?!? Well, by removing the phase, we have destroyed the timing information that, for instance, made the sharp attack possible (mathematically, note that by creating a zero-phase spectrum we did obtain a symmetric signal in the time domain!).\n", "\n", "If we play the waveform, we can hear that the pitch and some of the timbral quality have been preserved (after all, the magnitude spectrum is the same), but the typical piano-like envelope has been lost." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "IPython.display.Audio(prepare(xzp, mv), rate=Fs)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can amuse ourselves with even more brutal phase mangling: let's for instance set a random phase for each DFT component. The only tricky thing here is that we need to preserve the Hermitian symmetry of the DFT in order to have a real-valued time-domain signal:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# we know the signal is even-length so we need to build\n", "# a phase vector of the form [0 p1 p2 ... pM -pM ... -p2 -p1]\n", "# where M = len(x)/2\n", "ph = np.random.rand(int(len(x) / 2) ) * TWOPI * 1j\n", "# tricky but cute Python slicing syntax...\n", "ph = np.concatenate(([0], ph, -ph[-2::-1]))\n", "\n", "# now let's add the phase offset and take the IDFT\n", "xrp = np.fft.ifft(X * np.exp(ph))\n", "\n", "# always verify that the imaginary part is only roundoff error\n", "print (max(np.imag(xrp))/max(np.abs(xrp)))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "xrp = np.real(xrp)\n", "plt.plot(xrp)\n", "plt.show()\n", "\n", "IPython.display.Audio(prepare(xrp, mv), rate=Fs)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Pretty bad, eh? So, in conclusion, phase is very important to the temporal aspects of the sound, but not so important for sustained sounds. In fact, the brain processes the temporal and spectral cues of sound very differently: when we concentrate on attacks and sound envelope, the brain uses time-domain processing, whereas for pitch and timbre, it uses primarily the magnitude of the spectrum!" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { - "display_name": "Python 3", + "display_name": "Python", "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.6.9" + "version": "3.8.10" } }, "nbformat": 4, "nbformat_minor": 4 }