For循环内For循环：钟摆建模

Question

我是编程新手，正在研究一个相对复杂的问题。我正在模拟一个线性受迫摆，需要弄清楚摆动运动的振幅（弧度）是如何依赖于摩擦力（q值）和驱动力的频率（ωD）。所以，我想在for循环中需要一个for循环，因为我需要绘制3个q值的ωD与振幅的关系图。因此，对于q，迭代3次，对于ΩD，迭代次数更多。然而，我写的是每个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

你的问题是缩进。这部分代码是作为“outer”for循环的一部分运行的，这意味着当内部“for”循环结束时，它只运行“Amp”的最后一个值：

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