#################################################################### ##### CAPITOLO I ##### #################################################################### # Prog_1_1.py # Calcola la quota al variare del # tempo di un grave in caduta libera. # RICHIESTA DATI da console quota_in = float(input('quota iniziale in m: ')) # dt incremento della variabile t dt = float(input('dt in secondi: ')) # INIZIALIZZAZIONE di g e variabili g = 9.81 # m/sec^2 t = 0 # sec quota = quota_in # m # CICLO di CALCOLO e STAMPA while quota >= 0: quota = quota_in - (1 / 2) * g * (t**2) if quota >= 0: print('tempo: {:.2f} s quota: {:.2f} m'.format(t,quota)) t = t + dt #################################################################### ##### CAPITOLO II ##### #################################################################### # Prog_2_1.py # Traccia la curva oraria dalla tabella oraria. import numpy as np import matplotlib.pyplot as plt # crea array t=0,2,4,...28,30 t = np.arange(0,32,2) # Lista con i valori di s s = [0.0,255.9,622.3,845.4,1199.5,1418.0,1788.1,1950.3, 2000.0,1875.8,1662.6,1401.5,961.7,665.4,273.2,0.0] plt.plot(t, s) # genera il file grafico plt.xlabel('tempo [min]') # aggiunge descrizione asse x plt.ylabel('posizione [m]') # aggiunge descrizione asse y plt.grid() # crea reticolo nel piano cartesiano plt.show() # visualizza grafico #################################################################### # Prog_2_2.py # Traccia le curve orarie dei moti uniformemente accelerati # generici A e B. Le condizioni richieste sono s0_A < s0_B # per le posizioni iniziali e a_A > a_B per le accelerazioni. # Il codice si arresta quando il punto A raggiunge il punto B. import matplotlib.pyplot as plt import sys # INPUT equazioni orarie s0_A = float(input('s0_A in m : ')) s0_B = float(input('s0_B in m : ')) v0_A = float(input('v0_A in m/s : ')) v0_B = float(input('v0_B in m/s : ')) a_A = float(input('a_A in m/s^2 : ')) a_B = float(input('a_B in m/s^2 : ')) # TEST pos. iniziali e accel. dei due moti A e B if (s0_A >= s0_B) or (a_A <= a_B): sys.exit('s0_A >= s0_B o a_A <= a_B') dt = 0.01 # passo di calcolo per t t = 0.0 # istante iniziale t_valori = [t] # lista valori t s_A = [s0_A] # lista posizioni moto A s_B = [s0_B] # lista posizioni moto B sA = s0_A sB = s0_B # CICLO di CALCOLO tempo-posizione while sA < sB: t = t + dt sA = 0.5 * a_A * t**2 + v0_A * t + s0_A s_A.append(sA) sB = 0.5 * a_B * t**2 + v0_B * t + s0_B s_B.append(sB) t_valori.append(t) # GRAFICI plt.plot(t_valori,s_A,label='moto A') plt.plot(t_valori,s_B,label='moto B') plt.legend(loc='best') plt.xlabel('tempo [sec]') plt.ylabel('posizione [m]') plt.grid() plt.show() #################################################################### # Prog_2_3.py # moto verticale e parabolico import numpy as np import matplotlib.pyplot as plt v0 = 8.25 # velocità  di lancio in m/s alzo = 20 # angolo di lancio in gradi h0 = 10 # altezza iniziale in m tv = 1.45 # tempo di volo in sec x0 = 10.3 # posizione del grave moto verticale g = 9.81 # in m/sec^2 theta = alzo * (np.pi / 180) # trasformazione in radianti vx = v0 * np.cos(theta) # componente orizzontale velocità  vy = v0 * np.sin(theta) # componente verticale velocità  t = np.arange(0,tv, 0.05) xt = vx * t # spostamento orizzontale yt = h0 + vy * t - 0.5 * g * t**2 # spostamento verticale x0 = 0 * t + x0 # ascissa grave moto verticale plt.plot(xt,yt,'g:o',label='grave lanciato') plt.xlabel('x [m]',fontsize=20) plt.ylabel('y [m]',fontsize=20) plt.plot(x0,yt,'vr',label='moto verticale') plt.legend(loc=0) plt.grid() plt.xticks(fontsize=15) plt.yticks(fontsize=15) plt.show() #################################################################### # Prog_2_4.py # moto circolare e moto armonico import numpy as np import matplotlib.pyplot as plt r = 1 # raggio della circonferenza omega = 1 # velocità  angolare T = 6.4 # durata del moto t = np.arange(0, T, 0.3) # array degli istanti di tempo xt = r * np.cos(omega * t) # calcolo delle x di P yt = r * np.sin(omega * t) # calcolo delle y di P plt.axis('equal') # stessa unità  di misura assi plt.plot(xt,yt,'g:o',label='moto circolare uniforme(x,y)') plt.plot(t,yt,'r:o',label='moto armonico(t,x)') plt.legend(loc='upper right') plt.grid() plt.xticks(fontsize=16) plt.yticks(fontsize=16) plt.show() #################################################################### # Prog_2_5.py # moto armonico: grafici spostamento, velocità  e accelerazione import numpy as np import matplotlib.pyplot as plt r = 1.0 # ampiezza omega = 1.0 # velocità angolare T = 6.5 # durata del moto dt = 0.025 # passo di calcolo per t NUM_passi = int(T / dt) # numero iterazioni t = [0.0] * NUM_passi # lista per i valori di t x = [0.0] * NUM_passi # lista per i valori posizione v = [0.0] * NUM_passi # lista per i valori velocità  a = [0.0] * NUM_passi # lista per i valori accelerazione i = 1 while i < NUM_passi: t[i] = (i-1) * omega * dt # incremento calcolo x[i] = r * np.cos(t[i]) # calcolo posizione v[i] = (x[i] - x[i-1]) / dt # calcolo velocità  a[i] = (v[i] - v[i-1]) / dt # calcolo accelerazione i = i + 1 # grafica plt.plot(t[1:NUM_passi],x[1:NUM_passi],'g',label='x') plt.plot(t[2:NUM_passi],v[2:NUM_passi],'r',label='v') plt.plot(t[3:NUM_passi],a[3:NUM_passi],'b',label='a') plt.ylabel('a - v - x', fontsize=24) plt.xlabel('t',fontsize=24) plt.legend(loc='best') plt.grid() plt.xticks(fontsize=16) plt.yticks(fontsize=16) plt.show() #################################################################### # Prog_2_6.py # traiettoria elicoidale in 3D import numpy as np import matplotlib.pyplot as plt r = 1.0 # raggio orbita T = 0.12 # periodo rotazione omega = 2 * np.pi / T # velocità  angolare tempo_totale = 0.28 v = 1.0 # velocità  pianeta t=np.linspace(0,tempo_totale,36) # vettore valori parametro t z = v * t # Z pianeta # grafica figura = plt.figure() ax = figura.add_subplot(projection='3d') # traiettoria satellite ax.plot(r * np.cos(omega * t), r * np.sin(omega * t),z ,'b:.') # traiettoria pianeta ax.plot(0 * t, 0 * t, z, 'y.') plt.xlabel('X') plt.ylabel('Y') plt.show() #################################################################### ##### CAPITOLO III ##### #################################################################### # Prog_3_1.py # # Calcolo della posizione y e della velocità  v nel # caso di una forza costante (caduta di un grave). # import matplotlib.pyplot as plt g = 9.81 # accelerazione di gravità  in m/sec^2 dt = 0.01 # passo di calcolo t = 0.0 y = float(input('y0 in m : ')) # quota iniziale vy = float(input('vy in m/s : ')) # velocità  iniziale t_valori = [t] # lista valori t y_A = [y] # lista valori quota v_B = [vy] # lista valori velocità  # Ciclo di calcolo while y > 0: y = y + vy * dt # aggiorna posizione y_A.append(y) vy = vy - g * dt # aggiorna velocità  v_B.append(vy) t = t + dt t_valori.append(t) plt.subplot(2, 1, 1) # grafico quota plt.plot(t_valori, y_A) plt.ylabel('y') plt.grid() plt.subplot(2, 1, 2) # grafico velocità  plt.plot(t_valori, v_B) plt.ylabel('vel [m / sec]') plt.xlabel('t [sec]') plt.grid() plt.show() #################################################################### # Prog_3_2.py # # Calcolo della posizione x e della velocità  v nel # caso di forza elastica (oscillatore armonico). # import matplotlib.pyplot as plt x_m = 1.3 # spostamento iniziale [m] x = [x_m] t = [0] k = 1.7 # costante elastica [N/m] m = 1 # massa [kg] T = 7 # tempo totale [sec] dt = 0.0001 # passo di calcolo [sec] NUM_passi = int(T / dt) # numero cicli di calcolo v_1 = -(k / m) * x_m * (dt / 2) # primo passo calcolo v v = [v_1] for i in range(1, NUM_passi): t.append(dt * (i - 1)) a = -(k / m) * x[i - 1] x.append(x[i - 1] + v[i - 1] * dt) v.append(v[i - 1] + a * dt) plt.grid() plt.plot(t, x, 'g',label='pos [m]') plt.plot(t, v, 'r',label='vel [m / sec]') plt.legend(loc='best') plt.xlabel('t [sec]',fontsize=20) plt.xticks(fontsize=15) plt.yticks(fontsize=15) plt.show() #################################################################### # Prog_3_3.py # # Calcola l'energia cinetica, potenziale e totale # dell'oscillatore armonico. # import matplotlib.pyplot as plt x_m = 1.3 # spostamento iniziale [m] x = [x_m] t = [0] k = 1.7 # costante elastica [N/m] m = 1 # massa [kg] T = 5 # tempo totale [sec] dt = 0.0001 # passo di calcolo [sec] E_cin = [0.0] E_pot = (1 / 2) * k * x_m**2 E_tot = [E_pot] E_pot = [E_pot] NUM_passi = int(T / dt) # numero cicli di calcolo v_1 = -(k / m) * x_m * (dt / 2) # primo passo calcolo v v = [v_1] for i in range(1, NUM_passi): t.append(dt * (i - 1)) a = -(k / m) * x[i-1] x.append(x[i-1] + v[i-1] * dt) v.append(v[i-1] + a * dt) E_cin.append((1 / 2) * m * v[i-1]**2) E_pot.append((1 / 2) * k * x[i-1]**2) E_tot.append((1 / 2) * (m * v[i-1]**2 + k * x[i-1]**2)) plt.plot(t,E_cin,'g',label='E_cin') plt.plot(t,E_pot,'r',label='E_pot') plt.plot(t,E_tot,'b',label='E_tot') plt.legend(loc='best') plt.xlabel('t [sec]',fontsize=20) plt.ylabel('E [J]',fontsize=20) plt.grid() plt.xticks(fontsize=15) plt.yticks(fontsize=15) plt.show() #################################################################### # Prog_3_4.py # # Calcola e salva in un file i valori dell'angolo teta e quelli # della velocità  angolare omega di un pendolo al variare di t. # import numpy as np teta = 0.5 # angolo iniziale in rad omega = 0 # velocità  angolare iniziale [rad/s] g = 9.81 # accelerazione di gravità  [m/sec^2] L = 1 # lunghezza del pendolo [m] T = 2.05 # tempo totale [sec] dt = 0.0003 # passo di calcolo [sec] NUM_passi = int(T / dt) # numero cicli di calcolo # array contenente g e L gL = np.array([g,L]) # salva array gL nel file val_gL.npy np.save('val_gL',gL) # apre file per la scrittura dati=open('t_teta_omega.txt', 'w') for i in range(0,NUM_passi): # omega a t + dt: omega = omega - (g / L) * np.sin(teta) * dt # teta a t + dt: teta = teta + omega * dt t = dt * i dati.write('\t\t{:f}\t\t{:f}\t\t{:f}\n'\ .format(t, teta, omega)) # chiude il file t_teta_omega.txt dati.close() #################################################################### # Prog_3_5.py # # Traccia il grafico dei valori di teta e di omega # generati da Prog_3_4.py. # import matplotlib.pyplot as plt from numpy import loadtxt # carica il file con (t,teta,omega) in dati: dati = loadtxt('t_teta_omega.txt',float) # separa le colonne: t = dati[:,0] teta = dati[:,1] omega = dati[:,2] plt.figure(1) plt.plot(t,teta) plt.xlabel('t [sec]',fontsize=20) plt.ylabel('$\\theta$ [rad]',fontsize=20) plt.grid() plt.xticks(fontsize=15) plt.yticks(fontsize=15) plt.figure(0) plt.plot(teta,omega) plt.xlabel('$\\theta$ [rad]',fontsize=20) plt.ylabel('$\omega$ [rad/sec]',fontsize=20) plt.grid() plt.xticks(fontsize=15) plt.yticks(fontsize=15) plt.show() #################################################################### # Prog_3_6.py # # Calcola l'energia cinetica, potenziale e totale # del pendolo. # import matplotlib.pyplot as plt import numpy as np # massa del pendolo [kg] m = 1.0 # carica (g,L): valori = np.load('val_gL.npy') g = valori[0] L = valori[1] # carica (t,teta,omega): dati = np.loadtxt('t_teta_omega.txt',float) t = dati[:,0] teta = dati[:,1] omega = dati[:,2] E_cin = (1 / 2) * m * L**2 * omega**2 E_pot = m * g * L * (1 - np.cos(teta)) E_tot = E_cin + E_pot plt.plot(t,E_cin,'g',label='E_cin') plt.plot(t,E_pot,'r',label='E_pot') plt.plot(t,E_tot,'b',label='E_tot') plt.legend(loc='best') plt.xlabel('tempo [sec]',fontsize=20) plt.ylabel('Energia [J]',fontsize=20) plt.grid() plt.xticks(fontsize=15) plt.yticks(fontsize=15) plt.savefig('E_tot.png') #################################################################### ##### CAPITOLO IV ##### #################################################################### # Prog_4_1.py # Calcola la posizione di un pianeta e le componenti della # sua velocità  orbitale attorno al Sole e salva i dati nel # file 'i_x_y_vx_vy.txt'. import numpy as np As = 365.25 * 24 * 3600 # Anno sidereo [sec] UA = 149.59e9 # unità  astronomica [m] G = 6.673e-11 # costante G [N (m / kg)^2] M = 1.9891e30 # massa del Sole [kg] GM = G * M * (As**2 / UA**3) # valore normalizzato di G*M x0 = 1.5 # x iniziale [UA] y0 = 0 # y iniziale [UA] vx0 = 0.0 # valore iniziale v_x [UA / As] vy0 = 15.0e3 * As / UA # valore iniziale v_y [UA / As] T = 0.85 # tempo totale [As] dt = 0.00001 NUM_passi = int(T / dt) x = np.zeros(NUM_passi) y = np.zeros(NUM_passi) r = np.zeros(NUM_passi) vx = np.zeros(NUM_passi) vy = np.zeros(NUM_passi) ax = np.zeros(NUM_passi) ay = np.zeros(NUM_passi) x[0] = x0 y[0] = y0 vx[0] = vx0 vy[0] = vy0 # apre file dati=open('i_x_y_vx_vy.txt', 'w') # CICLO di CALCOLO e scrittura file for i in range(NUM_passi-1): r[i] = np.sqrt((x[i])**2 + (y[i])**2) ax[i] = - GM * ((x[i]) / (r[i]**3)) ay[i] = - GM * ((y[i]) / (r[i]**3)) vx[i+1] = vx[i] + ax[i] * dt vy[i+1] = vy[i] + ay[i] * dt x[i+1] = x[i] + vx[i+1] * dt y[i+1] = y[i] + vy[i+1] * dt # scrive file dati.write('\t\t{:}\t\t{:f}\t\t{:f}\t\t{:f}\t\t{:f}\n'\ .format(i,x[i],y[i],vx[i],vy[i])) # chiude file dati.close() #################################################################### # Prog_4_2.py # Carica il file 'i_x_y_vx_vy.txt' # e traccia l'orbita del pianeta. import matplotlib.pyplot as plt import numpy as np dati=np.loadtxt('i_x_y_vx_vy.txt',float) x=dati[:,1] y=dati[:,2] plt.xlabel('X (UA)',fontsize=20) plt.ylabel('Y (UA)',fontsize=20) plt.plot(x,y) plt.axis('equal') # visualizza posizione del Sole plt.plot(0.0, 0.0, 'yo', markersize=20) plt.grid() plt.xticks(fontsize=16) plt.yticks(fontsize=16) plt.show() #################################################################### # Prog_4_3.py # Valuta l'energia totale del sistema Sole-pianeta. import numpy as np As = 365.25 * 24 * 3600 UA = 149.59e9 G = 6.673e-11 M = 1.9891e30 C = (G * M * As**2) / (UA**3) dati=np.loadtxt('i_x_y_vx_vy.txt',float) x=dati[:,1] y=dati[:,2] vx=dati[:,3] vy=dati[:,4] r = np.sqrt(x**2 + y**2) v = np.sqrt(vx**2 + vy**2) E = 0.5 * v**2 - C / r Em=np.mean(E) vr=abs((max(E) - min(E)) / Em) print('E = {:.2f} con variazione relativa = {:.5f}'.format(Em,vr)) #################################################################### # Prog_4_4.py # Verifica I legge di Keplero (forma ellittica dell'orbita). import numpy as np dati=np.loadtxt('i_x_y_vx_vy.txt',float) x=dati[:,1] y=dati[:,2] x_afelio = x[0] x_perielio = abs(min(x)) x_Fuoco = x_afelio - x_perielio N = len(x) r_S = [0] * N r_F = [0] * N Somma_rS_rF = [0] * N for k in range(N): r_S[k] = np.sqrt((x[k])**2 + (y[k])**2) r_F[k] = np.sqrt((x[k] - x_Fuoco)**2 + (y[k])**2) Somma_rS_rF[k] = r_S[k] + r_F[k] print('k = {:.0f} r_Sole + r_Fuoco = {:.4f} UA' .format(k, Somma_rS_rF[k])) #################################################################### # Prog_4_5.py # Verifica II legge di Keplero (velocità  areolare costante) e # visualizza posizioni del pianeta a intervalli di tempo uguali. import numpy as np import matplotlib.pyplot as plt dati=np.loadtxt('i_x_y_vx_vy.txt',float) x=dati[:,1] y=dati[:,2] vx=dati[:,3] vy=dati[:,4] N = len(x) r = [0] * N v = [0] * N vxr = [0] * N vel_areolare = [0] * N for k in range(N): r[k] = [x[k],y[k],0] v[k] = [(vx[k],vy[k],0)] vxr[k] = np.cross(r[k],v[k]) vel_areolare[k] = np.linalg.norm(vxr[k]) print('velocità  areolare = {:.4f}' .format(vel_areolare[k])) # intervallo in giorni delta_giorni = 5 passo = int((N * delta_giorni) / 324) # visualizza posizioni for k in range(0,N,passo): plt.plot(x[k], y[k], 'b.') plt.grid() plt.xticks(fontsize=16) plt.yticks(fontsize=16) plt.xlabel('X (UA)',fontsize=20) plt.ylabel('Y (UA)',fontsize=20) plt.axis('equal') plt.plot(0.0, 0.0, 'yo', markersize=20) plt.show() #################################################################### # Prog_4_6.py # Calcola e traccia la traiettoria di una particella alfa # dovuta alla forza repulsiva coulombiana del nucleo atomico. import numpy as np import matplotlib.pyplot as plt Q = 79 # carica nucleo atomo oro [u.a.] k = 1 # costante di Coulomb [u.a.] q = 2 # carica particella alfa [u.a.] m = 7295 # massa particella alfa [u.a.] x0 = 0.002 # x iniziale [u.a.] y0 = 0.0002 # y iniziale [u.a.] vx0 = -7 # vx iniziale particella alfa [u.a.] vy0 = 0.0 # vy iniziale particella alfa [u.a.] T = 0.0009 # tempo di volo totale [u.a.] dt = 0.000001 NUM_passi = int(T / dt) x = np.zeros(NUM_passi) y = np.zeros(NUM_passi) vx = np.zeros(NUM_passi) vy = np.zeros(NUM_passi) v = np.zeros(NUM_passi) r = np.zeros(NUM_passi) ax = np.zeros(NUM_passi) ay = np.zeros(NUM_passi) x[0] = x0 y[0] = y0 vx[0] = vx0 vy[0] = vy0 dati=open('i_x_y_r.txt', 'w') # apre file for i in range(NUM_passi-1): r[i] = np.sqrt((x[i])**2 + (y[i])**2) ax[i] = (k * Q * q / m) * ((x[i]) / (r[i]**3)) ay[i] = (k * Q * q / m) * ((y[i]) / (r[i]**3)) vx[i+1] = vx[i] + ax[i] * dt vy[i+1] = vy[i] + ay[i] * dt x[i+1] = x[i] + vx[i+1] * dt y[i+1] = y[i] + vy[i+1] * dt dati.write('\t\t{:}\t\t{:f}\t\t{:f}\t\t{:f}\n' .format(i,x[i],y[i],r[i])) # scrittura file dati.close() # chiude file # visualizza nucleo atomo plt.plot(0.0, 0.0,'mo', markersize=15) plt.grid() plt.axis('equal') plt.plot(x,y) plt.xticks(fontsize=14) plt.yticks(fontsize=14) plt.xlabel('x ($a_0$)',fontsize=22) plt.ylabel('y ($a_0$)',fontsize=22) plt.show() #################################################################### # Prog_4_7.py # Risoluzione con odeint() dell'equazione dell'oscillatore # armonico smorzato. Traccia l'andamento della posizione y(t), # della velocità  v(t) e di (y,v). import matplotlib.pyplot as plt import numpy as np from scipy.integrate import odeint b = 0.45 # coefficiente di attrito [kg/sec] k = 1.7 # costante elastica [N/m] m = 1 # massa [kg] t0 = 0 # valore iniziale di t [sec] tf = 25 # valore finale di t [sec] y0 = 1 # spostamento iniziale [m] v0 = 0.0 # velocità  iniziale [m/sec] y0_v0 = [y0,v0] # condizioni iniziali t = np.linspace(t0,tf,500) # definizione di f1 e f2 def f1_f2(y_v,t): y, v = y_v # separazione di y e v f1 = v f2 = - (1 / m) * (k * y + b * v) return f1,f2 # integrazione equazione y_v_sol = odeint(f1_f2, y0_v0, t) y = y_v_sol[:,0] # y(t) v = y_v_sol[:,1] # v(t) plt.subplot(2,1,1) plt.plot(t,y,'g') plt.ylabel('y [m]',fontsize=20) plt.xticks(fontsize=16) plt.yticks(fontsize=16) plt.grid() plt.subplot(2,1,2) plt.plot(t,v,'r') plt.xlabel('tempo [sec]',fontsize=20) plt.ylabel('vel [m/sec]',fontsize=20) plt.xticks(fontsize=16) plt.yticks(fontsize=16) plt.grid() plt.figure(2) plt.plot(y,v) plt.xlabel('y [m]',fontsize=20) plt.ylabel('v [m/sec]',fontsize=20) plt.xticks(fontsize=16) plt.yticks(fontsize=16) plt.grid() plt.show() #################################################################### ##### CAPITOLO V ##### #################################################################### # Prog_5_1.py # Simulazione del processo di mescolamento di due tipi di # molecole (blu e magenta), dotate di moto casuale, che # all'istante iniziale sono confinate in parti opposte di # una scatola. Con NUM_molecole si indica il numero di # molecole di ciascun colore. import random import pylab as plt NUM_molecole = 10000 NUM_passi=1001 dt = 0.002 # Coordinate posizioni iniziali: -1<= Xmagenta <0, 0<= Xblu <1 # e -1<= Y <1. # Componenti spostanento: -1<= dX <1, -1<= dY <1. Xb = [random.random() for i in range(NUM_molecole)] Yb = [2 * (random.random() - 0.5) for i in range(NUM_molecole)] dXb = [2 * (random.random() - 0.5) for i in range(NUM_molecole)] dYb = [2 * (random.random() - 0.5) for i in range(NUM_molecole)] Xm= [-(random.random()) for i in range(NUM_molecole)] Ym = [2 * (random.random() - 0.5) for i in range(NUM_molecole)] dXm = [2 * (random.random() - 0.5) for i in range(NUM_molecole)] dYm = [2 * (random.random() - 0.5) for i in range(NUM_molecole)] # CICLO PRINCIPALE di calcolo for n in range(NUM_passi): # CICLO su ciascuna delle molecole for i in range(NUM_molecole): # TEST contenimento molecola_b direzione X if abs(Xb[i] + dt * dXb[i]) > 1: # inversione spostamento direzione X dXb[i] = -dXb[i] # TEST contenimento molecola_b direzione Y if abs(Yb[i] + dt * dYb[i]) > 1: # inversione spostamento direzione Y dYb[i] = -dYb[i] # CALCOLO nuova posizione Xb[i] = Xb[i] + dt * dXb[i] Yb[i] = Yb[i] + dt * dYb[i] # TEST contenimento molecola_m direzione X if abs(Xm[i] + dt * dXm[i]) > 1: # inversione spostamento direzione X dXm[i] = -dXm[i] # TEST contenimento molecola_m direzione Y if abs(Ym[i] + dt * dYm[i]) > 1: # inversione spostamento direzione Y dYm[i] = -dYm[i] # CALCOLO nuova posizione Xm[i] = Xm[i] + dt * dXm[i] Ym[i] = Ym[i] + dt * dYm[i] # SELEZIONA multipli di 200 contatore passi if n % 200 == 0: print ('Eseguiti {:.0f} passi di calcolo'.format(n)) plt.figure() plt.title('stato del sistema dopo {:.0f} passi' .format(n),fontsize=20) plt.plot(Xb,Yb,'b.',markersize=1.0) plt.plot(Xm,Ym,'m.',markersize=1.0) plt.axis((-1,1,-1,1)) plt.xticks([]) plt.yticks([]) plt.show() #################################################################### # Prog_5_2.py # Evoluzione del confinamento e della pressione di un gas # ideale inizialmente posto nella metà  di una scatola cubica. import numpy as np import random import matplotlib.pyplot as plt NUM_molecole = 20000 NUM_passi = 501 dt = 1e-7 vmod = 2 * 1e5 # modulo vettore velocità  NUM_urti = [0 for i in range(NUM_passi)] # coordinate iniziali (X,Y,Z) delle molecole # con 0 <= X < 0.5, 0 <= Y < 1, 0 <= Z < 1 X = [0.5*(random.random()) for i in range(NUM_molecole)] Y = [(random.random()) for i in range(NUM_molecole)] Z = [(random.random()) for i in range(NUM_molecole)] # componenti vettori spostamento (dX,dY,dZ) iniziali # con -0.5<=dX<0.5, -0.5<=dY<0.5, -0.5<=dZ<0.5 dX = [(random.random() - 0.5) for i in range(NUM_molecole)] dY = [(random.random() - 0.5) for i in range(NUM_molecole)] dZ = [(random.random() - 0.5) for i in range(NUM_molecole)] # array spostamento (dX,dY,dZ) dXn = np.zeros(NUM_molecole) dYn = np.zeros(NUM_molecole) dZn = np.zeros(NUM_molecole) # normalizzazione delle velocità  for i in range(NUM_molecole): dXn[i]=(vmod * dX[i]) / (np.sqrt((dX[i]**2) + (dY[i]**2) + (dZ[i]**2))) dYn[i]=(vmod * dY[i]) / (np.sqrt((dX[i]**2) + (dY[i]**2) + (dZ[i]**2))) dZn[i]=(vmod * dZ[i]) / (np.sqrt((dX[i]**2) + (dY[i]**2) + (dZ[i]**2))) # CICLO PRINCIPALE di calcolo for n in range(NUM_passi): # CICLO PRINCIPALE di calcolo # CICLO sulle molecole for i in range(NUM_molecole): if (X[i] + dt * dXn[i]) > 1: NUM_urti[n] = NUM_urti[n] + 1 # urti sulla parete X=1 # gestione confinamento molecole if (X[i] + dt * dXn[i]) > 1 or (X[i] + dt * dXn[i]) < 0: dXn[i] = -dXn[i] # inversione spostamento direzione X if (Y[i] + dt * dYn[i]) > 1 or (Y[i] + dt * dYn[i]) < 0: dYn[i] = -dYn[i] # inversione spostamento direzione Y if (Z[i] + dt * dZn[i]) > 1 or (Z[i] + dt * dZn[i]) < 0: dZn[i] = -dZn[i] # inversione spostamento direzione Z # posizione successiva X[i] = X[i] + dt * dXn[i] Y[i] = Y[i] + dt * dYn[i] Z[i] = Z[i] + dt * dZn[i] if n % 50 == 0: # seleziona multipli di 50 contatore print ("{:.0f}% dell'elaborazione".format(100 * n / NUM_passi)) # grafica figura = plt.figure() ax = figura.add_subplot(projection='3d') plt.xlim(0,1) plt.ylim(0,1) ax.plot(X,Y,Z,'b.',markersize=1.1) plt.xlabel('X') plt.ylabel('Y') fig = plt.figure() plt.xlabel('tempo',fontsize=20) plt.ylabel('numero urti',fontsize=20) u = np.arange(0,NUM_passi) plt.plot(u * dt,NUM_urti) plt.grid() plt.xticks(fontsize=15) plt.yticks(fontsize=15) plt.show() #################################################################### # Prog_5_3.py # Calcola, scelto il numero n di moli, con l'equazione # di stato dei gas perfetti pV = nRT, il valore di una # grandezza tra p, V e T, avendo inserito le altre due. import sys R = 8.3144 n = float(input('INSERIRE n (numero di moli) : ') or 0) if not n: sys.exit('NON ci sono molecole!!!') else: print('INSERIRE solo 2 valori tra p, V e T') p = float(input('pressione (in Pa) : ') or 0) V = float(input('volume (in m^3) : ') or 0) T = float(input('temperatura (in K) : ') or 0) if not p: p = (n * R * T) / V elif not V: V = (n * R * T) / p elif not T: T = (p * V) / (n * R) print('p = {:.3e} Pa' .format(p)) print('V = {:.3e} m^3' .format(V)) print('T = {:.3e} K' .format(T)) #################################################################### # Prog_5_4.py # Visualizza l'andamento della pressione P al variare della # temperatura T e del volume V secondo l'equazione di van der # Waals. Vengono tracciate anche le curve di alcune isoterme. import numpy as np import matplotlib.pyplot as plt from matplotlib import cm R = 8.3144 # costante dei gas nel S.I. a = 0.3640; b = 0.00004267 # a e b del biossido di carbonio Tmin = 260; TMax = 345 # temperature iniz. e fin. in K Vmin = 0.000065; VMax = 0.00025 # volumi iniz. e fin. in m^3 v=np.linspace(Vmin,VMax,20) # array per valori di v t=np.linspace(Tmin,TMax,20) # array per valori di t V,T=np.meshgrid(v, t) # griglia dei valori (v,t) # equazione di van der Waals P = R * T / (V - b) - (a / V**2) # crea e visualizza superficie van der Waals superficie = plt.figure() ax = superficie.add_subplot(projection='3d') ax.plot_surface(V * 1000, T, P / 1000, cmap=cm.copper) plt.xlabel('Volume [$10^{-3}m^3$]',fontsize=16) plt.ylabel('Temperatura [K]',fontsize=16) plt.show() # seleziona e traccia grafici isoterme Visot = np.linspace(Vmin,VMax,500) fig= plt.figure() deltaT = 25 # incremento temperatura delle isoterme # CICLO tracciamento isoterme for n in range(0,3): P = R * (Tmin + deltaT * n) / (Visot - b) - a / Visot**2 plt.plot(Visot * 1000, P / 1000, label=r'T=${}$K'.format(Tmin + deltaT * n)) plt.legend(loc='upper right') plt.xlabel('Volume [$10^{-3}m^3$]',fontsize=19) plt.ylabel('pressione [kPa]',fontsize=19) plt.grid(axis = 'y') # traccia linee orizzontali della griglia plt.xticks(fontsize=16) plt.yticks(fontsize=16) plt.show() #################################################################### # Prog_5_5.py # Traccia il grafico e calcola il lavoro di un ciclo costituito da # due trasformazioni isoterme (T2>T1) e da due isocore (V_B>V_A). import numpy as np import matplotlib.pyplot as plt # Costante dei gas [J/(mol*K)] R = 8.31446 # quantità  di gas [mol] n = 0.25 # Volumi isobare [m^3] e temperature isoterme [K] V_A, V_B = 0.00045, 0.0006; T1, T2 = 310, 380 # numero di valori del volume per grafici isoterme NUM_v = 200 v = np.linspace(V_A, V_B, NUM_v) p1v = (n * R * T1) / v # equazione isoterma T1 p2v = (n * R * T2) / v # equazione isoterma T2 Lav_2 = np.trapz(p2v,v) # lavoro espansione isoterma T2 Lav_1 = - np.trapz(p1v,v) # lavoro compressione isoterma T1 Lavoro_ciclo = Lav_2 + Lav_1 # isoterme T1 e T2 in kPa py1 = p1v / 1000 py2 = p2v / 1000 v_litri = v * 1000 # trasformazione da [m^3] a [litri] plt.plot(v_litri, py1, color = 'b', lw = 3, label = 'isoterma T1') plt.plot(v_litri, py2, color = 'r', lw = 3, label = 'isoterma T2') plt.fill_between(v_litri, py2, py1, facecolor = 'yellow') plt.vlines(x = V_A * 1000, ymin = py1[0], ymax = py2[0], color = 'k', lw = 3, label = 'isocora V_A') plt.vlines(x = V_B * 1000, ymin = py1[199], ymax = py2[199], color = 'm', lw = 3, label = 'isocora V_B') plt.grid() plt.xticks(fontsize=14) plt.yticks(fontsize=14) plt.xlabel('volume [litri]',fontsize=18) plt.ylabel('pressione [kPa]',fontsize=18) plt.legend(loc = 'upper right') plt.title('Lavoro compiuto = {:.2f} {}'.format (Lavoro_ciclo,'joule'),fontsize=20) plt.show() #################################################################### # Prog_5_6 # Calcola area cerchio tramite sommatoria # di rettangoli con valutazione dell'errore. import numpy as np # equazione semicirconferenza di raggio 1 def f(x): f = np.sqrt(1 - x**2) return f # numero dei rettangoli n = 100000 # base dei rettangoli delta_x = 1 / float(n) integrale = 0.0 for i in range(1,n): integrale = integrale + delta_x * f(i * delta_x) area= 4 * integrale errore = abs(area - np.pi) print('area = {:.5f}'.format(area)) print('errore = {:.5e}'.format(errore)) #################################################################### ##### CAPITOLO VI ##### #################################################################### # Prog_6_1.py # Calcola il campo elettrico generato da due cariche # Q1 e Q2 poste nelle posizioni (xQ1,yQ1) e (xQ2,yQ2) # e ne traccia le linee di forza. import numpy as np import matplotlib.pyplot as plt # griglia quadrata e coordinate N = 28 X = np.arange(0, N, 1) Y = np.arange(0, N, 1) X, Y = np.meshgrid(X, Y) # array per le componenti di E1, E2 e E E1x, E1y = np.zeros((N, N)), np.zeros((N, N)) E2x, E2y = np.zeros((N, N)), np.zeros((N, N)) Ex , Ey = np.zeros((N, N)), np.zeros((N, N)) # input valore cariche e posizioni Q1 = -1 Q2 = +4 xQ1, yQ1 = 6, 14 xQ2, yQ2 = 21, 14 k = 1 # calcolo del campo elettrico for i_y in range(N): for j_x in range(N): den1 = np.sqrt(((j_x - xQ1) ** 2 + (i_y - yQ1) ** 2)** 3) den2 = np.sqrt(((j_x - xQ2) ** 2 + (i_y - yQ2) ** 2)** 3) if (den1 != 0 and den2 != 0): # componenti campo generato da Q1 E1: E1x[i_y, j_x] = k * Q1 * (j_x - xQ1) / den1 E1y[i_y, j_x] = k * Q1 * (i_y - yQ1) / den1 # componenti campo generato da Q2 E2: E2x[i_y, j_x] = k * Q2 * (j_x - xQ2) / den2 E2y[i_y, j_x] = k * Q2 * (i_y - yQ2) / den2 # componenti campo totale E: Ex[i_y, j_x] = E1x[i_y, j_x] + E2x[i_y, j_x] Ey[i_y, j_x] = E1y[i_y, j_x] + E2y[i_y, j_x] # Plot ed effetti grafici fig = plt.figure(facecolor='k') effetto = (np.sqrt(Ex**2 + Ey**2))**(1/3) plt.plot(xQ1,yQ1, 'bo') plt.plot(xQ2,yQ2, 'ro') plt.streamplot(X, Y, Ex, Ey, color=effetto, linewidth=0.75, density=2.0, arrowstyle='->') plt.axis('equal') plt.axis('off') plt.show() #################################################################### # Prog_6_2.py # Dai valori misurati della tensione V ( in volt) e della # corrente I (in ampere) per un dato conduttore si determina # l'equazione della retta y = m * x + q che meglio esprime la # dipendenza tra y = V e x = I. Si ottiene q = 0 e posto m = R # si verifica la prima legge di Ohm V = R * I, dove R è la # resistenza (in ohm). import numpy as np import matplotlib.pyplot as plt # DATI: V = np.array([0.026, 1.0, 2.3, 3.15, 3.7, 4.5, 4.94, 6.41, 7.65, 8.25, 8.99, 9.8]) I = np.array([0.00017, 0.00066, 0.0015, 0.0021, 0.0024, 0.0030, 0.0032, 0.0042, 0.0051, 0.0055, 0.0059, 0.0065]) # m e q del fit di 1° grado dei dati (I,V) m_q = np.polyfit(I, V, 1) # m e q della retta m = m_q[0] q = m_q[1] Imax = max(I) t = np.linspace(0, Imax, 2) y = m * t + q plt.plot(I, V, 'o') plt.plot(t, y, '-',label = 'I legge di Ohm') plt.xlabel('I [ampere]',fontsize=18) plt.ylabel('V [volt]',fontsize=18) plt.grid() plt.xticks(fontsize=15) plt.yticks(fontsize=15) plt.title('R = {:.0f} {}'.format(m_q[0],'ohm'),fontsize=18) plt.show() print('q=', q) #################################################################### # Prog_6_3.py # Risoluzione con odeint() dell'equazione del circuito # costituito da un condensatore piano C, in serie con # un resistore R e un generatore di tensione continua V0. # Traccia l'andamento della carica Q(t). import matplotlib.pyplot as plt import numpy as np from scipy.integrate import odeint V_0 = 20 # ddp generatore (volt) R = 60000 # resistenza (ohm) C = 1000 * 10**-6 # capacità  (farad) t0 = 0 # valore iniziale di t (sec) tf = 400 # valore finale di t Q0 = 0 # valore iniziale di Q (coulomb) t = np.linspace(t0,tf,250) # definizione di f1 def f1(Q,t): f1 = (1 / R) * (V_0 - Q / C) return f1 # integrazione equazione circuito Q = odeint(f1, Q0, t) plt.plot(t, Q) plt.xlabel('tempo [sec]',fontsize=20) plt.ylabel('Q [C]',fontsize=20) plt.grid() plt.xticks(fontsize=15) plt.yticks(fontsize=15) plt.show() #################################################################### # Prog_6_4.py # Risoluzione con odeint() dell'equazione del circuito costituito # da un condensatore piano C, in serie con un resistore R, un # induttore L e un generatore di tensione continua V0. Traccia # l'andamento della carica Q(t) e della corrente I(t). import matplotlib.pyplot as plt import numpy as np from scipy.integrate import odeint V_0 = 10 # ddp generatore (volt) R = 25 # resistenza (ohm) C = 10**-5 # capacità  (farad) L = 0.1 # induttanza (henry) t0 = 0 # valore iniziale di t (sec) tf = 0.05 # valore finale di t (sec) Q0 = 0 # valore iniziale di Q(t) (coulomb) I0 = 0.0 # valore iniziale di I(t) (ampere) Q0_I0 = [Q0,I0] # condizioni iniziali t=np.linspace(t0,tf,500) # definizione di f1 e f2 def f1_f2(Q_I,t): Q, I = Q_I # separazione di Q e I f1 = I f2 = - (1 / L) * (R * I + Q / C - V_0) return f1,f2 # integrazione equazione circuito Q_I_sol = odeint(f1_f2, Q0_I0, t) Q = Q_I_sol[:,0] # Q(t) I = Q_I_sol[:,1] # I(t) plt.subplot(2,1,1) plt.plot(t,Q) plt.ylabel('Q [C]',fontsize=18) plt.grid() plt.xticks(fontsize=14) plt.yticks(fontsize=14) plt.subplot(2,1,2) plt.plot(t,I,'r') plt.xlabel('tempo [sec]',fontsize=14) plt.ylabel('I [A]',fontsize=18) plt.grid() plt.xticks(fontsize=14) plt.yticks(fontsize=14) plt.show() #################################################################### # Prog_6_5.py # Si calcola la velocità  di una gocciolina di olio, al variare # delle cariche elementari in essa presenti, che si muove # nell'aria in presenza di un campo elettrico e della gravità. # Otteniamo così la simulazione del fenomeno osservabile # nell'apparato di Millikan. import numpy as np import matplotlib.pyplot as plt from scipy.integrate import odeint g = 9.81 # accelerazione di gravità  [m/sec^2] r = 1.64e-6 # raggio gocciolina di olio [m] densità_olio = 851 # [kg/m^3] V = 500 # differenza di potenziale [V] d = 2.5e-3 # distanza armature [m] E = V / d # campo elettrico [N/C] e = 1.602e-19 # carica elementare [C] volume_gocciolina = 4/3 * np.pi * (r**3) NUM_e_max = 5 segno = -1 # -1 o +1 per cariche negative o positive densità_aria = 1.292 # [kg/m^3] eta = 1.82e-5 # coefficiente di viscosità  aria [kg/m*sec] # coefficiente per la forza di attrito viscoso c = 6 * np.pi * eta * r m = volume_gocciolina * densità_olio # massa gocciolina # spinta di Archimede spinta = volume_gocciolina * densità_aria * g # condizioni iniziali: y = 0, v = 0 a t = 0 y0_v0 = [0,0] t = np.linspace(0, 0.0001, 100) # definizione di f1 e f2 def f1_f2(y_v, t): y,v = y_v # separazione di y e v f1 = v f2 = g - (1 / m) * (spinta + c * v - segno * N * e * E) return f1, f2 # integrazione equazione del moto for N in range(0,NUM_e_max + 1): y_v = odeint(f1_f2, y0_v0, t) plt.plot(1000 * t, 1000 * y_v[:,1]) plt.xlabel('tempo [ms]',fontsize=18) plt.ylabel('velocità  [mm / s]',fontsize=18) plt.grid() plt.xticks(fontsize=15) plt.yticks(fontsize=15) plt.show() #################################################################### ##### CAPITOLO VII ##### #################################################################### # Prog_7_1.py # Genera immagine dello spettro di luce visibile # (intervallo lunghezza d'onda: [380;740]nm). import matplotlib.pyplot as plt # dizionario delle lunghezze d'onda e dei colori spettro = {(380,390):'#ae00ff',(390,400):'#9700ff', (400,410):'#8000ff',(410,420):'#6000ff',(420,430):'#3f00ff', (430,440):'#3200ff',(440,450):'#2600ff',(450,460):'#1900ff', (460,470):'#0d00ff',(470,480):'#0000ff',(480,490):'#003fff', (490,500):'#007fff',(500,510):'#00bfff',(510,520):'#00ffff', (520,530):'#00ffbf',(530,540):'#00ff80',(540,550):'#00ff40', (550,560):'#00ff00',(560,570):'#33ff00',(570,580):'#66ff00', (580,590):'#99ff00',(590,600):'#ccff00',(600,610):'#ffff00', (610,620):'#ffd500',(620,630):'#ffaa00',(630,640):'#ff8000', (640,650):'#ff5600',(650,660):'#ff2b00',(660,670):'#ff0000', (670,680):'#e90000',(680,690):'#d30000',(690,700):'#bd0000', (700,710):'#a70000',(710,720):'#8f0000',(720,730):'#760000', (730,740):'#610000'} # traccia bande verticali dei colori for lung_onda, colore in spettro.items(): plt.fill_between(lung_onda, 2, 1, facecolor = colore) plt.xlabel('$\lambda \;\;[10^{-9}m]$',fontsize=22) plt.xticks(fontsize=16) plt.xlim(380,740) plt.ylim(1,2) plt.yticks([]) plt.show() #################################################################### # Prog_7_2.py # Traccia due onde monocromatiche di ampiezza unitaria # e l'onda ottenuta dalla loro sovrapposizione. import matplotlib.pyplot as plt import numpy as np # lunghezze d'onda del blu e del rosso in nm Lunghezza_onda_Colore = {470:'#0000ff',665:'#ff0000'} A = 1 # ampiezza onda s = 0 x = np.arange(0, 2000, 0.5) for lunghezza_onda,colore in Lunghezza_onda_Colore.items(): k = 2 * np.pi / lunghezza_onda y = A * np.sin(k * x) plt.subplot(2,1,1) plt.plot(x, y,color = colore, label=r'${}$nm'.format(lunghezza_onda)) plt.legend(loc='best') plt.ylabel('$y_{\lambda 1},y_{\lambda 2}$',fontsize=18) plt.xticks(fontsize=14) plt.yticks(fontsize=14) plt.ylim(-2.0,2.0) s = s + y plt.grid() plt.subplot(2,1,2) plt.plot(x,s,'m') plt.ylabel('$y_{\lambda 1}+y_{\lambda 2}$',fontsize=18) plt.grid() plt.ylim(-2.0,2.0) plt.xticks(fontsize=14) plt.yticks(fontsize=14) plt.xlabel('r [nm]',fontsize=16) plt.show() #################################################################### # Prog_7_3.py # Figura di interferenza prodotta da onde circolari emesse # da due sorgenti puntiformi S1 e S2 (doppia fenditura). import numpy as np import matplotlib.pyplot as plt lunghezza_onda1 = 6 # cm lunghezza_onda2 = 6 # cm # numero d'onda in cm^-1: k1 = 2 * np.pi / lunghezza_onda1 k2 = 2 * np.pi / lunghezza_onda2 # ampiezza onda in cm: A1 = 2 A2 = 2 distanza_S1_S2 = 25 # cm lato_ondoscopio = 250 # cm # numero punti per ripartizione lato NUM_punti = 1000 delta_lato = lato_ondoscopio / NUM_punti # ampiezza intervalli # coordinate sorgenti S1 e S2 xS1 = 0 yS1 = lato_ondoscopio / 2 + (distanza_S1_S2 / 2) xS2 = 0 yS2 = lato_ondoscopio / 2 - (distanza_S1_S2 / 2) # array per contenere i valori onda1 + onda2 x_y = np.empty([NUM_punti,NUM_punti],float) # Doppio ciclo calcolo interferenza for n in range(1,NUM_punti): y = delta_lato * n for m in range(1,NUM_punti): x = delta_lato * m r1 = np.sqrt((x - xS1)**2 + (y - yS1)**2) r2 = np.sqrt((x - xS2)**2 + (y - yS2)**2) x_y[n,m] = A1 * np.sin(k1 * r1) + A2 * np.sin(k2 * r2) # GRAFICA plt.plot(xS1,yS1,'yo') plt.plot(xS2,yS2,'yo') plt.imshow(x_y, extent=[0,lato_ondoscopio,0,lato_ondoscopio]) plt.xlabel('cm') plt.ylabel('cm') plt.colorbar(label='ampiezza onda (in cm)') plt.show() #################################################################### # Prog_7_4.py # Calcola con due metodi numerici la curva del decadimento # radioattivo e si confrontano i risultati con la legge del # fenomeno. Si verifica anche la relazione quantitativa tra il # tempo di dimezzamento dei nuclei radioattivi e la loro vita # media di sopravvivenza. I dati qui usati sono quelli di un # grammo dell'isotopo radio-226. import matplotlib.pyplot as plt import numpy as np from scipy.integrate import odeint lamb = 1.36e-11 # lambda in sec^-1 Anno = 365 * 24 * 60 * 60 # secondi in un anno t_fin = 7000 * Anno # tempo finale in sec N0 = 2.66e21 # numero iniziale di nuclei tau = 1 / lamb # vita media in sec n = 100 t = np.linspace(0,t_fin,n) delta_t = t_fin / n # definizione di f1 def f1(N,t): f1 = -lamb * N return f1 # e integrazione equazione Node = odeint(f1, N0, t) # soluzione approssimata nn = np.arange(n) N = N0 * (1 - lamb * delta_t)**nn # soluzione esatta N_esatta = N0 * np.exp(- lamb * t) # verifica periodo di dimezzamento di N n_periodi = 4 for k in range(1,n_periodi+1): plt.vlines(x = k * 0.693 * tau / Anno, ymin = 0, ymax = N0 * np.exp(-k * 0.693 * lamb * tau), lw = -0.15 * k + 1.2, color='r') plt.hlines(N0 * np.exp(-k * 0.693 * lamb * tau), xmin = 0, xmax = k * 0.693 * tau / Anno, lw = -0.15 * k + 1.2, label=r'N=${}$N$_0$' .format(1 / (2**k)), color='r') # GRAFICA plt.plot(t / Anno,N,label='N (soluz. appross.)') plt.plot(t / Anno, Node,label='N (odeint)') plt.plot(t / Anno,N_esatta,label='N$_0$e$^{-\lambda t}$') plt.legend(loc='best') plt.xlabel('t [anni]',fontsize=20) plt.ylabel('N(t)',fontsize=20) plt.grid() plt.xticks(fontsize=18) plt.yticks(fontsize=18) plt.xlim(0, t_fin / Anno) plt.ylim(0, N0) plt.show() #################################################################### # Prog_7_5.py # Risoluzione equazione non algebrica con fsolve. import numpy as np from scipy.optimize import fsolve # definizione di f def f(x): f = np.exp(x) - (1 / 2) return f # risoluzione equazione x_sol = fsolve(f,0) print(x_sol) #################################################################### # Prog_7_6.py # Si descrive la diffusione di una malattia in una popolazione # di N individui ripartita in S (suscettibili), I (infetti) # e R (guariti e divenuti immuni). L'evoluzione del processo # dipende da I0 e R0 (infetti e eventuali rimossi all'istante 0) # e dai parametri beta e gamma. import matplotlib.pyplot as plt import numpy as np from scipy.integrate import odeint N = 5000 # numero totale individui I0 = 1 # numero iniziale infetti R0 = 0 # numero iniziale rimossi S0 = N - I0 - R0 # numero iniziale suscettibili beta = 0.18 # grado di contagiosità  gamma = 0.08 # tasso medio di guarigione in gg^-1 durata_media = 1 / gamma # durata decorso malattia in giorni S0_I0_R0 = S0, I0, R0 # condizioni iniziali n = 250 t_fin = 180 # tempo finale in giorni t = np.linspace(0, t_fin, n) # definizione di f1, f2 e f3 def f1_f2_f3(S_I_R, t, beta, gamma): S, I, R = S_I_R f1 = -beta * S * I / N f2 = beta * S * I / N - gamma * I f3 = gamma * I return f1, f2, f3 # integrazione equazioni SIR S_I_R_sol = odeint(f1_f2_f3, S0_I0_R0, t, args=(beta, gamma)) S = S_I_R_sol[:,0] I = S_I_R_sol[:,1] R = S_I_R_sol[:,2] # GRAFICA plt.plot(t,S,'b', lw=2, label='Suscettibili') plt.plot(t,I,'r', lw=2, label='Infetti') plt.plot(t,R,'g', lw=2, label='Rimossi') plt.xlabel('tempo [giorni]',fontsize=18) plt.grid() plt.xticks(fontsize=18) plt.yticks(fontsize=18) plt.legend(loc='best') plt.title('N={}, I0={}, $\\beta$={:.2f}, durata media malattia\ ={:.1f} giorni'.format(N,I0,beta,durata_media),fontsize=18) plt.xlim(0,t_fin) plt.show() #################################################################### # Prog_7_7.py # Le equazioni dell'attrattore di Lorenz sono risolte. # La dipendenza molto sensibile dalle condizioni iniziali # delle soluzioni è posta in evidenza con le traiettorie # tracciate nello spazio delle fasi. import matplotlib.pyplot as plt import numpy as np from scipy.integrate import odeint # condizioni iniziali u0_v0_w0 = [1.0, 1.0, 1.0] u0_v0_w0p = [1.0, 1.0, 1.00001] t_fin = 50.0 dt = 0.001 t = np.arange(0, t_fin, dt) # parametri sigma = 10 beta = 8/3 ro = 28 # definizione di f1, f2 e f3 def f1_f2_f3(u_v_w, t, sigma, beta, ro): u, v, w = u_v_w f1 = sigma * (v - u) f2 = ro * u - v - u * w f3 = u * v - beta * w return f1, f2, f3 # integrazione equazioni u_v_w_sol = odeint(f1_f2_f3, u0_v0_w0, t, args=(sigma, beta, ro)) u_v_w_sol_p = odeint(f1_f2_f3, u0_v0_w0p, t, args=(sigma, beta, ro)) # GRAFICA fig = plt.figure() ax = fig.add_subplot(projection='3d') fig.subplots_adjust(left=0, right=1, bottom=0, top=1) plt.xlabel('u',fontsize=22) plt.ylabel('v',fontsize=22) plt.plot(u_v_w_sol[:, 0],u_v_w_sol[:, 1], u_v_w_sol[:, 2],color='b',alpha=0.5) plt.plot(u_v_w_sol_p[:, 0],u_v_w_sol_p[:, 1], u_v_w_sol_p[:, 2],color='r',alpha=0.5) plt.show() #################################################################### # Prog_7_8.py # Calcolo delle funzioni d'onda e delle energie di un # elettrone confinato in un segmento di lunghezza L. import matplotlib.pyplot as plt import numpy as np from scipy.integrate import odeint # richiede E di prova da console E = float(input ('IMMETTERE VALORE di E (in joule x 10^-19) = ')) E = E * 10**-19 # joule m = 9.1e-31 # massa elettrone in kg hbar = 1.0545e-34 # joule x sec x0 = 0 # inizio confinamento in m x_L = 1.0e-9 # inizio confinamento in m psi0_Dpsi0 = [0,1] # condizioni iniziali x=np.linspace(x0,x_L,1000) # x per rappresentazione grafica # definizione di f1 e f2 def f1_f2(psi_Dpsi,t): psi, Dpsi = psi_Dpsi # separazione di psi e Dpsi f1 = Dpsi f2 = -(2 * m / hbar**2) * E * psi return f1,f2 # integrazione equazione psi_Dpsi_sol = odeint(f1_f2, psi0_Dpsi0, x) psi = psi_Dpsi_sol[:,0] # funzione d'onda # normalizzazione della funzione d'onda area = np.trapz(psi**2,x) psi = psi / np.sqrt(area) # GRAFICA plt.plot(x,psi,lw=2) plt.xlabel('$x [m]$',fontsize=20) plt.ylabel('$\psi(x)$',fontsize=26) plt.grid() plt.xticks(fontsize=18) plt.yticks(fontsize=18) plt.show() ####################################################################