Simulate a coupled ordinary differential equation

Question

I want to write a program which turns a 2nd order differential equation into two ordinary differential equations but I don't know how I can do that in Python.

I am getting lots of errors, please help in writing the code correctly.

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)

Answer 1

Most general-purpose solvers do not do structured state objects. They just work with a flat array as representation of the state space points. From the construction of the initial point you seem to favor the state space ordering

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

This allows to simply split both and to assemble the derivatives vector from the velocity and acceleration arrays.

Let's keep things simple and separate functionality in small functions

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])

or keeping the theme of avoiding allocations

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

Now you only need to set

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

to rapidly get to a more interesting class of errors.