मैं दूसरे क्रम के ओडीई की प्रणाली को हल करने के लिए scipy.integrate.solve_ivp का उपयोग कर रहा हूं।

यहां उपलब्ध कोड को मानते हुए (ध्यान दें कि यह odeint का उपयोग करता है, लेकिन समान दृष्टिकोण का उपयोग करता है)। https: //scipy-cookbook.readthedocs.io/items/CoupledSpringMassSystem.html

अब, कुछ पैरामीटर (बी 1 और बी 2) समय और समाधान-निर्भर बनाते हुए, कहें: उदा।

b1 = 0.8 * x1

if (x2 >= 1.0):
    b2 = 20*x2
else:
    b2 = 1.0

Solution_ivp समाधान आसानी से मुद्रित और प्लॉट किए जा सकते हैं।

मैं परिभाषित फ़ंक्शन के भीतर print(b1) डालकर समय/समाधान-निर्भर पैरामीटर का हिस्सा प्रिंट करने में भी सक्षम हूं। लेकिन मैं इस आउटपुट से भ्रमित हूं। मैं जो चाहता हूं वह समाधान टाइमस्टैम्प पर हल किए जा रहे परिभाषित फ़ंक्शन के भीतर निर्धारित समय/समाधान-निर्भर पैरामीटर को प्रिंट और प्लॉट करना है: उदा। print(t, b1).

उम्मीद है कि यह समझ में आता है। कोई संकेत बहुत अच्छा होगा।

चियर्स

नमूना कोड (ओडीइंट और बी1 और बी2 के साथ गतिशील पैरामीटर के रूप में):

    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,L1,L2,b1,b2]
    """
    x1, y1, x2, y2 = w
    m1, m2, k1, k2, L1, L2, = p
    
    # Friction coefficients
    b1 = 0.8 * x1
    
    if (x2 >= 1.0):
        b2 = 20*x2
    else:
        b2 = 1.0
    
    # Create f = (x1',y1',x2',y2'):
    f = [y1,
         (-b1 * y1 - k1 * (x1 - L1) + k2 * (x2 - x1 - L2)) / m1,
         y2,
         (-b2 * y2 - k2 * (x2 - x1 - L2)) / m2]
    
    print('b1', b1)
    # print('b2', b2)
    
    return f

# Use ODEINT to solve the differential equations defined by the vector field
from scipy.integrate import odeint

# Parameter values
# Masses:
m1 = 1.0
m2 = 1.5
# Spring constants
k1 = 8.0
k2 = 40.0
# Natural lengths
L1 = 0.5
L2 = 1.0

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

# ODE solver parameters
abserr = 1.0e-8
relerr = 1.0e-6
stoptime = 10.0
numpoints = 101

# Create the time samples for the output of the ODE solver.
t = [stoptime * float(i) / (numpoints - 1) for i in range(numpoints)]

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

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

# Print Output
print(t)
# print(wsol)
# print(t, b1)
# print(t, b2)


# Plot the solution that was generated
from numpy import loadtxt
from pylab import figure, plot, xlabel, grid, legend, title, savefig
from matplotlib.font_manager import FontProperties

figure(1, figsize=(6, 4.5))

xlabel('t')
grid(True)
lw = 1

plot(t, wsol[:,0], 'b', linewidth=lw)
plot(t, wsol[:,1], 'g', linewidth=lw)
plot(t, wsol[:,2], 'r', linewidth=lw)
plot(t, wsol[:,3], 'c', linewidth=lw)
# plot(t, b1, 'k', linewidth=lw)

legend((r'$x_1$', r'$x_2$',r'$x_3$', r'$x_4$'), prop=FontProperties(size=16))
title('Mass Displacements for the\nCoupled Spring-Mass System')
# savefig('two_springs.png', dpi=100)
0
etienne 9 अक्टूबर 2020, 15:56

1 उत्तर

सबसे बढ़िया उत्तर

नया ओडीई सॉल्वर एपीआई

का उपयोग करने के लिए अपने कोड को अनुकूलित करना solve_ivp (जो ODE को हल करने के लिए आधुनिक scipy API है) हम सिस्टम को बिना दखल के हल कर सकते हैं:

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

numpoints = 1001
t0 = np.linspace(0, stoptime, numpoints)

def func(t, x0, beta):
    return vectorfield(x0, t, beta)

#sol1 = odeint(vectorfield, w0, t0, args=(p,), atol=abserr, rtol=relerr)
sol2 = solve_ivp(func, [t0[0], t0[-1]], w0, args=(p,), t_eval=t0, method='LSODA')

नया API पुराने API odeint. solve_ivp द्वारा लौटाई गई समाधान वस्तु कार्य निष्पादन के बारे में कई अतिरिक्त जानकारी प्रदान करती है:

  message: 'The solver successfully reached the end of the integration interval.'
     nfev: 553
     njev: 29
      nlu: 29
      sol: None
   status: 0
  success: True
        t: array([ 0.  ,  0.01,  0.02, ...,  9.98,  9.99, 10.  ])
 t_events: None
        y: array([[ 0.5       ,  0.50149705,  0.50596959, ...,  0.60389634,
         0.60370222,  0.6035087 ],
       [ 0.        ,  0.29903815,  0.59476327, ..., -0.01944383,
        -0.01937904, -0.01932493],
       [ 2.25      ,  2.24909287,  2.24669965, ...,  1.62461438,
         1.62435635,  1.62409906],
       [ 0.        , -0.17257692, -0.29938513, ..., -0.02583719,
        -0.02576507, -0.02569245]])
 y_events: None

गतिशील मापदंडों को पकड़ना

अगर मैं आपके प्रश्न को सही ढंग से समझता हूं, तो आप समय श्रृंखला के रूप में भी अपने घर्षण पैरामीटर रखना चाहते हैं।

सॉल्वर निष्पादन का पालन करें

इसे कॉलबैक फ़ंक्शन या पासिंग का उपयोग करके प्राप्त किया जा सकता है ODE फ़ंक्शन का मूल्यांकन करते समय एक अतिरिक्त डेटा संरचना और स्टोर गुणांक मान।

यदि आप सॉल्वर निष्पादन को ट्रैक करना चाहते हैं तो यह सहायक होगा, लेकिन परिणाम वांछित समाधान समय पैमाने के बजाय सॉल्वर एकीकरण ग्रिड पर निर्भर होगा।

वैसे भी, यदि आप ऐसा करना चाहते हैं, तो बस अपना कार्य अनुकूलित करें:

def vectorfield2(w, t, p, store):
   # ...
   store.append({'time': t, 'coefs': [b1, b2]})
   # ...

c = []
sol1 = odeint(vectorfield2, w0, t0, args=(p, c), atol=abserr, rtol=relerr)

जैसा कि आप देख सकते हैं, परिणाम समाधान समय पैमाने का पालन नहीं करता है बल्कि इसके बजाय सॉल्वर एकीकरण ग्रिड का पालन करता है:

[{'time': 0.0, 'coefs': [0.4, 45.0]},
 {'time': 3.331483023263848e-07, 'coefs': [0.4, 45.0]},
 {'time': 3.331483023263848e-07,
  'coefs': [0.4000000000026638, 44.999999999955605]},
 {'time': 6.662966046527696e-07,
  'coefs': [0.4000000000053275, 44.99999999991122]},
 {'time': 6.662966046527696e-07,
  'coefs': [0.40000000000799124, 44.99999999986683]},
 {'time': 0.0001385450759485236,
  'coefs': [0.4000002303394594, 44.99999616108455]}, ...]

राज्य समारोह के रूप में गुणांक Co

दूसरी ओर, आपके गुणांक केवल सिस्टम की स्थिति पर निर्भर करते हैं, इसका मतलब है कि हम इसकी गणना ODE फ़ंक्शन के अंदर कर सकते हैं (हर बार सॉल्वर इसे कॉल करता है) और बाद में जब समय पर निर्भर समाधान ज्ञात होते हैं।

मैं आपके सिस्टम को निम्नानुसार फिर से लिखने का सुझाव दूंगा:

def coefs(x):
    # Ensure signature:
    if isinstance(x, (list, tuple)):
        x = np.array(x)
    if len(x.shape) == 1:
        x = x.reshape(1, -1)
    # Compute coefficients
    b1 = 0.8 * x[:,0]
    b2 = 20. * x[:,2]
    q2 = (x[:,2] < 1.0)
    b2[q2] = 1.0
    return np.stack([b1, b2]).T

def system(t, w, p):
    x1, y1, x2, y2 = w
    m1, m2, k1, k2, L1, L2, = p
    b = coefs(w)
    return [
        y1, (-b[:,0] * y1 - k1 * (x1 - L1) + k2 * (x2 - x1 - L2)) / m1,
        y2, (-b[:,1] * y2 - k2 * (x2 - x1 - L2)) / m2,
    ]

sol3 = solve_ivp(system, [t0[0], t0[-1]], w0, args=(p,), t_eval=t0, method='LSODA')

ध्यान दें कि coefs विधि का इंटरफ़ेस सरल अवस्था (सूची या वेक्टर) और समय पर समाधान (मैट्रिक्स) पर संचालित करने के लिए डिज़ाइन किया गया है।

फिर, समय समाधान प्रदर्शित करना आसान है:

fig, axe = plt.subplots()
axe.plot(t0, sol3.y.T)
axe.set_title("Dynamic System: Coupled Spring Mass")
axe.set_xlabel("Time, $t$")
axe.set_ylabel("System Coordinates, $x_i(t)$")
axe.legend([r'$x_%d$' % i for i in range(sol3.y.T.shape[1])])
axe.grid()

enter image description here

और गुणांक:

fig, axes = plt.subplots(2,1, sharex=True, sharey=False)
axes[0].plot(t0, coefs(sol3.y.T)[:,0])
axes[1].plot(t0, coefs(sol3.y.T)[:,1])
axes[0].set_title("Dynamic System: Coupled Spring Mass")
axes[1].set_xlabel("Time, $t$")
axes[0].set_ylabel("$b_0(t)$")
axes[1].set_ylabel("$b_1(t)$")
for i in range(2):
    axes[i].grid()

enter image description here

जो आप हासिल करना चाहते थे, उससे अधिक संभावना है।

0
jlandercy 10 अक्टूबर 2020, 13:52