複数の変数を持つループの場合 - Python

Question

Verlet Algorithmと呼ばれる特定のアルゴリズムをコーディングしようとしています。アルゴリズムを実行するたびに、x、v、Eと呼ばれる特定の変数の値が基本的に更新されます。

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]) 

基本的には、特定の初期条件が与えられている間に、x、v、Eの更新値を空のリストに追加することです。 、J = 1：

N = 4,16,128：DT = 0.1,0.01

は今私の問題は、私はオーバーループするよ私の変数は、多くの異なる値を持っているという事実から来ていますN/2

及びIは、Iが作成した各配列をプロットすることになっていて、以下に示す例は：N = 4

プロット、DT = 0.1、J = 1

（私のN付= 4ループ）

一つの問題は、私が実際にしたい一方で、私は、X1と呼ばれる、1つのリストにN = 1,2,3,4のすべての値を格納していますということです4つの別々のリストを持ちます.1つはNの各値に対して1つ、つまり繰り返します。 x1_N1、x2_N2、x3_N3、x4_N4。

でも、私がそれをしてもかなり素早くかなり面倒になるでしょう。例えば、N = 128の場合、グラフはdt = 0.01、j = 1、dt = 0.01、j = 1、dt = 0.01、j = N/2、 j = N/2））に128個のリストがあります。したがって、N = 4のために8つのループを作成する必要があります。

だから、このforループ内で新しいリストを作成して追加する方法はありますか？その外に空のリストを定義する必要はありませんか？

編集：コメントありがとうございました！

私は多くの結合振動子の運動方程式を解くためにVerletアルゴリズムを使用しています：

Hは、私のコード内の変数Eです。 Nは発振器の数です。 Dv/Dx = x_n - x_N-1が力である。

実際のアルゴリズムは、以下により行われます

τは私のコードではdtのです

。

今、私は実際にこのコードでやっているのです：

正直に言うと、私は何であるかjを全くわからないんだけど！私はそれが位置変数xの初期条件を変更するために使用される変数だと思いますか？

もう一度ご協力いただきありがとうございます。

Answer 1

私はprintステートメントをコメントしましたが、多くの反復と組み合わせがあり、ターミナルがオーバーフローするため、提供された例では役に立ちません。リストを変更して1つだけの組み合わせを作成し、反復回数を少し減らすと、それらのコメントを外して各ステップを表示します。また、np.set_printoptions()のコメントを外すことを忘れないでください。

また、129凡例が正しく表示されるには多すぎるため、matplotlibの凡例にはコメントが残っています。 Eとvのすべてのxの両方をプロットしています。Eです。

j=0は、x_(N/2) = 1を除き、すべての初期条件が0である特殊なケースです。 jの残りは、あなたが書いた2番目の初期条件ケースとして機能します。 4つの第1 forループで4つのリストから抽出された全ての可能な変数の組み合わせとdict Sを有するlist呼ばcombinationsを作成

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() 

まずイム。私はスクリプトの始めに定数を使ってユーザーが編集できることを示したが、最後のlistはNに依存しており、3番目のループの中でしか知ることができないので、すべてループに残した。

次に、すべての組み合わせに対して次のループを繰り返し、より洗練された構文のためにdictから変数を抽出します。私もdtと計算しています。

x_nとv_nの場合のサブインデックスとEの場合は、最初のインデックスが時刻を表し、2番目のサブインデックスが2つの行列を作成しています。この行列/ベクトルの最初の瞬間はあなたが記述した特別な場合のj == 0という特別な場合を使って、記述したように開始されます。

ループを入力して残りのすべてのインスタントを繰り返し、最初にx_n(t),v_n(t)が計算され、E(t)と計算されます。これは、私が前の瞬間からuが計算されたことを間違えた部分です。この種の問題を手で解決する方法が分かっている場合、仮定が間違っている場合にこの部分を修正することができます。

（私たちは別のものの前のループ内変数tを使用していたように）我々は、両方を含むendにstartフロンnumber要素を持つベクトルを与えるnp.linspace(start, end, number)と終了時にX軸のtベクトルを作成しています。

私たちは、底部に再びx_nとE上部にあるとv_nとEを印刷するには2行のプロットを使用しています。