For循環for循環：擺錘建模

Question

我是新來編程和工作的一個相對複雜的問題。我正在模擬一個線性強迫擺，並且需要弄清擺動運動幅度（弧度）如何取決於摩擦力（q值）和驅動力頻率（Omega_D）。所以，我認爲我需要for循環內的for循環，因爲我需要做3個q值的Omega_D與Amplitude的關係圖。因此，對於q重複3次，並且在Omega_D中重複3次。但是，我寫的只是給每個q值一個幅度值。這是我的代碼;讓我知道你可能有什麼建議。

import numpy as np 
import matplotlib.pyplot as plt 
from ode import rk4_step 


def derivs(t, starting_values): 

    d0dt = starting_values[1] 
    d20dt2 = g/l * starting_values[0] - starting_values[3] * \ 
    starting_values[1] + starting_values[4] * np.sin(starting_values[5] * t) 

    dqdt = 0. 
    dFdt = 0. 
    dldt = 0. 
    d_Omega_dt = 0.     # defining these for later use in RK4 

    derivatives = np.array([d0dt, d20dt2, dqdt, dFdt, dldt, d_Omega_dt]) 
    return derivatives 



qs = [0.1, 1.0, 1.6] #Pick arbitrary q-values to run through 
for i in qs: 

    theta_0 = 10. #initial values, chosen at random 
    theta_v0 = 10. 
    l = 1. 
    Omega_D = np.linspace(0.5, 5, 100) 
    F_D = .3 
    for x in Omega_D: 
     starting_values = [theta_0, theta_v0, l, i, F_D, x] 
     solution = np.copy(starting_values) 
     last_values = np.zeros(solution.size) 

    dt = .01 
    g = -9.8 

    t = 0. 
    Amp = 0. 


    starttime = 150. 

    while Amp == 0. : #Amp==0 because it never actually WILL be zero. 
       #So, the while loop to find amplitude only needs to run \ 
       #until a nonzero amp is found 
     two_values_ago = np.copy(last_values) 
     last_values = np.copy(solution) 


     t = t + dt 
     solution = rk4_step(solution, derivs, t, dt) #take a step 


     if solution[1] == 0 and t > starttime: #if somehow we hit v=0 exactly 
      Amp == np.abs(solution[0]) 
      print Amp 

#This if statement interpolates to find the amp at the point where v=0 
     if solution[1] * last_values[1] < 0 and t > starttime: 
      fit_vs = np.array([two_values_ago[1], last_values[1]]) 
      fit_xs = np.array([two_values_ago[0], last_values[0]]) 
      v_interp = 0. 
      Amp = np.abs(np.interp(v_interp, fit_vs, fit_xs)) 


    w = np.sqrt(-g/l) # This is the natural frequency 

#Calculate the analytic soln 
    exact_solution = F_D/np.sqrt((w**2 - Omega_D**2)**2 + (i * Omega_D)**2) 

#plot num and exact solns together 
    plt.plot(Omega_D, exact_solution) 
    plt.plot(Omega_D, Amp) 
    plt.title('q = ') 
    plt.ylabel('Amplitude (radians)') 
    plt.xlabel('$\Omega_{D}$ (rad/s)') 

    print Amp 
plt.show()   

Answer 1

你的問題是與縮進。的代碼的該部分正被運行作爲「外」 for循環中，這意味着它是運行僅「AMP」的最後值被留下的一部分，當內部的「for」循環完成：

w = np.sqrt(-g/l) # This is the natural frequency 

#Calculate the analytic soln 
    exact_solution = F_D/np.sqrt((w**2 - Omega_D**2)**2 + (i * Omega_D)**2) 

#plot num and exact solns together 
    plt.plot(Omega_D, exact_solution) 
    plt.plot(Omega_D, Amp) 
    plt.title('q = ') 
    plt.ylabel('Amplitude (radians)') 
    plt.xlabel('$\Omega_{D}$ (rad/s)') 

    print Amp 

您需要縮進一個級別以便它作爲內部「for」循環的一部分運行。

此外，該行沒有做你想要什麼：

 Amp == np.abs(solution[0]) 

您試圖分配np.abs（溶液[0]）放大器，而是你是否np.abs測試（解決方案[0]）等於Amp（但隨後扔掉測試結果）。此行應爲：

 Amp = np.abs(solution[0])