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