График падения мяча с течением времени с питоном

Я пытаюсь написать код, который будет отображать симуляцию падения мяча с высоты h и строить график положения во времени с использованием кинематических уравнений y = y_0. Мой код таков:

из matplotlib.pylab import show, xlabel, ylabel, scatter, график из пустого импорта numpy

def drop():

    """
    This function calculates and creates arrays for the velocity at eac time interval as well as the position and plots it. Assuming no drag. 
    """
    #Define the constants in the problem
    h_0 = 10
    g = -9.8 #gravitational constant N
    dt = 0.1 #timestep

    #Now need to create arrays to hold the positins, time, and velocities
    vel = empty(1000,float)
    time = empty(1000,float)
    y = empty(1000,float)

    time[0] = 0
    vel[0] = 0
    y[0] = h_0

    #code for the kinematic equations for the calculation of time, velocity and position
    for i in range of (1000-1):
        time[i+1] = time[i] + dt
        vel[i+1] = vel[i] + (g * dt)
        y[i+1] = time[i] + (vel[i+1] * dt)

        if y[i] > 0:
        #ensures that the graph will not keep going when the ball hits the ground
            break


    plot(time,y, '.')
    xlabel("Time(s)")
    ylabel("Position")
    show()

Однако мой график отображает три точки по одной в каждом углу графика, когда он должен выглядеть как кривая, и мой график меняется каждый раз, когда он не должен, поскольку ни одна из переменных не меняется.


python arrays plot scatter kinematics
person Anny    schedule 08.01.2016    source источник

Ответы (1)


arrow_upward
1
arrow_downward

Хорошо, давайте уберем синтаксическую ошибку. for i in range of (1000-1) на самом деле for i in range(1000-1), но я предполагаю, что это была опечатка с вашей стороны, так как вы могли запустить код.

Ваши уравнения движения неверны.

y[i+1] = time[i] + (vel[i+1] * dt)

# should be
y[i+1] = y[i] + (vel[i] * dt)

Ваше условие выхода из симуляции также ошибочно.

if y[i] > 0:

# you have to stop when the height becomes until negative
if y[i+1] < 0:

Ваши ошибки до сих пор означают, что вы выйдете из цикла после одной итерации, фактически не изменив свой массив y. Последняя проблема начинается здесь. numpy.empty() создает массив без инициализации значений. Это означает, что исходными значениями будут все, что находится в памяти в этот момент. Если вы напечатаете y после разрыва цикла, вы можете заметить, что большинство значений равны 0, а некоторые очень малы, но не близки к 0, например. 3.18377034э-308. Поскольку они являются самыми высокими значениями в вашем массиве, они будут масштабировать ваш график в соответствии со своим диапазоном. Но поскольку это произвольные значения, каждый раз, когда вы запускаете код, он будет выдавать разные числа.

У вас есть два варианта исправить это. Либо используйте numpy.zeros(), либо нанесите только первые значения y[:i], которые являются точкой в ​​петле, которую вы прерываете при ударе о землю.

Поскольку у нас есть аналитическое решение уравнений в вашей задаче, вы можете избавиться от циклов и векторизовать все с помощью массивов. Мы можем решить уравнение смещения относительно t (квадратичное), чтобы узнать, когда мы коснемся земли. Затем мы инициализируем массив времени и используем его для вычисления смещения (скорость необязательна).

import numpy as np
import matplotlib.pyplot as plt

def time_to_hit_ground(y0, v0, a):
    discriminant = np.sqrt(v0**2 - 2*a*y0)
    t1 = (-v0 - discriminant) / a
    t2 = (-v0 + discriminant) / a
    if t1 >=0:
        return t1
    return t2

def drop(y0, v0=0.0, dt=0.1, g=-9.8):
    if y0 < 0:
        print('Object is underground.')
        return
    # if you also allow the user to change `dt` and `g` from the arguments,
    # you want to check they have the correct sign.

    t_ground = time_to_hit_ground(y0, v0, g)

    t = np.arange(0, t_ground+dt, dt)
    v = v0 + g * t
    y = y0 + v0 * t + g * t**2 / 2.

    plt.plot(t, y, '.')
    plt.axvline(t_ground, color='g')
    plt.axhline(0, color='g')
    plt.xlabel('Time (s)')
    plt.ylabel('Height (m)')
    plt.show()
person Reti43    schedule 08.01.2016