Saya dapat mensimulasikan sistem pegas-massa di bawah osilasi teredam. Namun, saya ingin menambahkan subplot posisi vs waktu dan subplot lain kecepatan vs posisi (jalur fase) sehingga saya akan memiliki tiga animasi. Bagaimana saya bisa menambahkannya? Kode sumber yang saya gunakan ditunjukkan di bawah ini

Pembaruan: Saya mencoba menambahkan posisi subplot pertama vs waktu tetapi saya tidak bisa mendapatkan kurva yang diinginkan untuk itu.

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation


#Constants
w = np.pi #angular frequency
b = np.pi * 0.2  #damping parameter

#Function that implements rk4 integration
def rk4(t, x, f, dt):
    dx1 = f(t, x)*dt
    dx2 = f(t+0.5*dt, x+0.5*dx1)*dt
    dx3 = f(t+0.5*dt, x+0.5*dx2)*dt
    dx4 = f(t+dt, x+dx3)*dt
    return x+dx1/6.0+dx2/3.0+dx3/3.0+dx4/6.0

#Function that returns dX/dt for the linearly damped oscillator
def dXdt(t, X):
    x = X[0]
    vx = X[1]
    ax = -2*b*vx - w**2*x
    return np.array([vx, ax])

#Initialize Variables
x0 = 5.0 #Initial x position
vx0 = 1.0 #Initial x Velocity

#Setting time array for graph visualization
dt = 0.05
N = 200
x = np.zeros(N)
vx = np.zeros(N)
y = []
# integrate equations of motion using rk4;
# X is a vector that contains the positions and velocities being integrated
X = np.array([x0, vx0])


for i in range(N):
    x[i] = X[0]
    vx[i] = X[1]
    y.append(0)
    # update the vector X to the next time step
    X = rk4(i*dt, X, dXdt, dt)

fig, (ax1, ax2) = plt.subplots(2,1, figsize=(8,6))
fig.suptitle(r' Damped Oscillation with $\beta$$\approx$' + str(round(b,2)) + r' and $\omega$$\approx$'
             + str(round(w,2)), fontsize=16)
line1, = ax1.plot([], [], lw=10,c="blue",ls="-",ms=50,marker="s",mfc="gray",fillstyle="none",mec="black",markevery=2)
line2, = ax2.plot([], [], lw=2, color='r')
time_template = '\nTime = %.1fs'
time_text = ax1.text(0.1, 0.9, '', transform=ax1.transAxes)

for ax in [ax1, ax2]:
    ax1.set_xlim(1.2*min(x), 1.2*max(x))
    ax2.set_ylim(1.2*min(x), 1.2*max(x),1)
    ax2.set_xlim(0, N*dt)
    ax1.set_yticklabels([])

def init():
    line1.set_data([], [])
    line2.set_data([], [])
    time_text.set_text('')
    return line1, line2, time_text

def animate(i):
    thisx1 = [x[i],0]
    thisy1 = [y[i],0]
    thisx2 = [i*dt,0]
    thisy2 = [x[i],0]
    line1.set_data(thisx1, thisy1)
    line2.set_data(thisx2, thisy2)
    time_text.set_text(time_template % (i*dt))
    return line1, line2, time_text  

ani = animation.FuncAnimation(fig, animate, np.arange(1, N),
                              interval=50, blit=True, init_func=init,repeat=False)
2
Lawrence Calixto 12 Mei 2021, 09:02

1 menjawab

Jawaban Terbaik

Setelah beberapa anak di bawah umur mengubah kode awal Anda, yang paling penting adalah:

thisx1 = [x[i],0]
thisy1 = [y[i],0]
thisx2 = [i*dt,0]
thisy2 = [x[i],0]
line1.set_data(thisx1, thisy1)
line2.set_data(thisx2, thisy2)

# should be written like this
line1.set_data([x[i],0], [y[i],0])
line2.set_data(t[:i], x[:i])
line3.set_data(x[:i], vx[:i])

Versi kerja, dengan plot ruang fase berwarna hijau, adalah sebagai berikut:

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import animation


#Constants
w = np.pi #angular frequency
b = np.pi * 0.2  #damping parameter

#Function that implements rk4 integration
def rk4(t, x, f, dt):
    dx1 = f(t, x)*dt
    dx2 = f(t+0.5*dt, x+0.5*dx1)*dt
    dx3 = f(t+0.5*dt, x+0.5*dx2)*dt
    dx4 = f(t+dt, x+dx3)*dt
    return x+dx1/6.0+dx2/3.0+dx3/3.0+dx4/6.0

#Function that returns dX/dt for the linearly damped oscillator
def dXdt(t, X):
    x = X[0]
    vx = X[1]
    ax = -2*b*vx - w**2*x
    return np.array([vx, ax])

#Initialize Variables
x0 = 5.0 #Initial x position
vx0 = 1.0 #Initial x Velocity

#Setting time array for graph visualization
dt = 0.05
N = 200
t = np.linspace(0,N*dt,N,endpoint=False)
x = np.zeros(N)
vx = np.zeros(N)
y = np.zeros(N)
# integrate equations of motion using rk4;
# X is a vector that contains the positions and velocities being integrated
X = np.array([x0, vx0])

for i in range(N):
    x[i] = X[0]
    vx[i] = X[1]
    # update the vector X to the next time step
    X = rk4(i*dt, X, dXdt, dt)

fig, (ax1, ax2, ax3) = plt.subplots(3,1, figsize=(8,12))
fig.suptitle(r' Damped Oscillation with $\beta$$\approx$' + str(round(b,2)) + r' and $\omega$$\approx$'
             + str(round(w,2)), fontsize=16)
line1, = ax1.plot([], [], lw=10,c="blue",ls="-",ms=50,marker="s",mfc="gray",fillstyle="none",mec="black",markevery=2)
line2, = ax2.plot([], [], lw=1, color='r')
line3, = ax3.plot([], [], lw=1, color='g')
time_template = '\nTime = %.1fs'
time_text = ax1.text(0.1, 0.9, '', transform=ax1.transAxes)


ax1.set_xlim(1.2*min(x), 1.2*max(x))
ax1.set_xlabel('x')
ax1.set_ylabel('y')

ax2.set_ylim(1.2*min(x), 1.2*max(x),1)
ax2.set_xlim(0, N*dt)
ax2.set_xlabel('t')
ax2.set_ylabel('x')

ax3.set_xlim(1.2*min(x), 1.2*max(x),1)
ax3.set_ylim(1.2*min(vx), 1.2*max(vx),1)
ax3.set_xlabel('x')
ax3.set_ylabel('vx')

def init():
    line1.set_data([], [])
    line2.set_data([], [])
    line3.set_data([], [])
    time_text.set_text('')
    return line1, line2, line3, time_text

def animate(i):
    line1.set_data([x[i],0], [y[i],0])
    line2.set_data(t[:i], x[:i])
    line3.set_data(x[:i], vx[:i])
    time_text.set_text(time_template % (i*dt))
    return line1, line2, line3, time_text  

ani = animation.FuncAnimation(fig, animate, np.arange(1, N),
                              interval=50, blit=True, init_func=init,repeat=False)

ani.save('anim.gif')

Dan memberikan:

anim.gif

2
Yacola 12 Mei 2021, 14:12