使用scipy二階ODE積分

我想使用'scipy.integrate.odeint'積分二階微分方程。我的等候如下使用scipy二階ODE積分

m*x[i]''+x[i]'= K/N*sum(j=0 to N)of sin(x[j]-x[i])

我已經轉換爲兩個一階ODE，如下所示。在下面的代碼中，yinit是初始值x（0）和x'（0）的數組。我的問題是什麼應該是x（0）和x'（0）的值？

x'[i]=y[i] 
y'[i]=(-y[i]+K/N*sum(j=0 to N)of sin(x[j]-x[i]))/m 

from numpy import * 
from scipy.integrate import odeint 

N = 50 

def f(theta, t): 
    global N 
    x, y = theta 
    m = 0.95 
    K = 1.0 
    fx = zeros(N, float) 
    for i in range(N): 
     s = 0.0 
     for j in range(i+1,N): 
      s = s + sin(x[j] - x[i]) 
     fx[i] = (-y[i] + (K*s)/N)/m 
    return array([y, fx]) 

t = linspace(0, 10, 100, endpoint=False)

均勻發熱隨機數

theta = random.uniform(-180, 180, N)

集成功能f使用odeint

yinit = array([x(0), x'(0)]) 
y = odeint(f, yinit, t)[:,0] 
print (y)

來源

2014-01-10 ADK

我猜你的意思是'yinit = array（[x（0），x'（0）]'，因爲這些都是初始值，初始值是使方程組唯一解的唯一邊界條件。 – Nabla

是的，你是正確的初始值是x（0）和x'（0）。我已經糾正了他們的問題。如果我把m = 0，然後在一階ODE中，我從Uniform distribution ，將是x（0）。而在2階ODE的情況下，我知道x（0）= theta，但我怎麼能找到x'（0）？ – ADK

你可以選擇任何你想要的初始條件。

對於您的情況，您決定對所有振盪器使用x的隨機初始條件。你可以使用'y'的隨機初始條件，我猜也是如此。

上面的代碼中有一些錯誤，主要是關於如何從theta解壓縮x,y以及如何在最後重新打包它們（請參閱下面的concatenate中已更正的代碼）。請參閱concatenate的yinit。

其餘的是時尚/微小的變化。

from numpy import concatenate, linspace, random, mod, zeros, sin 
from scipy.integrate import odeint 

Nosc = 20 
assert mod(Nosc, 2) == 0 

def f(theta, _): 
    N = theta.size/2 
    x, y = theta[:N], theta[N:] 
    m = 0.95 
    K = 1.0 
    fx = zeros(N, float) 
    for i in range(N): 
     s = 0.0 
     for j in range(i + 1, N): 
      s = s + sin(x[j] - x[i]) 
     fx[i] = (-y[i] + (K * s)/N)/m 
    return concatenate(([y, fx])) 

t = linspace(0, 10, 50, endpoint=False) 

theta = random.uniform(-180, 180, Nosc) 
theta2 = random.uniform(-180, 180, Nosc) #added initial condition for the velocities of the oscillators 

yinit = concatenate((theta, theta2)) 

res = odeint(f, yinit, t) 
X = res[:, :Nosc].T 
Y = res[:, Nosc:].T

繪製系統隨時間的變化，你可以使用像

import matplotlib.pylab as plt 

fig, ax = plt.subplots() 
for displacement in X: 
    ax.plot(t, displacement) 
ax.set_xlabel('t') 
ax.set_ylabel('x') 
fig.show()

你在造型？起初，等式看起來有點像kuramoto振盪器，但後來我注意到你也有一個x[i]''期限。

請注意，在你的模型，你不必在方程中的春季學期，像在LHS的一個術語x(t)的x值收斂爲任意值：

enter image description here

來源

2014-01-10 14:42:49 gg349

非常感謝您糾正我的代碼。我仍然不確定使用速度theta2的隨機序列是否正確。是的，實際上它是Kuramoto方程，但是慣性項x「'。我刪除了方程式右邊的自然頻繁項。只是爲了讓我的程序看起來更簡單。讓我看看它是否有效。 – ADK

爲什麼你會確定使用隨機序列作爲位移的初始條件而不是速度是正確的？ – gg349

我的意思是theta2（在t = 0時）是theta在（t = 0）時的導數。那麼，如果theta2不依賴theta？不管怎麼說，還是要謝謝你。 – ADK

使用scipy二階ODE積分

回答

相關問題