Моделируйте связанное обыкновенное дифференциальное уравнение

Я хочу написать программу, которая превращает дифференциальное уравнение 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)

python numpy python-2.7 differential-equations tightly-coupled-code
person user68207    schedule 07.03.2021    source источник
comment
Добро пожаловать в SO. Замечательный вопрос. Возможно, вы захотите рассмотреть возможность обновления до Python 3 (текущая версия 3.9.2). Это может даже помочь с вашим кодированием. Кроме того, не могли бы вы указать фактическую ошибку, которую вы получаете, в своем вопросе.   -  person VirtualScooter    schedule 07.03.2021
comment
@VirtualScooter: я вижу, что a не определен перед использованием, w0 не является плоским массивом, присвоение из кортежа одной позиции с плавающей запятой в a[p], возвращаемый тип accel не является плоским массивом, если последнее исправлено не полностью.   -  person Lutz Lehmann    schedule 07.03.2021

Ответы (1)


arrow_upward
1
arrow_downward

Большинство решателей общего назначения не создают объекты структурированного состояния. Они просто работают с плоским массивом как представление точек пространства состояний. Судя по построению начальной точки, вам кажется, что вы предпочитаете упорядочение в пространстве состояний.

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

чтобы быстро перейти к более интересному классу ошибок.

person Lutz Lehmann    schedule 07.03.2021
comment
w0 - начальное условие для дифференциального уравнения - person user68207; 07.03.2021
comment
Как обеспечить начальное состояние. - person user68207; 07.03.2021
comment
Короче w0=np.concatenate([np.zeros(N+2),np.sin(3*np.pi*np.linspace(0,1,N+2))]). - person Lutz Lehmann; 07.03.2021