2

我想编写一个程序，将二阶微分方程转换为两个常微分方程，但我不知道如何在 Python 中做到这一点。

我收到很多错误，请帮助正确编写代码。

from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
import numpy as np

N = 30          # Number of coupled oscillators.
alpha=0.25
A = 1.0 

# Initial positions.
y[0] = 0        # Fix the left-hand side at zero.
y[N+1] = 0      # Fix the right-hand side at zero.
# The range(1,N+1) command only prints out [1,2,3, ... N].
for p in range(1, N+1):    # p is particle number.
    y[p] = A * np.sin(3 * p * np.pi /(N+1.0))

####################################################
# Initial velocities.
####################################################
v[0]   = 0               # The left and right boundaries are
v[N+1] = 0               # clamped and don't move.
# This version sets them all the particle velocities to zero.
for p in range(1, N+1):
    v[p] = 0

w0 = [v[p], y[p]]
def accel(t,w):
    v[p], y[p] = w
    global a
    a[0] = 0.0
    a[N+1] = 0.0

    # This version loops explicitly over all the particles.
    for p in range(1,N+1):
        a[p] = [v[p], y(p+1)+y(p-1)-2*y(p)+ alpha * ((y[p+1] - y[p])**2 - (y[p] - y[p-1])**2)] 
    return a

   
duration = 50
t = np.linspace(0, duration, 800) 
abserr = 1.0e-8
relerr = 1.0e-6

solution = solve_ivp(accel, [0, duration], w0, method='RK45', t_eval=t, 
    vectorized=False, dense_output=True, args=(),  atol=abserr, rtol=relerr)
4

1 回答 1

1

大多数通用求解器不处理结构化状态对象。它们只是使用平面数组作为状态空间点的表示。从初始点的构造来看，您似乎更喜欢状态空间排序

  [ v[0], v[1], ... v[N+1], y[0], y[1], ..., y[N+1] ]

这允许简单地拆分两者并从速度和加速度数组中组装导数向量。

让我们在小函数中保持简单和独立的功能

a = np.zeros(N+2)
def accel(y):
    global a ## initialized to the correct length with zero, avoids repeated allocation
    a[1:-1] = y[2:]+y[:-2]-2*y[1:-1] + alpha*((y[2:]-y[1:-1])**2-(y[1:-1]-y[:-2])**2)
    return a

def derivs(t,w):
    v,y = w[:N+2], w[N+2:]
    return np.concatenate([accel(y), v])

或保持避免分配的主题

dwdt = np.zeros(2*N+4)
def derivs(t,w):
    global dwdt
    v,y = w[:N+2], w[N+2:]
    dwdt[:N+2] = accel(y)
    dwdt[N+2:] = v
    return dwdt

现在你只需要设置

w0=np.concatenate([v,y])

快速找到更有趣的错误类别。

于 2021-03-07T14:38:09.900 回答