标题：结构受激振动-阻尼系数的测定
这一实际工作的目的是确定阻尼系数，以防止过度振动的结构受到循环载荷。
初始值问题：
通过求解牛顿第二定律给出的微分方程，得到位置、速度、加速度，选择每周期至少500点的计算步长，用小组指定的方法求解初值问题。获得期望阻尼系数：
定义一个函数“f（x）”，其中“x”表示“C”的值，“f（x）”表示所需的振幅。函数的值通过运行初始值问题算法获得。阻尼系数“C”的值应使用三个有效的正确数字获得。报告应遵循以下结构：B.导言，包括目标和工作总结（最多1页）C.工作进展情况;酌情包括图表D.结论，特别是关于阻尼系数值的结论E.附件一：源代码-您可以附加Python文件或提供到开发它的协作平台的链接。
将报告和代码上传到校园平台上的相应任务。
以下是定义的Bashford方法的方程：
x1c 0d1x的数据
价值观：

• K = 615 N/mm
• m = 9.7吨
• Fm = 850 N
• xadm = 1.25 mm
• P = 0.7890929760279801 s

产品特点：
该结构被赋予一个从0到P/10的力，然后9 P/10 s没有任何力，并在P上再次施加。如果xMax计算< Xadm返回C. C应使用正割法计算。但在这一刻，我只是试图使不同的图形，我得到了这个结果与此代码：

from decimal import Decimal
from typing import List

import numpy as np
import matplotlib.pyplot as plt

m = Decimal(9700)
K = Decimal(615)
F = Decimal(850)
x_adm = Decimal(0.00125)
P = Decimal(0.7890929760279801)
c = Decimal(1) #test value

def force_periodic(x, P, F):  # gives the value of F in x time
    normalized_x = Decimal(x) % P  # Normalize x to the range [0, P)
    if x < 0: return Decimal(0)
    if normalized_x < P / 10:
        return F
    else:
        return Decimal(0)

class Vibrations:
    def __init__(self, m, K, x_adm, P, F):
        self.m = m
        self.K = K
        self.x_adm = x_adm
        self.P = P
        self.F = F

    t_points = np.linspace(0, float(2 * P), 1000)
    F_points = [force_periodic(xi, P, F) for xi in t_points]
    x_points: List[Decimal] = [Decimal(0)]
    v_points: List[Decimal] = [Decimal(0)]
    a_points = [0]

    def f(self, i):
        if i < 0: return 0
        if i == 0: return F / m
        return self.F_points[i] / m - K * self.x_points[i] - c * self.v_points[i]

    def Adams_BashfordSpeed(self, i):
        if i <= 0:
            self.v_points.append((P / 1000) * (3 * self.f(i - 1) - self.f(i - 2)))
        else:
            self.v_points.append(self.v_points[i - 1] + (P / 1000) * (3 * self.f(i - 1) - self.f(i - 2)))

    def Adam_BashfordXpos(self, i):
        if i <= 0:
            self.x_points.append((P / 1000) * (3 * self.v_points[i - 1] - self.v_points[i - 2]))
        else:
            self.x_points.append(self.x_points[i - 1] + (P / 1000) * (3 * self.v_points[i - 1] - self.v_points[i - 2]))

    def calculate(self):
        for x in range(len(self.t_points) - 1):
            self.Adams_BashfordSpeed(x)
            self.Adam_BashfordXpos(x)

    def graph(self):
        plt.plot(self.t_points, self.F_points, label='Force')
        plt.xlabel('Time')
        plt.ylabel('Force')
        plt.legend()
        plt.show()
        plt.plot(self.t_points, self.v_points, color='blue', label='V(t)',drawstyle='steps')
        plt.plot(self.t_points, self.x_points, color='Red', label='X(t)',drawstyle='steps')
        plt.legend()
        plt.show()

if __name__ == '__main__':
    vibrations = Vibrations(m, K, x_adm, P, F)
    vibrations.calculate()
    print(len(vibrations.t_points))
    print(len(vibrations.x_points))
    print(len(vibrations.v_points))
    vibrations.graph()

字符串

F(t)

的

V(t) & X(t)

的
在V（T）图上，我不知道为什么它会产生蓝色区域。这是主要的问题，我不知道它是没有定义一个准时的值，还是绘制了一个糟糕的图。
我尝试了所有的代码改变变量，我不知道是什么wroken逻辑或代码。