Разработка уравнения

Question

Я пытаюсь решить дифференциальное уравнение численно, и я пишу уравнение, которое даст мне массив решения для каждой точки времени.

import numpy as np 
import matplotlib.pylab as plt 

pi=np.pi 
sin=np.sin 
cos=np.cos 
sqrt=np.sqrt 
alpha=pi/4 
g=9.80665 
y0=0.0 
theta0=0.0 

sina = sin(alpha)**2 
second_term = g*sin(alpha)*cos(alpha) 

x0 = float(raw_input('What is the initial x in meters?')) 
x_vel0 = float(raw_input('What is the initial velocity in the x direction in m/s?')) 
y_vel0 = float(raw_input('what is the initial velocity in the y direction in m/s?')) 
t_f = int(raw_input('What is the maximum time in seconds?')) 

r0 = x0 
vtan = sqrt(x_vel0**2+y_vel0**2) 
dt = 1000 
n = range(0,t_f) 
r_n = r0*(n*dt) 
r_nm1 = r0((n-1)*dt) 
F_r = ((vtan**2)/r_n)*sina-second_term 
r_np1 = 2*r_n - r_nm1 + dt**2 * F_r 
data = [r0] 

for time in n: 
    data.append(float(r_np1)) 
print data 

Я не уверен, как решить уравнение для r_np1 каждый раз в диапазоне n. Я все еще новичок в Python и хотел бы, чтобы какая-то помощь поняла, как сделать что-то подобное.

Answer 1

Первый вопрос:

n = range(0,t_f) r_n = r0*(n*dt)

Здесь определяется п в виде списка и попытаться умножить список п с целым DT. Это не будет работать. Pure Python не является векторизованным языком, например NumPy или Matlab, где вы можете делать такое векторное умножение. Вы могли бы сделать эту линию работы с

n = np.arange(0,t_f) r_n = r0*(n*dt),

, но вы не должны. Вместо этого вы должны переместить все в цикле for, чтобы выполнять вычисления на каждый временной интервал. На данный момент вы делаете расчет один раз, затем добавляете тот же результат только t_f раз в список data.

Конечно, вы должны оставить свои начальные условия (которые являются ключевой частью решения ODE) OUTSIDE цикла, потому что они влияют только на первый шаг решения, а не на все из них.

Итак:

# Initial conditions 
r0 = x0 
data = [r0] 

# Loop along timesteps 
for n in range(t_f): 

    # calculations performed at each timestep 
    vtan = sqrt(x_vel0**2+y_vel0**2) 
    dt = 1000 
    r_n = r0*(n*dt) 
    r_nm1 = r0*((n-1)*dt) 
    F_r = ((vtan**2)/r_n)*sina-second_term 
    r_np1 = 2*r_n - r_nm1 + dt**2 * F_r 

    # append result to output list 
    data.append(float(r_np1)) 

# do something with output list 
print data 
plt.plot(data) 
plt.show() 

я не добавлял любой кусок кода, только переставить свои линии. Обратите внимание на то, что часть:

n = range(0,t_f) 
for time in n: 

Может быть упрощена:

for time in range(0,t_f): 

Однако использовать n в качестве временной переменной при вычислении (ранее - и ошибочно - определяется как список вместо одного номер). Таким образом, вы можете написать:

for n in range(0,t_f): 

Примечание 1: Я не знаю, если этот код является правильным математически, так как я даже не знаю, уравнение вы решаете. Теперь код запускается и дает результат - вы должны проверить, хорош ли результат.

Примечание 2: Pure Python - не лучший инструмент для этой цели. Вы должны попробовать некоторые высоко оптимизированные встроенные модули SciPy для решения ODE, так как вы уже получили подсказки в комментариях.