for循环中的for循环:钟摆建模

For loop within for loop: pendulum modeling

我是编程新手,正在处理一个相对复杂的问题。我正在模拟一个线性受力摆,需要弄清楚摆动幅度(以弧度为单位)如何取决于摩擦力(q 值)和驱动力频率(Omega_D)。所以,我想我需要一个 for 循环内的 for 循环,因为我需要为 3 个 q 值绘制 Omega_D vs Amplitude 的图。因此,对于 q 迭代 3 次,对于 Omega_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()           

您的问题与缩进有关。这部分代码被 运行 作为 "outer" 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(solution[0]) 分配给 Amp,但您正在测试 np.abs(solution[0]) 是否等于 Amp(但随后抛出测试结果)。这一行应该是:

        Amp = np.abs(solution[0])