python第八章課後習題

明天11234發表於2024-11-13
Python

8.4 求微分方程組的數值解 x'=-x3-y,x(0)=1,y'=x-y3,y(0)=0.5,0<=t<=30，要求畫出x（t）和y（t）的解曲線圖形，再相平面上畫出軌線

點選檢視程式碼 
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

def system(t, state):
    x, y = state
    dxdt = -x**3 - y
    dydt = x - y**3
    return [dxdt, dydt]

t_span = (0, 30)
y0 = [1, 0.5]

sol = solve_ivp(system, t_span, y0, t_eval=np.linspace(t_span[0], t_span[1], 1000))

x = sol.y[0]
y = sol.y[1]
t = sol.t

# 繪製 x(t) vs t 和 y(t) vs t
plt.figure(figsize=(12, 6))

plt.subplot(1, 2, 1)
plt.plot(t, x, label='x(t)')
plt.xlabel('t')
plt.ylabel('x(t)')
plt.title('x(t) vs t')
plt.legend()
plt.grid(True)

plt.subplot(1, 2, 2)
plt.plot(t, y, label='y(t)')
plt.xlabel('t')
plt.ylabel('y(t)')
plt.title('y(t) vs t')
plt.legend()
plt.grid(True)

plt.figure(figsize=(6, 6))
plt.plot(x, y, label='Phase Trajectory')
plt.xlabel('x')
plt.ylabel('y')
plt.title('Phase Plane (x, y)')
plt.legend()
plt.grid(True)
plt.axis('equal')  

plt.tight_layout()  
plt.show()


print("學號：2023310143028")

8.5 求微分方程組(豎直加熱板的自然對流)的數值解。(d3f)/(dm3)+3f(df)/(dm2)-2((df)/(dm))2+T=0,(d2T)/(dm2)+2.1f(dT)/(dm)=0，已知當m=0時,f=0,(df)/(dm)=0,(d2f)/(dm^2)=0.68,T=1,(dT)/(dm)=-0.5。要求在[0,10]區間上，畫出f(a)、T(a)的解曲線

點選檢視程式碼 
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

def model(t, y):
    f, df_dm, d2f_dm2, T, dT_dm = y
    d3f_dm3 = -3 * f * d2f_dm2 + 2 * (df_dm)**2 - T  
    d2T_dm2 = -2.1 * f * dT_dm  # 注意透過 y 引用 f
   
    return [df_dm, d2f_dm2, d3f_dm3, dT_dm, d2T_dm2]

y0 = [0, 0, 0.68, 1, -0.5]

t_span = (0, 10)
t_eval = np.linspace(t_span[0], t_span[1], 1000)


sol = solve_ivp(model, t_span, y0, t_eval=t_eval, method='RK45')

f = sol.y[0]
T = sol.y[3]

plt.figure(figsize=(12, 6))

plt.subplot(2, 1, 1)
plt.plot(sol.t, f, label='f(t)')  
plt.xlabel('t') 
plt.ylabel('f(t)')
plt.title('Solution for f(t)')
plt.grid(True)

plt.subplot(2, 1, 2)
plt.plot(sol.t, T, label='T(t)', color='orange')  
plt.xlabel('t')  
plt.ylabel('T(t)')
plt.title('Solution for T(t)')
plt.grid(True)

plt.tight_layout()
plt.show()


print("學號：2023310143028")

相關文章