Отклонение в решениях (дифференциальные уравнения) с использованием odeint против Runge-kutta-4th ⇐ Python

[Anonymous]

Я моделирую связанную пружинную систему: два объекта с массами M1 и M2, и связаны с помощью пружин с констант пружины K1 и K2, с демпфирующим D1 и D2.# Use ODEINT to solve the differential equations defined by the vector field
from scipy.integrate import odeint
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np

def vectorfield(w, t, p):
"""
Defines the differential equations for the coupled spring-mass system.

Arguments:
w : vector of the state variables:
w = [x1,y1,x2,y2]
t : time
p : vector of the parameters:
p = [m1,m2,k1,k2,b1,b2]
"""
x1, y1, x2, y2 = w
m1, m2, k1, k2, b1, b2 = p

f = [y1,
(-d1 * y1 - k1 * x1 + k2 * (x2 - x1 )+ d2 * (y2-y1)) / m1,
y2,
(-d2 * (y2-y1) - k2 * (x2 - x1 )) / m2]

return f

# Parameter values
tau=1.02
f_1=0.16
f_2=tau*f_1

m1 = 2000000 # mass1
m2 = 20000 # mass2

# Spring constants
k1 = m1*pow(2*np.pi*f_1,2)
k2 = m2*pow((2*np.pi*f_2),2)

# damping

d1 = (0.04/2/np.pi)*2*pow(k1*m1,0.5)

d_d1=6000
l_p=9.81/pow(2*np.pi*f_2,2)
b=0.3
d2=d_d1*pow((l_p-b)/l_p,2)

# Initial conditions
# x1 and x2 are the initial displacements; y1 and y2 are the initial velocities
x1 = 0.5
y1 = 0.0
x2 = 0.25
y2 = 0.0

# ODE solver parameters
abserr = 1.0e-8
relerr = 1.0e-6
stoptime = 100.0
numpoints = 5001

# Create the time samples for the output of the ODE solver.
# I use a large number of points, only because I want to make
# a plot of the solution that looks nice.
t = [stoptime * float(i) / (numpoints - 1) for i in range(numpoints)]

# Pack up the parameters and initial conditions:
p = [m1, m2, k1, k2, d1, d2]
w0 = [x1, y1, x2, y2]

# Call the ODE solver.
wsol = odeint(vectorfield, w0, t, args=(p,),
atol=abserr, rtol=relerr)

# convert zip object to tuple object
temp=tuple(zip(t, wsol[:,0],wsol[:,1],wsol[:,2],wsol[:,3]))
# convert tulpe object to list then to dataframe
df=pd.DataFrame(list(map(list,temp)))

# Create plots with pre-defined labels.
fig, ax = plt.subplots()
ax.plot(df.loc[:,0], df.loc[:,1],label='displacement object1')

legend = ax.legend(loc='upper right', shadow=None, fontsize='small')
ax.set_xlabel('time ', fontdict=None, labelpad=None, loc='center')
ax.set_ylabel('pos [m] or V [m/s]', fontdict=None, labelpad=None, loc='center')
< /code>
methode 2: Использование Runge-Kutta-4-й порядок (запрограммирован сам) для решения дифференциальных уравнений < /p>
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import matplotlib.pyplot as plt

# Parameter values
tau=1.02
f_1=0.16
f_2=tau*f_1

m1 = 2000000 # mass1
m2 = 20000 # mass2

# Spring constants
k1 = m1*pow(2*np.pi*f_1,2)
k2 = m2*pow((2*np.pi*f_2),2)

# damping
d1 = (0.04/2/np.pi)*2*pow(k1*m1,0.5)

d_d1=6000
l_p=9.81/pow(2*np.pi*f_2,2)
b=0.3
d2=d_d1*pow((l_p-b)/l_p,2)

def system1(x1,y1,x2,y2):
return (-d1 * y1 - k1 * x1 + k2 * (x2 - x1 )+ d2 * (y2-y1)) / m1

def system2(x1, y1, x2, y2):
return (-d2 * (y2-y1) - k2 * (x2 - x1)) / m2

def runge_kutta_4_sys1(f, x1, y1, x2, y2, h):
k1 = f(x1,y1,x2,y2)
k2 = f(x1,y1+k1*h/2,x2,y2)
k3 = f(x1,y1+k2*h/2,x2,y2)
k4 = f(x1,y1+k3*h,x2,y2)

temp1= y1
temp2= y1+k1*h/2
temp3= y1+k2*h/2
temp4= y1+k3*h

displace_solution = x1 + (temp1 + 2 * temp2 + 2 * temp3 + temp4)*h / 6
velocity_solution = y1 + (k1 + 2 * k2 + 2 * k3 + k4)*h / 6
acc_solution = f(x1,y1,x2,y2)

return displace_solution, velocity_solution, acc_solution

def runge_kutta_4_sys2(f, x1, y1, x2, y2, h):
k1 = f(x1,y1,x2,y2)
k2 = f(x1,y1,x2,y2+k1*h/2)
k3 = f(x1,y1,x2,y2+k2*h/2)
k4 = f(x1,y1,x2,y2+k3*h)

temp1= y2
temp2= y2+k1*h/2
temp3= y2+k2*h/2
temp4= y2+k3*h

displace_solution = x2 + (temp1 + 2 * temp2 + 2 * temp3 + temp4)*h / 6
velocity_solution = y2 + (k1 + 2 * k2 + 2 * k3 + k4)*h / 6
acc_solution = f(x1,y1,x2,y2)

return displace_solution, velocity_solution, acc_solution

x1 = 0.5
y1 = 0.0
x2 = 0.25
y2 = 0.0

h = 0.02
numpoints = 5000
time = 0

temp=0

df = pd.DataFrame(index=range(1+numpoints),columns=range(7))

df.iloc[0]=[time,x1,y1,temp,x2,y2,temp]

for i in range(numpoints):
x1, y1, acc1 = runge_kutta_4_sys1(system1, x1, y1, x2, y2, h)
x2, y2, acc2 = runge_kutta_4_sys2(system2, x1, y1, x2, y2, h)
time=time+h
df.iloc[i+1]=[time,x1,y1,acc1,x2,y2,acc2]
df.iloc[i,3]=acc1
df.iloc[i,6]=acc2

# Create plots with pre-defined labels.
fig, ax = plt.subplots()
ax.plot(df.loc[:,0], df.loc[:,1],label='displacement object1')

legend = ax.legend(loc='upper right', shadow=None, fontsize='small')
ax.set_xlabel('time ', fontdict=None, labelpad=None, loc='center')
ax.set_ylabel('pos [m] or V [m/s]', fontdict=None, labelpad=None, loc='center')
< /code>
Явно имеется отклонение в растворах TH между методом 1 и методом 2.

< /p>
Я предполагаю, что Deviatia inaccurtate-in. Мое предположение правильно? Как установить или улучшить сопротивление для Runge-kutta-methode?
Большое спасибо!>

Подробнее здесь: https://stackoverflow.com/questions/794 ... -kutta-4th