生成所有实验报告所需的图片 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, kindlinear) y_linear f_linear(x_dense) tck2 interpolate.splrep(x, y, k2) y_bs2 interpolate.splev(x_dense, tck2) tck3 interpolate.splrep(x, y, k3) y_bs3 interpolate.splev(x_dense, tck3) fig, ax plt.subplots(figsize(10, 6)) ax.plot(x, y, ko, markersize8, labelOriginal data points) ax.plot(x_dense, y_true, k-, alpha0.4, labelTrue sin(x)) ax.plot(x_dense, y_linear, r--, labelLinear interpolation) ax.plot(x_dense, y_bs2, g-., labelQuadratic B-spline) ax.plot(x_dense, y_bs3, b-, labelCubic B-spline) ax.legend(); ax.set_title(Interpolation Methods Comparison) ax.set_xlabel(x); ax.set_ylabel(y); ax.grid(True, alpha0.3) fig.tight_layout() fig.savefig(f{OUT}/exp1_interp1.png, dpi150, bbox_inchestight) 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), methodcubic) 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, levels20, cmapviridis) ax1.scatter(xr, yr, cred, s5, alpha0.5) ax1.set_title(Interpolation Result (100 points)) fig.colorbar(c1, axax1) c2 ax2.contourf(XI, YI, ZI_true, levels20, cmapviridis) ax2.set_title(True Function) fig.colorbar(c2, axax2) fig.tight_layout() fig.savefig(f{OUT}/exp1_interp2.png, dpi150, bbox_inchestight) 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-, lw2, labely cos(x)) ax.plot(x2, np.arccos(x2), r--, lw2, labely arccos(x)) ax.plot(x3, x3, k:, lw1, labely x) ax.axhline(y0, colorgray, alpha0.3); ax.axvline(x0, colorgray, alpha0.3) ax.set_xlim([-0.2, np.pi0.2]); ax.set_ylim([-0.2, np.pi0.2]) ax.grid(True, alpha0.3); ax.legend(fontsize12) ax.set_title(cos(x) and arccos(x) are symmetric about yx) ax.set_aspect(equal) fig.tight_layout() fig.savefig(f{OUT}/exp1_graph1.png, dpi150, bbox_inchestight) 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-, lw2) ax.axhline(y0, colorgray, alpha0.3); ax.axvline(x0, colorgray, alpha0.3) ax.grid(True, alpha0.3); ax.axis(equal) ax.set_title(Ellipse: x2cos(t), ysin(t)) ax.set_xlabel(x); ax.set_ylabel(y) fig.tight_layout() fig.savefig(f{OUT}/exp1_graph2.png, dpi150, bbox_inchestight) 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-, lw2) ax1.plot([w1, w2], [P1, P2], ro, markersize8) ax1.set_xlabel(Weight w (g)); ax1.set_ylabel(Total Price P (yuan)) ax1.set_title(Total Price vs Weight) ax1.grid(True, alpha0.3) ax2.plot(w, p, r-, lw2) 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, alpha0.3) fig.tight_layout() fig.savefig(f{OUT}/exp2_pack.png, dpi150, bbox_inchestight) 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 (ux0), 0), (Tailwind (ux-2), -2), (Headwind (ux4), 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-, lw2) ax.set_xlabel(Speed v (m/s)); ax.set_ylabel(Total Rain Q) ax.set_title(title); ax.grid(True, alpha0.3) idx np.argmin(Q) ax.plot(v_range[idx], Q[idx], ro, markersize6) fig.tight_layout() fig.savefig(f{OUT}/exp2_rain.png, dpi150, bbox_inchestight) 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, cred, s80, zorder5, labelOlympic data) ax.plot(w_smooth, k1*w_smooth**(2/3), b--, lw2, labelfModel 1: L{k1:.2f}w^(2/3), R²{r2_1:.4f}) ax.plot(w_smooth, power_law(w_smooth, k2, alpha), g-, lw2, labelfModel 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, alpha0.3) fig.tight_layout() fig.savefig(f{OUT}/exp2_lift.png, dpi150, bbox_inchestight) plt.close(fig) # 实验三 print(生成实验三图片...) # 3a 仙鹤种群 N0, years 100, 50 t np.arange(0, years1) 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(years1); Na[0] N0 for i in range(years): Na[i1] Na[i] * np.exp(-0.0324) 5 fig, (ax1, ax2) plt.subplots(1, 2, figsize(14, 5)) ax1.plot(t, N_good, g-, lw2, labelGood (r0.0194)) ax1.plot(t, N_medium, orange, lw2, labelMedium (r-0.0324)) ax1.plot(t, N_bad, r-, lw2, labelPoor (r-0.0382)) ax1.axhline(y0, colorgray, alpha0.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, alpha0.3) ax2.plot(t, N_medium, orange, lw2, labelMedium (no aid)) ax2.plot(t, Na, b-, lw2, labelMedium (5 released/year)) ax2.axhline(y0, colorgray, alpha0.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, alpha0.3) fig.tight_layout() fig.savefig(f{OUT}/exp3_crane.png, dpi150, bbox_inchestight) 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_step0.1, dense_outputTrue) 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_step0.1, eventsfind_apex, dense_outputTrue) 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-, lw1.5, labelPowered flight) ax.plot(t2, d2, r-, lw1.5, labelCoasting) ax.axvline(xt_burn, colorgray, linestyle--, alpha0.5, labelfBurnout t{t_burn:.0f}s) ax.set_ylabel(ylabel); ax.set_title(title) ax.legend(); ax.grid(True, alpha0.3) axes[2].set_xlabel(Time t (s)) fig.tight_layout() fig.savefig(f{OUT}/exp3_rocket.png, dpi150, bbox_inchestight) 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, methodNelder-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, sR/5, cR, cmapYlOrRd, edgecolorsblack, alpha0.8, zorder5) for i, (xi, yi, ri) in enumerate(zip(x, y, R)): ax.annotate(f{i1}({ri}), (xi, yi), xytext(5, 5), textcoordsoffset points, fontsize8) ax.scatter(best_pos[0], best_pos[1], cblue, marker*, s400, zorder10, edgecolorsblack, lw2, labelOptimal Location) ax.scatter(x0_w, y0_w, cgreen, marker^, s150, zorder10, labelWeighted 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, alpha0.3) fig.colorbar(sc, axax, labelPopulation R) ax.axis(equal) fig.tight_layout() fig.savefig(f{OUT}/exp4_location.png, dpi150, bbox_inchestight) 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, mu4*sigma, 300) ax1.plot(xs, norm.pdf(xs, mu, sigma), b-, lw2, labelfSingle part N(20,1.5²)) ax1.fill_between(xs, norm.pdf(xs, mu, sigma), where(xs19)(xs22), alpha0.3, colorblue) ax1.axvline(xx0, colorred, linestyle--, labelfx0{x0:.2f} (10% quantile)) ax1.set_xlabel(Size (cm)); ax1.set_ylabel(Probability Density) ax1.set_title(fSingle Part: P(19≤X≤22){p1*100:.1f}%) ax1.legend() xb np.linspace(mu-4*sigma_xbar, mu4*sigma_xbar, 300) ax2.plot(xb, norm.pdf(xb, mu, sigma_xbar), g-, lw2, labelfSample mean N(20,{sigma_xbar:.2f}²)) ax2.fill_between(xb, norm.pdf(xb, mu, sigma_xbar), where(xb19)(xb22), alpha0.3, colorgreen) ax2.set_xlabel(Sample Mean (cm)); ax2.set_ylabel(Probability Density) ax2.set_title(f25 Samples Mean: P(19≤X̄≤22)≈{p3*100:.1f}%) ax2.legend() fig.tight_layout() fig.savefig(f{OUT}/exp5_normal.png, dpi150, bbox_inchestight) 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, width15, alpha0.7, edgecolorblack, labelEmpirical dist.) ax1.axvline(xQ_opt, colorred, linestyle--, lw2, labelfQ*≈{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-, lw2) ax2.axvline(xQ_best, colorred, linestyle--, labelfQ*{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, alpha0.3) fig.tight_layout() fig.savefig(f{OUT}/exp5_newsboy.png, dpi150, bbox_inchestight) plt.close(fig) print(f 最优订货量: Q≈{Q_opt:.0f}, 期望利润: {max(profits):.2f}) print(\n所有图片生成完毕)