Как исправить ошибку индекса в дифференциальном уравнении?

Question

Я пытаюсь создать программу, которая решает систему массовых пружин-демпферов с использованием обратного дифференцирования, единственная проблема заключается в том, что я сталкиваюсь с ошибкой индекса, которую я не знаю, как ее решить:

import numpy as np 
import matplotlib.pyplot as plt

def MSD_Solver(m,b,K):
    #input: m = mass, b = damping ratio, K = spring constant
    #output: (t,x) time vs position

    tinitial = 0
    tfinal = 15
    step = .005

    t = np.linspace(tinitial,tfinal,step)

    x = np.zeros_like(t)

    x[0]=0
    x[1]=0
    for k in range (len(t)-1):            # extra element so subtract by 1
        x[k] = (t**2)/((m+b)*t+(t**2)*k) + (x[k-2](-m))/((m+b)*t+(t**2)*k) + (x[k-1]((2*m)+(b*t)))/((m+b)*t+(t**2)*k)
    return plt.plot(t,x)

print(MSD_Solver(1,.5,5)),MSD_Solver(1,1,5),MSD_Solver(1,2,5)
plt.show()

Answer 1

Вы хотите, вероятно (?), Использовать коэффициенты разности первого и второго порядка для дискретизации

m*x''(t) + b*x'(t) + K*x(t) = 1

до

m*(x[j+1]-2*x[j]+x[j-1]) + 0.5*dt*b*(x[j+1]-x[j-1]) + dt^2*K*x[j] = dt**2

так что

x[j+1] = ( dt**2 + (2*m-K*dt**2)*x[j] - (m-0.5*dt*b)*x[j-1] ) / (m+0.5*dt*b)

В коде

def MSD_Solver(m,b,K):
    #input: m = mass, b = damping ratio, K = spring constant
    #output: (t,x) time vs position

    tinitial = 0
    tfinal = 15
    step = .005

    t = np.arange(tinitial,tfinal,step)
    x = np.zeros_like(t)

    dt = t[1]-t[0]  # use the actual time step
    x[0:2] = [ 0, 0]
    for j in range(1,len(t)-1):
        x[j+1] = ( dt**2 + (2*m-K*dt**2)*x[j] - (m-0.5*dt*b)*x[j-1] ) / (m+0.5*dt*b)
    return t,x

t,x = MSD_Solver(1,.5,5)        
plt.plot(t,x); plt.show();

Answer 2

Документ linspace показывает, что третьим аргументом является количество элементов, а не шаг. Ваше значение step было усечено до 0, поэтому возвращаемый массив для t был пуст. В результате x не имеет элементов, а x[0] находится вне диапазона.

Попробуйте это:

tinitial = 0
tfinal = 15
step = .005
num = (tfinal - tinitial) / step + 1

t = np.linspace(tinitial,tfinal,num)

Это приведет вас к семантическим ошибкам в ваших сложных вычислениях.