2016-05-02 9 views
1

私はこれらのソリューションの初期条件を変更するだけで、同じグラフに2つのソリューションをプロットしようとしています。私は、初期条件をy = 0とy = -pi/2の1つの値にします。私はそれらを同じグラフ上にどのようにして、すべてをコピーして貼り付けたり、新しい配列を作成したりすることなく、どうしたらいいのか分かりません。1つのグラフで同じドライブ強度を持つ2つのソリューションをプロットする

import numpy as np 
import matplotlib.pyplot as plt 

t = 0.0 
y = 0.0 
u = 0.0 
F = 1.073 
Wd = 2*np.pi 
w0 = 1.5*Wd 
b = w0/4 

ts =[] 
ys =[] 

h= 0.001 

while (t <= 30): 
    m1 = u 
    k1 = (-w0**2)*np.sin(y) + u*(-2*b) + F*(w0**2)*np.cos(Wd*t) 
    m2 = u + (h/2.) * k1 
    t_2 = t + (h/2.) 
    y_2 = y +(h/2.) * m1 
    u_2 = m2 
    k2 = (-w0**2)*np.sin(y_2) + u_2*(-2*b) + F*(w0**2)*np.cos(Wd*t_2) 
    m3 = u + (h/2.) * k2 
    t_3 = t + (h/2.) 
    y_3 = y + (h/2.) * m2 
    u_3 = m3 
    k3 = (-w0**2)*np.sin(y_3) + u_3*(-2*b) + F*(w0**2)*np.cos(Wd*t_3) 
    m4 = u + h * k3 
    t_4 = t + h 
    y_4 = y + h * m3 
    u_4 = m4 
    k4 = (-w0**2)*np.sin(y_4) + u_4*(-2*b) + F*(w0**2)*np.cos(Wd*t_4) 
    t = t + h 
    y = y + (h/6.) * (m1 + (2 * m2) + (2 * m3) + m4) 
    u = u + (h/6.) * (k1 + (2 * k2) + (2 * k3) + k4) 
    ts.append(t) 
    ys.append(y) 

plt.plot(ts,ys) 
plt.xlabel('$t$',fontsize=18) 
plt.ylabel('$\phi$',fontsize=18) 
plt.axhline(y=0, color = 'black') 
plt.show() 

答えて

1
import numpy as np 
import matplotlib.pyplot as plt 


def plotit(y): 
    t = 0.0 
    #y = 0.0 
    u = 0.0 
    F = 1.073 
    Wd = 2*np.pi 
    w0 = 1.5*Wd 
    b = w0/4 

    ts =[] 
    ys =[] 

    h= 0.001 

    while (t <= 30): 
     m1 = u 
     k1 = (-w0**2)*np.sin(y) + u*(-2*b) + F*(w0**2)*np.cos(Wd*t) 
     m2 = u + (h/2.) * k1 
     t_2 = t + (h/2.) 
     y_2 = y +(h/2.) * m1 
     u_2 = m2 
     k2 = (-w0**2)*np.sin(y_2) + u_2*(-2*b) + F*(w0**2)*np.cos(Wd*t_2) 
     m3 = u + (h/2.) * k2 
     t_3 = t + (h/2.) 
     y_3 = y + (h/2.) * m2 
     u_3 = m3 
     k3 = (-w0**2)*np.sin(y_3) + u_3*(-2*b) + F*(w0**2)*np.cos(Wd*t_3) 
     m4 = u + h * k3 
     t_4 = t + h 
     y_4 = y + h * m3 
     u_4 = m4 
     k4 = (-w0**2)*np.sin(y_4) + u_4*(-2*b) + F*(w0**2)*np.cos(Wd*t_4) 
     t = t + h 
     y = y + (h/6.) * (m1 + (2 * m2) + (2 * m3) + m4) 
     u = u + (h/6.) * (k1 + (2 * k2) + (2 * k3) + k4) 
     ts.append(t) 
     ys.append(y) 

    plt.plot(ts,ys) 
    plt.xlabel('$t$',fontsize=18) 
    plt.ylabel('$\phi$',fontsize=18) 
    plt.axhline(y=0, color = 'black') 



plotit(y=0) 
plotit(y=-np.pi/2.0) 
plt.show() 
+0

は完璧に動作します! – Cosmoman

関連する問題