第八章習題

倒头睡不醒a發表於2024-11-18

學號後四位:3018
8.4:

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


# 定義微分方程組
def differential_equations(state, t):
    x, y = state
    dxdt = -x ** 3 - y
    dydt = x - y ** 3
    return [dxdt, dydt]


# 設定初始條件
initial_conditions = [1, 0.5]

# 定義時間範圍
t = np.linspace(0, 30, 1000)

# 求解微分方程組
solution = odeint(differential_equations, initial_conditions, t)
x_solution = solution[:, 0]
y_solution = solution[:, 1]

# 繪製x(t)的解曲線
plt.subplot(2, 1, 1)
plt.plot(t, x_solution)
plt.xlabel('t')
plt.ylabel('x(t)')
plt.title('Solution Curve of x(t)')

# 繪製y(t)的解曲線
plt.subplot(2, 1, 2)
plt.plot(t, y_solution)
plt.xlabel('t')
plt.ylabel('y(t)')
plt.title('Solution Curve of y(t)')

# 在相平面上繪製軌線
plt.figure()
plt.plot(x_solution, y_solution)
plt.xlabel('x')
plt.ylabel('y')
plt.title('Phase Plane Trajectory')

plt.show()
print("xuehao3018")

8.5:

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

# 定義微分方程組
def equations(state, eta):
    f, df_deta, T, dT_deta = state
    # 用中心差分近似計算二階導數
    h = 1e-4
    df2_deta2_approx = (state[1] - odeint(equations, [state[0]+h, df_deta+h, state[2], dT_deta], [eta])[0][1]) / h**2
    d2T_deta2 = -2.1 * f * dT_deta
    return [df_deta, df2_deta2_approx, dT_deta, d2T_deta2]

# 初始條件
f0 = 0
df_deta0 = -0.5
T0 = 1
dT_deta0 = 0
initial_conditions = [f0, df_deta0, T0, dT_deta0]

# 定義 eta 的範圍
eta = np.linspace(0, 10, 1000)

# 求解微分方程組
solution = odeint(equations, initial_conditions, eta)
f_solution = solution[:, 0]
T_solution = solution[:, 2]

# 繪製 f(η) 的解曲線
plt.plot(eta, f_solution, label='f(η)')
plt.xlabel('η')
plt.ylabel('f')
plt.title('Solution Curve of f(η)')

# 繪製 T(η) 的解曲線
plt.plot(eta, T_solution, label='T(η)')
plt.xlabel('η')
plt.ylabel('T')
plt.title('Solution Curve of T(η)')

plt.legend()
plt.show()
print("xuehao3018")

相關文章