第八章

• 8.4代码

点击查看代码

import numpy as np
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

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("学号后两位：15")

• 代码结果
  

• 8.5

点击查看代码

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
# 定义微分方程组
def system(state, eta):
    f, df, d2f, T, dT = state
    d3f = -3 * d2f + 2 * df ** 2 - T
    d2T = -2.1 * df * dT
    return [df, d2f, d3f, dT, d2T]
# 初始条件
f0 = 0
df0 = 0
d2f0 = 0.68
T0 = 1
dT0 = -0.5
state0 = [f0, df0, d2f0, T0, dT0]

# 时间范围
eta = np.linspace(0, 10, 1000)

# 求解微分方程组
solution = odeint(system, state0, eta)
f_sol = solution[:, 0]
T_sol = solution[:, 3]

# 绘制解曲线
plt.figure()
plt.plot(eta, f_sol, label='f(eta)')
plt.plot(eta, T_sol, label='T(eta)')
plt.xlabel('eta')
plt.ylabel('f, T')
plt.title('Solution of the system of differential equations')
plt.legend()
plt.show()

• 8.8

点击查看代码

#P8.8
 
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
 
# 设置matplotlib以使用支持中文的字体
matplotlib.rcParams['font.sans-serif'] = ['SimHei']  # 使用黑体
matplotlib.rcParams['axes.unicode_minus'] = False  # 正确显示负号
 
# 参数设置
a = 1 - 0.2 ** (1 / 12)
m = 1.109e5
w3 = 17.86
w4 = 22.99
 
# 初始化变量
X = []
Z = []
K = np.arange(0, 0.875, 0.001)
 
# 循环计算
for k in K:
    # 计算x1, x2, x3, x4
    x1 = 1.22e11 * (1 - 1 / (m * (1 - a - 0.42 * k))) ** 8 * (1 - a) ** 24 * (0.5 + (1 - a - k) ** 8 * (1 - a) ** 4 / (1 - (1 - a - k) ** 8 * (1 - a) ** 4))
    x2 = (1 - a) ** 12 * x1
    x3 = (1 - a) ** 12 * x2
    x4 = (1 - a - 0.42 * k) ** 8 * (1 - a) ** 4 / (1 - (1 - a - k) ** 8 * (1 - a) ** 4) * x3
    
    # 将结果添加到X列表中，注意要转换为列向量形式
    X.append([x1, x2, x3, x4])
    
    # 计算n（虽然未在MATLAB代码中直接使用，但为保持一致性，这里也计算出来）
    n = m * (1 - a - 0.42 * k) ** 8 * (1 - a) ** 24 * (0.5 + (1 - a - k) ** 8 * (1 - a) ** 4 / (1 - (1 - a - k) ** 8 * (1 - a) ** 4)) * x1
    
    # 计算z
    z = 0.42 * k * w3 * (1 - (1 - a - 0.42 * k) ** 8) / (a + 0.42 * k) * x3 + \
        k * w4 * (1 - (1 - a - k) ** 8) / (a + k) * x4
    
    # 将结果添加到Z列表中
    Z.append(z)
 
# 将X和Z转换为NumPy数组，方便后续处理
X = np.array(X).T  # 转置为列向量形式
Z = np.array(Z)
 
# 查找最大值
mz = np.max(Z)
ind = np.argmax(Z)
k4 = K[ind]
k3 = 0.42 * k4  # 最优捕捞强度
 
# 最优捕捞强度下的各年龄组鱼群数
xx = X[:, ind]
 
# 绘图
plt.plot(K, Z)
plt.xlabel('$k$', fontsize=14)
plt.ylabel('$Z$', fontsize=14, rotation=0)
plt.title('捕捞强度与总收益的关系')
 
plt.grid(True)
plt.show()
print("15")

• 代码结果