2017-10-16 3 views
1

Je tente de coder un algorithme particulier, appelé l'algorithme de Verlet, qui met à jour les valeurs de certaines variables, appelées x, v et E, chaque fois que l'algorithme s'exécute.Pour les boucles avec plusieurs variables - Python

import matplotlib.pyplot as plt 
import numpy as np 


# set constants 
k = 1 
m = 1 

#Set the number of iterations to be used in the for loop 
no_of_iterations=1001 



# make arrays to define the data 
t = np.zeros(no_of_iterations) 
x = np.zeros(no_of_iterations) 
v = np.zeros(no_of_iterations) 
E = np.zeros(no_of_iterations) 

# set initial conditions 
t[0] = 0 
x[0] = 0 
v[0] = 1 
E[0] = 0 


# make arrays to store data 
x1=[] 
v1=[] 
E1=[] 




'''N = 4 Loops''' 
#loop for dt = 0.1,N = 4,j=1 
for N in range(1,5,1): 
    for i in range(1,no_of_iterations): 
      j = 1 
      x[0]= np.sin((np.pi*j*k)/(N+1)) 
      t_max = 100.0 
      dt = t_max/no_of_iterations #time step 
      t[i] = dt * i 
      v[i] = v[i-1] - ((dt *k/m*(x[N]-x[N-1]))/2) 
      x[i] = x[i-1] + dt * v[i] 
      E[i] = ((((v[i])**2)/2) + (x[N]-x[N-1])) 
      x1.append(x[i]) 
      v1.append(v[i]) 
      E1.append(E[i]) 

Essentiellement ce que cela fait est annexant les valeurs mises à jour de x, v et E pour vider les listes, alors que dans certaines conditions initiales.

Maintenant, mon problème vient du fait que mes variables que je bouclez ont beaucoup de valeurs différentes:

N = 4,16,128: dt = 0.1,0.01: j = 1, N/2

et je suis censé tracer chaque matrice créer, et l'exemple de ce qui est indiqué ci-dessous:

Terrain pour N = 4, dt = 0,1, j = 1

enter image description here

(Avec mon N = 4 boucles)

Un problème est que je stocke toutes les valeurs de N = 1,2,3,4 dans une liste, appelée x1, alors que je veux réellement avoir 4 listes séparées, une pour chaque valeur de N que je parcoure, ie. x1_N1, x2_N2, x3_N3, x4_N4. Mais même si je le fais, les choses deviendront plutôt désordonnées plutôt rapidement. Par exemple, avec N = 128, je devrais avoir 4 graphiques ((dt = 0.1 et j = 1), (dt = 0.01 et j = 1), (dt = 0.01 et j = N/2), (dt = 0.01 et j = N/2)) chacun avec 128 listes. Donc je devrais créer 8 boucles pour juste N = 4!

Alors, existe-t-il un moyen de créer une nouvelle liste dans cette boucle for et de l'ajouter, plutôt que de devoir définir une liste vide en dehors de celle-ci?

Edit: Merci les gars de commenter!

J'utilise l'algorithme Verlet pour résoudre l'équation de mouvement de plusieurs oscillateurs couplés:

Equation of motion in the form of a Hamiltonian

H est la variable E dans mon code. N est le nombre d'oscillateurs. Dv/Dx = x_n - x_N-1 est la force.

L'algorithme réel est fait par ce qui suit:

Verlet Algoithm Implementation

où τ est dt dans mon code.

Maintenant, ce que je suis en train de faire avec ce code est:

What I'm doing with this code

Pour être honnête, je ne suis pas tout à fait sûr de ce que j est! Je pense que c'est juste une variable utilisée pour changer la condition initiale de la variable de position x?

Encore une fois merci pour votre aide!

+0

Pouvez-vous expliquer ce que signifie chaque variable? 'N' est le nombre de séries à tracer,' dt' est l'incrément de temps, mais qu'est-ce que 'j'? – Adirio

+0

Je pense aussi que vous avez un problème d'algorithme lorsque vous calculez 'v [i]' car vous utilisez 'x [N]' pour le calculer. – Adirio

+0

S'il vous plaît faire une section purement sur la physique du problème. Quelle est la fonction d'énergie, la force, l'accélération, pourquoi est-ce conservateur. Il semble que vous souhaitiez implémenter la méthode de Leapfrog Verlet, mais il semble que ce soit mal fait. Avec une initialisation erronée, vous obtenez l'ordre 1 méthode Euler symplectique. Quelles sont exactement les différentes courses, que comparez-vous? – LutzL

Répondre

0

J'ai laissé les instructions d'impression commentées mais elles ne sont pas utiles dans l'exemple fourni car il y a beaucoup d'itérations et de combinaisons et le terminal débordera. Si vous modifiez la liste pour avoir une seule combinaison possible et réduisez un peu le nombre d'itérations, les décommentera affichera chaque étape. Ne pas oublier également de décommenter np.set_printoptions() alors. J'ai aussi laissé les légendes de matplotlib commentées car 129 légendes sont aussi trop nombreuses pour s'afficher correctement. Je suis en train de tracer tous les x s avec E et tous les v avec E.

j=0 est le cas particulier où toutes les conditions initiales sont 0 sauf x_(N/2) = 1. Le reste de j agira comme le deuxième cas de condition initiale que vous avez écrit.

import matplotlib.pyplot as plt 
import numpy as np 


#np.set_printoptions(precision=3, suppress=True) 


combinations = [] 
for no_iterations in [1000]: 
    for t_end in [10, 100]: 
     for N in [4, 16, 128]: 
      for j in [0, 1, N//2]: 
       combinations.append({ 
        "no_iterations": no_iterations, 
        "t_end": t_end, 
        "N": N, 
        "j": j, 
       }) 


for vars_dict in combinations: 
    ''' Extract the vars from the dict ''' 
    no_iterations = vars_dict["no_iterations"] 
    t_end = vars_dict["t_end"] 
    dt = t_end/no_iterations 
    N = vars_dict["N"] 
    j = vars_dict["j"] 
    ''' Print the variables ''' 
    print(" N = {:3} | dt = {:.2f} | j = {:2}".format(N, dt, j)) 
    #print("-"*80) 
    ''' Create the arrays for the different values ''' 
    v = np.zeros((no_iterations + 1, N)) 
    x = np.zeros((no_iterations + 1, N)) 
    E = np.zeros(no_iterations + 1) 
    ''' Initial conditions (t=0) ''' 
    if j == 0: 
     x[0][N//2 - 1] = 1 
    else: 
     for k in range(1, N+1): 
      x[0][k-1] = np.sin(np.pi*j*k/(N+1)) 
    E[0] = 0.5 * (sum(map(lambda x: x**2, v[0])) + sum(map(lambda x, y: x - y, x[0][1:], x[0][:-1]))) 
    #print("t={:4} x={} v={} E={:+.4f}".format(0, x[0], v[0], E[0])) 
    ''' Loop (t=1..no_iterations) ''' 
    for t in range(1, no_iterations+1): 
     for k in range(1, N+1): 
      if k == 1: 
       u = x[t-1][0] - x[t-1][1] 
      elif k == N: 
       u = x[t-1][N-1] - x[t-1][N-2] 
      else: 
       u = 2* x[t-1][k-1] - x[t-1][k-2] - x[t-1][k] 
      v[t][k-1] = v[t-1][k-1] - dt * u/2 
      x[t][k-1] = x[t-1][k-1] + dt * v[t][k-1] 
     E[t] = 0.5 * (sum(map(lambda x: x**2, v[t])) + sum(map(lambda x, y: x - y, x[t][1:], x[t][:-1]))) 
     #print("t={:4} x={} v={} E={:+.4f}".format(t, x[t], v[t], E[t])) 
    ''' The x axis for plotting ''' 
    t = np.linspace(0, t_end, no_iterations+1) 
    fig, (ax1, ax2) = plt.subplots(2, 1) 
    for k in range(1, N+1): 
     ax1.plot(t, x[:,k-1], label="x_{}".format(k)) 
     ax2.plot(t, v[:,k-1], label="v_{}".format(k)) 
     pass 
    ax1.plot(t, E, label="Energy") 
    ax2.plot(t, E, label="Energy") 
    # ax1.legend() 
    # ax2.legend() 
    plt.show() 
    print("-"*80) 
    print() 

d'abord la création d'un im list appelé combinations qui a un dict s avec toutes les combinaisons possibles variables extraites des 4 listes dans les quatre premières boucles for. J'ai été tenté d'utiliser des constantes au début du script pour montrer qu'elles peuvent être éditées par l'utilisateur mais la dernière list dépend de N et cela ne peut être connu qu'à l'intérieur de la troisième boucle donc je les ai toutes laissées dans les boucles.

Ensuite, la boucle suivante parcourt toutes les combinaisons et extrait les variables du dict pour une syntaxe plus propre. Je calcule également dt.

Je crée deux où le premier matrices indice représente l'instant de temps et le second sous-indice en cas de x_n et v_n et un vecteur pour E. Les premiers instants de ce matrices/vecteurs sont initiés comme vous l'avez décrit, en utilisant le cas particulier de j == 0 pour le cas particulier que vous avez mentionné.

Une boucle est alors entré à parcourir tous les instants restants et tous les x_n(t), v_n(t) sont calculées d'abord pour que nous puissions calculer E(t). C'est la partie qui peut être erronée car j'ai fait l'hypothèse que u est calculée à partir de l'instant précédent. Si vous savez comment résoudre ce genre de problèmes à la main, vous pourrez peut-être corriger cette partie au cas où cette hypothèse serait fausse.

Nous créons alors le vecteur t pour l'axe x à la fin (comme nous utilisions la variable t dans la boucle précédente pour une autre chose) avec np.linspace(start, end, number) qui donne un vecteur avec number éléments fron start à end, y compris les deux.

Nous utilisons un 2 rangées complotent pour imprimer x_n et E dans la partie supérieure et v_n et E à nouveau dans la partie inférieure.

+0

Edité pour inclure les nouvelles informations données dans le PO – Adirio

+0

Merci beaucoup pour votre aide! Il y a beaucoup de code vraiment cool là-dedans! J'ai passé un certain temps à le déchiffrer et j'ai tellement appris! J'apprécie sincèrement toute votre aide. –

+0

Vous auriez pu demander des explications et j'aurais pu vous aider là-bas. Éditera pour l'inclure. – Adirio