Я пытаюсь написать программу, использующую систему связанных пружинных масс. Решить, используя ODE

Question

Я пытаюсь написать программу, используя систему связанных пружинных масс.Решить, используя ODE's

import numpy as np
import matplotlib.pyplot as plt

kitchen1 = 2.0
kitchen2 = 2.0
money1 = 2.0
money2 = 2.0
walk1 = 5.0
walk2 = 5.0
Lenon1 = 2.0
Lenon2 = 2.0

def f1(x1, v1, x2, v2, t):
return v1

def f2 (x1, v1, x2, v2, t):
    return ((-kitchen1/money1)*(x1-Lenon1)) + (kitchen2/money1*(x2-x1-walk1- Lenon2))

def f3(x1, v1, x2, v2, t):
return v2

def f4(x1, v1, x2, v2, t):
     return (-kitchen2/money2) * (x2-x1-walk1-Lenon2)

def rk4_gen(x1_arr, v1_arr, x2_arr, v2_arr, t_arr, h):
x1, v1, x2, v2, t  = t_arr[0], x1_arr[0], v1_arr[0], x2_arr[0], v2_arr[0]

for i in range(1,len(t)):
    t, h = t[i-1], t[i]-t[i-1]  

    k11 = h*f1(x1, v1, x2, v2, t)
    k12 = h*f2(x1, v1, x2, v2, t)
    k13 = h*f3(x1, v1, x2, v2, t)
    k14 = h*f4(x1, v1, x2, v2, t)


    k21 = h*f1(x1+k11/2.0, v1+k12/2.0,x2+k13/2.0,v2+k14/2.0, t+h/2.0)
    k22 = h*f2(x1+k11/2.0, v1+k12/2.0,x2+k13/2.0,v2+k14/2.0, t+h/2.0)
    k23 = h*f3(x1+k11/2.0, v1+k12/2.0,x2+k13/2.0,v2+k14/2.0, t+h/2.0)
    k24 = h*f4(x1+k11/2.0, v1+k12/2.0,x2+k13/2.0,v2+k14/2.0, t+h/2.0)

    k31 = h*f1(x1+k21/2.0, v1+k22/2.0,x2+k23/2.0,v2+k24/2.0, t+h/2.0)
    k32 = h*f2(x1+k21/2.0, v1+k22/2.0,x2+k23/2.0,v2+k24/2.0, t+h/2.0)
    k33 = h*f3(x1+k21/2.0, v1+k22/2.0,x2+k23/2.0,v2+k24/2.0, t+h/2.0)
    k34 = h*f4(x1+k21/2.0, v1+k22/2.0,x2+k23/2.0,v2+k24/2.0, t+h/2.0)

    k41 = h*f1(x1+k31, v1+k32,x2+k33,v2+k34, t+h)
    k42 = h*f2(x1+k31, v1+k32,x2+k33,v2+k34, t+h)
    k43 = h*f3(x1+k31, v1+k32,x2+k33,v2+k34, t+h)
    k44 = h*f4(x1+k31, v1+k32,x2+k33,v2+k34, t+h)


    x1 = x1 + (k11 + 2*k12 + 2*k13 + k14)/6.0
    v1 = v1 + (k21 + 2*k22 + 2*k23 + k24)/6.0
    x2 = x2 + (k31 + 2*k32 + 2*k33 + k34)/6.0
    v2 = v2 + (k41 + 2*k42 + 2*k43 + k44)/6.0

    x1_arr[i], v1_arr[i], x2_arr[i], v2_arr[i] = x1, v1, x2, v2  

return x1_arr, v1_arr , x2_arr , v2_arr, t_arr 

---

def grafik_ciz(t,y,x_baslik,y_baslik,baslik):
plt.figure(figsize = [12, 9])
plt.plot(t,y)
plt.xlabel(x_baslik)
plt.ylabel(y_baslik)
plt.title(baslik)
plt.grid()
plt.show()

ti = 0 
tf = 20 
x1 = 2 
v1 = 0
x2 = 15 
v2 = 0
h = 0.5

t = np.arange(ti,tf+h,h) 
x1 = np.zeros(len(t))
v1 = np.zeros(len(t))
x2 = np.zeros(len(t))
v2 = np.zeros(len(t))


x1[0] = x1i
v1[0] = v1i
x2[0] = x2i
v2[0] = v2i

x1 = rk4_gen(x1_arr, v1_arr, x2_arr, v2_arr, t_arr, h
print x1

---

Используя метод Рунге-Кутты 4-го порядка

Мне нужно построить график, показывающий x и y как функцию времениот t = 0 до t = 20.

ошибка m.ValueError: слишком много значений для распаковки

Вот мой код, но пока ничего не отображается на графике.

Я думаю, эта ошибка t = np.arange(ti,tf+h,h).
Мне нужно построить график, показывающий как x1, так и v1 как функцию времени от t = 0 до t = 20.

Answer 1

Помимо множества ошибок при правильной реализации алгоритма, использовании согласованных интерфейсов функций, использовании правильных имен функций и т. Д.

Ваша проблема в том, что вы не распаковываете массивы, переданные rk4_gen при вычислении отдельных шагов RK4. Что вы можете сделать без слишком большого количества изменений, это переименовать параметры, как в

def rk4_gen(t_arr, x1_arr, v1_arr, x2_arr, v2_arr, h):

(почему вы передаете h? На каждом шаге вы можете вычислить правильное h=t[i]-t[i-1]). Затем инициализируйте работающие переменные

    t, x1, v1, x2, v2 = t_arr[0], x1_arr[0], v1_arr[0], x2_arr[0], v2_arr[0]

выполнить (исправленные) шаги RK4 и в конце каждого шага сохранить вычисленное значение в массиве

def rk4_gen(x1_arr, v1_arr, x2_arr, v2_arr, t_arr, h):
    t, x1, v1, x2, v2 = t_arr[0], x1_arr[0], v1_arr[0], x2_arr[0], v2_arr[0]

    for i in range(1,len(t_arr)):
        # Do the RK4 step, it is from i-1 to i
        t, h = t_arr[i-1], t_arr[i]-t_arr[i-1]
        k11 = h*f1(...)
        ...
        v2 = v_2+(k14+2*k24+2*k34+k44)/6.0
        x1_arr[i], v1_arr[i], x2_arr[i], v2_arr[i] = x1, v1, x2, v2

    return x1_arr, v1_arr , x2_arr , v2_arr, t_arr 

Также обратите внимание, что в вашем конечном расчете значения k расположены в транспонированном порядке.

---

Вы можете решить эту проблему гораздо более компактным способом, используя scipy.integrate методы, как в

from scipy.integrate import solve_ivp

def spring_ode(t,u):
    x1,v1,x2,v2 = u
    a1 = ((-kitchen1/money1)*(x1-Lenon1)) + (kitchen2/money1*(x2-x1-walk1- Lenon2))
    a2 = (-kitchen2/money2) * (x2-x1-walk1-Lenon2)
    return [v1,a1,v2,a2]

res = solve_ivp(spring_ode, [ti,tf], [x1i,v1i,x2i,v2i], dense_output=True, rtol=1e-12, atol=1e-14)

h=0.02
t = np.arange(ti,tf+h/2,h) 
u = res.sol(t)
x1,v1,x2,v2 = u
plt.figure(figsize = [8, 5])
plt.plot(t,x1,t,x2); plt.legend(["$x_1$","$x_2$"])
plt.grid(); plt.show()

с результирующим сюжетом