Modes propres d'oscillations des systèmes discrets

Information générales du micro-contenu.[1]

FondamentalDiaporama du micro-contenu (rafraichir la page si redimensionnée)

Lien vers les sources H5P du diaporama [zip]

Texte légalPolycopié du micro-contenu

Sources LaTeX du polycopié [zip]

SimulationProgramme Python

1
# -*- coding: utf-8 -*-
2
"""
3
Discrete Laplacian
4
O. Thual, 21/08/2021
5
"""
6
7
#  clear all
8
for iglob in list(globals().keys()):
9
    if(iglob[0] != '_'):
10
        exec('del {}'.format(iglob))
11
# import libraries
12
import numpy as np
13
import matplotlib.pyplot as plt
14
import os
15
16
17
def inifig(xaxe=0,yaxe=0,xlab='x',ylab='y'):
18
    plt.figure(2)
19
    plt.axvline(xaxe)
20
    plt.axhline(yaxe)
21
    plt.xticks(fontsize=12)
22
    plt.yticks(fontsize=12)
23
    plt.xlabel(xlab,fontsize=16 )
24
    plt.ylabel(ylab,fontsize=16)
25
    
26
def zfi(x,le=2):
27
    miss=le-len(str(x))
28
    a='0'*miss+str(x)
29
    return a
30
        
31
def anim(Vecp,m,N):
32
    a=np.linspace(0,1,N+1)
33
    acont=np.linspace(0,1,100)
34
    sc=.05; pi=np.pi; 
35
    n=1+m; print("n=",n)
36
    xim=np.concatenate(([0],Vecp[:,m],[0]))/np.max(Vecp)
37
    cor=xim[1]/np.sin(n*pi*a[1])
38
    Nt=12
39
    for i in range(0,Nt):
40
        title="N="+zfi(N)+" n="+zfi(n)+" i="+zfi(i)
41
        titlefig="N="+zfi(N,1)+" n="+zfi(n,1)
42
        print(title)
43
        fig=plt.figure(N,figsize=(7,4))
44
        plt.xlabel('a',fontsize=16 )
45
        plt.ylabel(r'Déplacement $\xi$',fontsize=16)
46
        plt.title(titlefig,fontsize=16)
47
        la=np.sin(2*pi*i/Nt)
48
        xi=la*xim
49
        xicont=la*np.sin(n*pi*acont)*cor
50
        plt.scatter(a,xi,marker='o',color='blue',s=40)
51
        plt.scatter(a+sc*xi,0*xi,marker='o',color='red',s=40)
52
        plt.scatter(a,0*xi,marker='o',color='green',s=40)
53
        plt.plot(acont,xicont,color='black')
54
        plt.plot(a,0*xi,color='black')
55
        plt.ylim(-1,1)
56
        basN='ResN='+zfi(N)
57
        basNn=basN+'n='+zfi(n)
58
        basNni=basNn+'i='+zfi(i)+'.png'; 
59
        #box=([0,0],[3,3]); plt.savefig(basNni,bbox_inches=Bbox(box))
60
        plt.savefig(basNni)
61
        #plt.show()
62
        plt.close()
63
        gifanim="Anim"+basNn+'.gif'
64
    os.system('/opt/local/bin/convert -set delay 20 '+basNn+'* '+gifanim);
65
    os.system('rm '+basNn+'* ');   
66
    return
67
68
def disclapl(Nmax,figvp):  
69
    ########## All #############
70
    for N in range(2,Nmax+1):
71
        print("======== N=",N,"==========")
72
        # Define the matrix AN
73
        AN=np.zeros((N-1,N-1))
74
        AN[0,0]=-2; 
75
        if N > 2:
76
            AN[0,1]=1; AN[N-2,N-3]=1; AN[N-2,N-2]=-2;
77
            for i in range(1,N-2):
78
                AN[i,i-1]=1; AN[i,i]=-2; AN[i,i+1]=1;
79
        #print("AN="); print(AN)
80
        # Compute eigenvalues and eigenmodes
81
        Valp, Vecp=np.linalg.eig(AN)     
82
        Omega=N*(-Valp)**(.5)/np.pi
83
        #sort_perm = Valp.argsort(); Valp.sort(); Vecp=Vecp[:, sort_perm]
84
        sort_perm = Omega.argsort(); 
85
        Valp=Valp[sort_perm] 
86
        Vecp=Vecp[:, sort_perm]
87
        #print("Valp=",Valp)
88
        #print("N/pi sqrt(Valp)=",Omega)
89
        #print("Vecp",Vecp)
90
        if figvp:
91
            Num=np.linspace(N,N,N-1)
92
            ymin=0; ymax=10
93
            plt.yticks(np.linspace(0,10,11))
94
            plt.ylim(ymin,ymax) 
95
            plt.grid(color='black', axis='y', linestyle='-', linewidth=1)        
96
            plt.scatter(Num,Omega,marker='o',color='red',s=20)
97
            plt.xlabel(r'$N$',fontsize=16 )
98
            plt.ylabel(r'Pulsations $\Omega_n$',fontsize=16)
99
            plt.title(r'Pulsations propres normalisées',fontsize=16)
100
101
            namepng='Omega-'+str(N)+'.png'; plt.savefig(namepng)
102
        else:    
103
            AnigResN='AnigResN='+zfi(N)
104
            AnimResN='AnimResN='+zfi(N)
105
            for m in range(N-1):
106
                #print("m=",m," Valp[m]=",Valp[m], " Vecp[:,m]=",Vecp[:,m])
107
                anim(Vecp,m,N)
108
            os.system('/opt/local/bin/convert '+AnimResN+'*.gif '+AnigResN+'.gif');
109
110
            
111
    # Figure des valeurs propres   
112
    if figvp:     
113
        plt.savefig('Omega.pdf')
114
        plt.show(); plt.close()
115
        for n in range(N,N+10):
116
            os.system('cp '+namepng+' Omega-'+str(n)+'.png')
117
        os.system("/opt/local/bin/convert Omega-* NimOmega.gif");
118
        os.system('rm Omega-*');
119
    else:
120
        os.system('/opt/local/bin/convert AnigResN=*.gif AnigResALL.gif');
121
122
123
# Main 
124
F=False; T=True
125
if F: Nmax=20; figvp=True; disclapl(Nmax,figvp)  
126
if F: Nmax=5; figvp=False; disclapl(Nmax,figvp)
127
if T: Nmax=10; figvp=False; disclapl(Nmax,figvp)
128
# N=4 
129
if F:
130
    s1=-2+np.sqrt(2); s2=-2; s3=-2-np.sqrt(2)
131
    Omega1=4/np.pi*np.sqrt(-s1); print('Omega1=',Omega1)
132
    Omega2=4/np.pi*np.sqrt(-s2); print('Omega2=',Omega2)
133
    Omega3=4/np.pi*np.sqrt(-s3); print('Omega3=',Omega3)
134
    
135
136