"""生成所有实验报告所需的图片""" import numpy as np from scipy import interpolate, integrate, optimize from scipy.stats import norm import matplotlib matplotlib.use('Agg') import matplotlib.pyplot as plt OUT = '/Users/chenhuidediannao/Desktop/数学建模实验报告' plt.rcParams['font.sans-serif'] = ['Arial Unicode MS', 'SimHei', 'DejaVu Sans'] plt.rcParams['axes.unicode_minus'] = False # ==================== 实验一 ==================== print("生成实验一图片...") # 1a 插值对比 x = np.array([0,1,2,3,4,5,6,7,8,9,10]) y = np.sin(x) x_dense = np.linspace(0, 10, 200) y_true = np.sin(x_dense) f_linear = interpolate.interp1d(x, y, kind='linear') y_linear = f_linear(x_dense) tck2 = interpolate.splrep(x, y, k=2) y_bs2 = interpolate.splev(x_dense, tck2) tck3 = interpolate.splrep(x, y, k=3) y_bs3 = interpolate.splev(x_dense, tck3) fig, ax = plt.subplots(figsize=(10, 6)) ax.plot(x, y, 'ko', markersize=8, label='Original data points') ax.plot(x_dense, y_true, 'k-', alpha=0.4, label='True sin(x)') ax.plot(x_dense, y_linear, 'r--', label='Linear interpolation') ax.plot(x_dense, y_bs2, 'g-.', label='Quadratic B-spline') ax.plot(x_dense, y_bs3, 'b-', label='Cubic B-spline') ax.legend(); ax.set_title('Interpolation Methods Comparison') ax.set_xlabel('x'); ax.set_ylabel('y'); ax.grid(True, alpha=0.3) fig.tight_layout() fig.savefig(f'{OUT}/exp1_interp1.png', dpi=150, bbox_inches='tight') plt.close(fig) # 1b 二元插值 np.random.seed(42) xr = np.clip(np.random.normal(0, 0.4, 100), -1, 1) yr = np.clip(np.random.normal(0, 0.4, 100), -1, 1) zr = (xr + yr) * np.exp(-5*(xr**2 + yr**2)) xi = np.linspace(-1, 1, 100) yi = np.linspace(-1, 1, 100) XI, YI = np.meshgrid(xi, yi) ZI = interpolate.griddata((xr, yr), zr, (XI, YI), method='cubic') ZI_true = (XI + YI) * np.exp(-5*(XI**2 + YI**2)) fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5)) c1 = ax1.contourf(XI, YI, ZI, levels=20, cmap='viridis') ax1.scatter(xr, yr, c='red', s=5, alpha=0.5) ax1.set_title('Interpolation Result (100 points)') fig.colorbar(c1, ax=ax1) c2 = ax2.contourf(XI, YI, ZI_true, levels=20, cmap='viridis') ax2.set_title('True Function') fig.colorbar(c2, ax=ax2) fig.tight_layout() fig.savefig(f'{OUT}/exp1_interp2.png', dpi=150, bbox_inches='tight') plt.close(fig) # 1c cos/arccos x1 = np.linspace(0, np.pi, 200) x2 = np.linspace(-1, 1, 200) x3 = np.linspace(0, np.pi, 200) fig, ax = plt.subplots(figsize=(8, 8)) ax.plot(x1, np.cos(x1), 'b-', lw=2, label='y = cos(x)') ax.plot(x2, np.arccos(x2), 'r--', lw=2, label='y = arccos(x)') ax.plot(x3, x3, 'k:', lw=1, label='y = x') ax.axhline(y=0, color='gray', alpha=0.3); ax.axvline(x=0, color='gray', alpha=0.3) ax.set_xlim([-0.2, np.pi+0.2]); ax.set_ylim([-0.2, np.pi+0.2]) ax.grid(True, alpha=0.3); ax.legend(fontsize=12) ax.set_title('cos(x) and arccos(x) are symmetric about y=x') ax.set_aspect('equal') fig.tight_layout() fig.savefig(f'{OUT}/exp1_graph1.png', dpi=150, bbox_inches='tight') plt.close(fig) # 1d 椭圆参数方程 t = np.linspace(0, 2*np.pi, 500) fig, ax = plt.subplots(figsize=(10, 6)) ax.plot(2*np.cos(t), np.sin(t), 'b-', lw=2) ax.axhline(y=0, color='gray', alpha=0.3); ax.axvline(x=0, color='gray', alpha=0.3) ax.grid(True, alpha=0.3); ax.axis('equal') ax.set_title('Ellipse: x=2cos(t), y=sin(t)') ax.set_xlabel('x'); ax.set_ylabel('y') fig.tight_layout() fig.savefig(f'{OUT}/exp1_graph2.png', dpi=150, bbox_inches='tight') plt.close(fig) # ==================== 实验二 ==================== print("生成实验二图片...") # 2a 包装成本模型 w1, P1 = 50, 1.50 w2, P2 = 120, 3.00 A = np.array([[w1, w1**(2/3)], [w2, w2**(2/3)]]) b_arr = np.array([P1, P2]) c, k = np.linalg.solve(A, b_arr) w = np.linspace(10, 300, 200) P = c*w + k*w**(2/3) p = P / w fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5)) ax1.plot(w, P, 'b-', lw=2) ax1.plot([w1, w2], [P1, P2], 'ro', markersize=8) ax1.set_xlabel('Weight w (g)'); ax1.set_ylabel('Total Price P (yuan)') ax1.set_title('Total Price vs Weight') ax1.grid(True, alpha=0.3) ax2.plot(w, p, 'r-', lw=2) ax2.set_xlabel('Weight w (g)'); ax2.set_ylabel('Unit Price p (yuan/g)') ax2.set_title('Unit Price vs Weight (decreasing)') ax2.grid(True, alpha=0.3) fig.tight_layout() fig.savefig(f'{OUT}/exp2_pack.png', dpi=150, bbox_inches='tight') plt.close(fig) # 2b 雨中行走 rho, S, D, a_p, b_p = 1.0, 1.0, 100.0, 0.5, 0.3 def rainfall(v, ux, uy, uz): T = D / v return S*abs(v-ux)*T + a_p*S*abs(uy)*T + b_p*S*abs(uz)*T v_range = np.linspace(1, 15, 200) scenarios = [('No wind (ux=0)', 0), ('Tailwind (ux=-2)', -2), ('Headwind (ux=4)', 4)] fig, axes = plt.subplots(1, 3, figsize=(14, 4)) for ax, (title, ux) in zip(axes, scenarios): Q = [rainfall(v, ux, 3, 5) for v in v_range] ax.plot(v_range, Q, 'b-', lw=2) ax.set_xlabel('Speed v (m/s)'); ax.set_ylabel('Total Rain Q') ax.set_title(title); ax.grid(True, alpha=0.3) idx = np.argmin(Q) ax.plot(v_range[idx], Q[idx], 'ro', markersize=6) fig.tight_layout() fig.savefig(f'{OUT}/exp2_rain.png', dpi=150, bbox_inches='tight') plt.close(fig) # 2c 举重模型 weight = np.array([54, 59, 64, 70, 76, 83, 91, 99, 108, 110]) total = np.array([287.5, 307.5, 335, 357.5, 367.5, 392.5, 402.5, 420, 430, 457.5]) k1 = np.mean(total / weight**(2/3)) def power_law(w, k, alpha): return k * w**alpha popt, _ = optimize.curve_fit(power_law, weight, total, p0=[10, 0.67]) k2, alpha = popt w_smooth = np.linspace(50, 115, 200) r2_1 = 1 - np.sum((total - k1*weight**(2/3))**2) / np.sum((total - np.mean(total))**2) r2_2 = 1 - np.sum((total - power_law(weight, k2, alpha))**2) / np.sum((total - np.mean(total))**2) fig, ax = plt.subplots(figsize=(10, 6)) ax.scatter(weight, total, c='red', s=80, zorder=5, label='Olympic data') ax.plot(w_smooth, k1*w_smooth**(2/3), 'b--', lw=2, label=f'Model 1: L={k1:.2f}w^(2/3), R²={r2_1:.4f}') ax.plot(w_smooth, power_law(w_smooth, k2, alpha), 'g-', lw=2, label=f'Model 2: L={k2:.2f}w^{alpha:.2f}, R²={r2_2:.4f}') ax.set_xlabel('Body Weight (kg)'); ax.set_ylabel('Total Score (kg)') ax.set_title('Weightlifting Performance vs Body Weight') ax.legend(); ax.grid(True, alpha=0.3) fig.tight_layout() fig.savefig(f'{OUT}/exp2_lift.png', dpi=150, bbox_inches='tight') plt.close(fig) # ==================== 实验三 ==================== print("生成实验三图片...") # 3a 仙鹤种群 N0, years = 100, 50 t = np.arange(0, years+1) N_good = N0 * np.exp(0.0194*t) N_medium = N0 * np.exp(-0.0324*t) N_bad = N0 * np.exp(-0.0382*t) Na = np.zeros(years+1); Na[0] = N0 for i in range(years): Na[i+1] = Na[i] * np.exp(-0.0324) + 5 fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5)) ax1.plot(t, N_good, 'g-', lw=2, label='Good (r=0.0194)') ax1.plot(t, N_medium, 'orange', lw=2, label='Medium (r=-0.0324)') ax1.plot(t, N_bad, 'r-', lw=2, label='Poor (r=-0.0382)') ax1.axhline(y=0, color='gray', alpha=0.3) ax1.set_xlabel('Year'); ax1.set_ylabel('Crane Population') ax1.set_title('Crane Population in Different Environments (No intervention)') ax1.legend(); ax1.grid(True, alpha=0.3) ax2.plot(t, N_medium, 'orange', lw=2, label='Medium (no aid)') ax2.plot(t, Na, 'b-', lw=2, label='Medium (5 released/year)') ax2.axhline(y=0, color='gray', alpha=0.3) ax2.set_xlabel('Year'); ax2.set_ylabel('Crane Population') ax2.set_title('Effect of Artificial Release in Medium Environment') ax2.legend(); ax2.grid(True, alpha=0.3) fig.tight_layout() fig.savefig(f'{OUT}/exp3_crane.png', dpi=150, bbox_inches='tight') plt.close(fig) # 3b 火箭发射 g, k_d, m0, m_fuel = 9.8, 0.4, 1400.0, 1080.0 m_dry = m0 - m_fuel burn_rate, thrust, t_burn = 18.0, 32000.0, m_fuel / 18.0 def rocket_ode(t, y): h, v = y m = m0 - burn_rate * t drag = k_d * v * abs(v) a = (thrust - m * g - drag) / m return [v, a] sol1 = integrate.solve_ivp(rocket_ode, [0, t_burn], [0, 0], max_step=0.1, dense_output=True) def coast_ode(t, y): h, v = y drag = k_d * v * abs(v) a = (-m_dry * g - drag) / m_dry return [v, a] def find_apex(t, y): return y[1] find_apex.terminal = True; find_apex.direction = -1 h_bo, v_bo = sol1.y[0, -1], sol1.y[1, -1] sol2 = integrate.solve_ivp(coast_ode, [t_burn, 500], [h_bo, v_bo], max_step=0.1, events=find_apex, dense_output=True) t_apex = sol2.t[-1]; h_apex = sol2.y[0, -1] t1 = np.linspace(0, t_burn, 200); y1 = sol1.sol(t1) # y1[0]=h, y1[1]=v # Compute acceleration from ODE def rocket_accel(t, h, v): m = m0 - burn_rate * t if t <= t_burn else m_dry drag = k_d * v * abs(v) if t <= t_burn: return (thrust - m * g - drag) / m else: return (-m_dry * g - drag) / m_dry a1 = np.array([rocket_accel(ti, hi, vi) for ti, hi, vi in zip(t1, y1[0], y1[1])]) t2 = np.linspace(t_burn, t_apex, 200); y2 = sol2.sol(t2) a2 = np.array([rocket_accel(ti, hi, vi) for ti, hi, vi in zip(t2, y2[0], y2[1])]) fig, axes = plt.subplots(3, 1, figsize=(12, 10)) data_pairs = [(y1[0], y2[0], 'h (m)', 'Altitude'), (y1[1], y2[1], 'v (m/s)', 'Velocity'), (a1, a2, 'a (m/s²)', 'Acceleration')] for ax, (d1, d2, ylabel, title) in zip(axes, data_pairs): ax.plot(t1, d1, 'b-', lw=1.5, label='Powered flight') ax.plot(t2, d2, 'r-', lw=1.5, label='Coasting') ax.axvline(x=t_burn, color='gray', linestyle='--', alpha=0.5, label=f'Burnout t={t_burn:.0f}s') ax.set_ylabel(ylabel); ax.set_title(title) ax.legend(); ax.grid(True, alpha=0.3) axes[2].set_xlabel('Time t (s)') fig.tight_layout() fig.savefig(f'{OUT}/exp3_rocket.png', dpi=150, bbox_inches='tight') plt.close(fig) a_bo = rocket_accel(t_burn, h_bo, v_bo) print(f" 引擎关闭: h={h_bo:.1f}m, v={v_bo:.1f}m/s, a={a_bo:.1f}m/s²") print(f" 最高点: h={h_apex:.1f}m, t={t_apex:.1f}s") # ==================== 实验四 ==================== print("生成实验四图片...") x = np.array([0, 8.20, 0.50, 5.70, 0.77, 2.87, 4.43, 2.58, 0.72, 9.76, 3.19, 5.55]) y = np.array([0, 0.50, 4.0, 5.0, 6.49, 8.76, 3.26, 9.32, 9.96, 3.16, 7.20, 7.88]) R = np.array([600, 1000, 800, 1400, 1200, 700, 600, 800, 1000, 1200, 1000, 1100]) def objective(pos): cx, cy = pos return np.sum(R * np.sqrt((x-cx)**2 + (y-cy)**2)) x0_w = np.sum(R*x)/np.sum(R); y0_w = np.sum(R*y)/np.sum(R) best_val, best_pos = float('inf'), None for start in [(x0_w, y0_w), (5,5), (3,6), (4,4), (2,7)]: res = optimize.minimize(objective, start, method='Nelder-Mead', options={'xatol': 1e-8, 'fatol': 1e-8}) if res.fun < best_val: best_val = res.fun; best_pos = res.x fig, ax = plt.subplots(figsize=(10, 9)) sc = ax.scatter(x, y, s=R/5, c=R, cmap='YlOrRd', edgecolors='black', alpha=0.8, zorder=5) for i, (xi, yi, ri) in enumerate(zip(x, y, R)): ax.annotate(f'{i+1}({ri})', (xi, yi), xytext=(5, 5), textcoords='offset points', fontsize=8) ax.scatter(best_pos[0], best_pos[1], c='blue', marker='*', s=400, zorder=10, edgecolors='black', lw=2, label='Optimal Location') ax.scatter(x0_w, y0_w, c='green', marker='^', s=150, zorder=10, label='Weighted Mean') ax.set_xlabel('x (km)'); ax.set_ylabel('y (km)') ax.set_title('Island Settlement Distribution & Optimal Service Center Location') ax.legend(); ax.grid(True, alpha=0.3) fig.colorbar(sc, ax=ax, label='Population R') ax.axis('equal') fig.tight_layout() fig.savefig(f'{OUT}/exp4_location.png', dpi=150, bbox_inches='tight') plt.close(fig) print(f" 最优位置: ({best_pos[0]:.4f}, {best_pos[1]:.4f})") # ==================== 实验五 ==================== print("生成实验五图片...") # 5a 正态分布 mu, sigma = 20, 1.5 n = 25; sigma_xbar = sigma/np.sqrt(n) p1 = norm.cdf(22, mu, sigma) - norm.cdf(19, mu, sigma) x0 = norm.ppf(0.10, mu, sigma) p3 = norm.cdf(22, mu, sigma_xbar) - norm.cdf(19, mu, sigma_xbar) fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5)) xs = np.linspace(mu-4*sigma, mu+4*sigma, 300) ax1.plot(xs, norm.pdf(xs, mu, sigma), 'b-', lw=2, label=f'Single part N(20,1.5²)') ax1.fill_between(xs, norm.pdf(xs, mu, sigma), where=(xs>=19)&(xs<=22), alpha=0.3, color='blue') ax1.axvline(x=x0, color='red', linestyle='--', label=f'x0={x0:.2f} (10% quantile)') ax1.set_xlabel('Size (cm)'); ax1.set_ylabel('Probability Density') ax1.set_title(f'Single Part: P(19≤X≤22)={p1*100:.1f}%') ax1.legend() xb = np.linspace(mu-4*sigma_xbar, mu+4*sigma_xbar, 300) ax2.plot(xb, norm.pdf(xb, mu, sigma_xbar), 'g-', lw=2, label=f'Sample mean N(20,{sigma_xbar:.2f}²)') ax2.fill_between(xb, norm.pdf(xb, mu, sigma_xbar), where=(xb>=19)&(xb<=22), alpha=0.3, color='green') ax2.set_xlabel('Sample Mean (cm)'); ax2.set_ylabel('Probability Density') ax2.set_title(f'25 Samples Mean: P(19≤X̄≤22)≈{p3*100:.1f}%') ax2.legend() fig.tight_layout() fig.savefig(f'{OUT}/exp5_normal.png', dpi=150, bbox_inches='tight') plt.close(fig) print(f" P(19≤X≤22)={p1*100:.2f}%, x0={x0:.4f}, P(19≤X̄≤22)={p3*100:.2f}%") # 5b 报童模型 a_n, b_n, c_n = 0.8, 1.0, 0.75 demand_bins = np.array([100, 120, 140, 160, 180, 200, 220, 240, 260, 280, 300]) days = np.array([3, 9, 13, 22, 32, 35, 20, 15, 8, 2]) probs = days / np.sum(days) demand_mid = (demand_bins[:-1] + demand_bins[1:]) / 2 cdf = np.cumsum(probs) cr = (b_n - a_n) / (b_n - c_n) opt_idx = np.argmax(cdf >= cr) Q_opt = demand_mid[opt_idx] def expected_profit(Q): profit = 0 for d, p in zip(demand_mid, probs): sales = min(Q, d) leftover = max(0, Q - d) profit += p * (b_n * sales + c_n * leftover - a_n * Q) return profit Q_range = np.arange(100, 300, 5) profits = [expected_profit(q) for q in Q_range] Q_best = Q_range[np.argmax(profits)] fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5)) ax1.bar(demand_mid, probs, width=15, alpha=0.7, edgecolor='black', label='Empirical dist.') ax1.axvline(x=Q_opt, color='red', linestyle='--', lw=2, label=f'Q*≈{Q_opt:.0f}') ax1.set_xlabel('Demand'); ax1.set_ylabel('Probability') ax1.set_title('Demand Distribution & Optimal Order Quantity'); ax1.legend() ax2.plot(Q_range, profits, 'b-', lw=2) ax2.axvline(x=Q_best, color='red', linestyle='--', label=f'Q*={Q_best}') ax2.set_xlabel('Order Qty Q'); ax2.set_ylabel('Expected Profit (yuan)') ax2.set_title('Expected Profit vs Order Quantity'); ax2.legend(); ax2.grid(True, alpha=0.3) fig.tight_layout() fig.savefig(f'{OUT}/exp5_newsboy.png', dpi=150, bbox_inches='tight') plt.close(fig) print(f" 最优订货量: Q≈{Q_opt:.0f}, 期望利润: {max(profits):.2f}") print("\n所有图片生成完毕！")
