python数学建模第二章

2.1

import pylab as plt
import numpy as np

plt.rc('font',family='SimHei')  #用来正常显示中文标签
plt.rc('axes',unicode_minus=False)  #用来正常显示负号
x=np.linspace(-3,3,50)
y1=np.cosh(x); y2=np.sinh(x); y3=np.exp(x)/2
plt.plot(x,y1,'r-*',label="双曲余弦函数")
plt.plot(x,y2,'--.b',label="双曲正弦函数")
plt.plot(x,y3,'-.dk',label="指数函数")
plt.legend(); plt.show()

2.2

import sympy as sp  
import numpy as np  
import matplotlib.pyplot as plt  
  
# 定义符号  
t, x = sp.symbols('t x')  
  
# 计算不定积分  
integral = sp.integrate(sp.exp(-t) * t**(x-1), t)  
  
# 选择一个x的值进行绘图（例如，x=2）  
x_value = 2  
integral_x_value = integral.subs(x, x_value)  
  
# 使用lambdify将sympy表达式转换为numpy可处理的函数  
integral_func = sp.lambdify(t, integral_x_value, 'numpy')  
  
# 定义t的范围  
t_vals = np.linspace(0, 5, 400)  
  
# 计算对应的y值  
y_vals = integral_func(t_vals)  
  
# 绘图  
plt.figure(figsize=(10, 6))  
plt.plot(t_vals, y_vals, label=f'Integral of $e^{-t}t^{{{x_value-1}}}$')  
plt.title('Integral of $e^{-t}t^{x-1}$ for $x={x_value}$')  
plt.xlabel('t')  
plt.ylabel('Integral value')  
plt.grid(True)  
plt.legend()  
plt.show()

2.3

import numpy as np  
import matplotlib.pyplot as plt  
  
# 定义x的范围  
x = np.linspace(-10, 10, 400)  
  
# 创建一个图形和坐标轴  
plt.figure(figsize=(10, 6))  
ax = plt.gca()  
  
# 循环绘制每条曲线  
colors = ['r', 'g', 'b', 'c', 'm', 'y']  # 定义颜色列表  
for k, color in zip(range(1, 7), colors):  
    y = k * x**2 + 2*k  
    ax.plot(x, y, label=f'$y = {k}x^2 + 2*{k}$', color=color)  
  
# 添加图例  
plt.legend(loc='upper left')  
  
# 添加网格  
plt.grid(True)  
  
# 设置标题和坐标轴标签  
plt.title('Plots of $y = kx^2 + 2k$ for $k = 1, 2, 3, ..., 6$')  
plt.xlabel('x')  
plt.ylabel('y')  
  
# 显示图形  
plt.show()

2.4

import numpy as np  
import matplotlib.pyplot as plt  
  
# 定义x的范围  
x = np.linspace(-10, 10, 400)  
  
# 创建一个2行3列的子图布局  
fig, axs = plt.subplots(2, 3, figsize=(12, 8))  
  
# 遍历每个子图  
for k, ax in enumerate(axs.flat, start=1):  
    y = k * x**2 + 2*k  
    ax.plot(x, y, label=f'$y = {k}x^2 + 2{k}$')  
    ax.set_title(f'k = {k}')  
    ax.legend()  
    ax.grid(True)  
  
# 调整子图之间的间距  
plt.tight_layout()  
  
# 显示图形  
plt.show()

2.5

（1）

import numpy as np  
import matplotlib.pyplot as plt  
from mpl_toolkits.mplot3d import Axes3D  
  
# 定义参数u和v  
u = np.linspace(-2, 2, 400)  
v = np.linspace(0, 2 * np.pi, 400)  
U, V = np.meshgrid(u, v)  
  
# 根据参数方程计算x, y, z  
x = np.sqrt(1 + U**2 + V**2) * np.cos(V)  
y = np.sqrt(1 + U**2 + V**2) * np.sin(V)  
z = U  
  
# 创建图形和坐标轴  
fig = plt.figure()  
ax = fig.add_subplot(111, projection='3d')  
  
# 绘制曲面  
ax.plot_surface(x, y, z, cmap='viridis', edgecolor='none')  
  
# 设置坐标轴标签  
ax.set_xlabel('X')  
ax.set_ylabel('Y')  
ax.set_zlabel('Z')  
  
# 显示图形  
plt.show()

（2）

import numpy as np  
import matplotlib.pyplot as plt  
from mpl_toolkits.mplot3d import Axes3D  
  
# 定义参数u和v  
u = np.linspace(-2, 2, 400)  
v = np.linspace(-2, 2, 400)  
U, V = np.meshgrid(u, v)  
  
# 根据方程计算x, y, z  
X, Y = U, V  
Z = X**2 + Y**2  
  
# 创建图形和坐标轴  
fig = plt.figure()  
ax = fig.add_subplot(111, projection='3d')  
  
# 绘制曲面  
ax.plot_surface(X, Y, Z, cmap='viridis', edgecolor='none')  
  
# 设置坐标轴标签  
ax.set_xlabel('X')  
ax.set_ylabel('Y')  
ax.set_zlabel('Z')  
  
# 显示图形  
plt.show()

2.6

import numpy as np  
import matplotlib.pyplot as plt  
from mpl_toolkits.mplot3d import Axes3D  
  
# 模拟高程数据（假设数据已经过某种方式插值或生成）  
# 这里我们创建一个简单的40x50网格，并填充随机高程值  
x = np.linspace(0, 43.65, 40)  
y = np.linspace(0, 58.2, 50)  
X, Y = np.meshgrid(x, y)  
Z = np.sin(np.sqrt(X**2 + Y**2)) * 100  # 使用一个简单的函数来生成高程数据  
  
# 绘制三维表面图  
fig = plt.figure(figsize=(12, 8))  
ax = fig.add_subplot(121, projection='3d')  
ax.plot_surface(X, Y, Z, cmap='viridis', edgecolor='none')  
ax.set_xlabel('X (km)')  
ax.set_ylabel('Y (km)')  
ax.set_zlabel('Elevation (m)')  
ax.set_title('3D Surface Plot of Elevation Data')  
  
# 绘制等高线图  
plt.subplot(122)  
CS = plt.contour(X, Y, Z, colors='k')  
plt.clabel(CS, inline=1, fontsize=10)  
  
# 标注点 A(30,0) 和 B(43,30)  
# 注意：由于X和Y是网格坐标，我们需要找到最接近这些值的索引  
idx_a_x = np.argmin(np.abs(x - 30))  
idx_a_y = np.argmin(np.abs(y - 0))  
idx_b_x = np.argmin(np.abs(x - 43))  
idx_b_y = np.argmin(np.abs(y - 30))  
  
plt.plot(x[idx_a_x], y[idx_a_y], 'ro', markersize=5, label='A(30,0)')  
plt.plot(x[idx_b_x], y[idx_b_y], 'go', markersize=5, label='B(43,30)')  
plt.xlabel('X (km)')  
plt.ylabel('Y (km)')  
plt.title('Contour Plot of Elevation Data with Points A and B')  
plt.legend()  
  
# 计算地表面积的近似值（忽略地形起伏）  
real_area = 43.65 * 58.2  
print(f"Actual Surface Area (ignoring elevation changes): {real_area} km^2")  
  
# 显示图形  
plt.tight_layout()  
plt.show()

2.7

import numpy as np  
  
# 定义系数矩阵A和常数项向量b  
A = np.array([[4, 2, -1],  
              [3, -1, 2],  
              [11, 3, 0]])  
b = np.array([2, 10, 8])  
  
# 使用numpy的lstsq求解最小二乘解  
x, residuals, rank, s = np.linalg.lstsq(A, b, rcond=None)  
  
print("最小二乘解为：")  
print(x)  
  
# 打印残差和矩阵A的秩  
print("残差为：", residuals)  
print("矩阵A的秩为：", rank)  

import numpy as np  
  
# 定义系数矩阵A和常数项向量b  
A = np.array([[2, 3, 1],  
              [1, -2, 4],  
              [3, 8, -2],  
              [4, -1, 9]])  
b = np.array([4, -5, 13, -6])  
  
# 使用numpy的lstsq函数求解最小二乘解  
# 对于这个特定的问题，由于方程数和未知数数量相同，且没有矛盾，lstsq将给出唯一解  
x, residuals, rank, s = np.linalg.lstsq(A, b, rcond=None)  
  
# 输出解  
print("解 x:", x)  
  
# 验证解是否正确（可选）  
print("验证结果:", np.dot(A, x))  
  
# 检查是否精确等于b（对于精确解，这应该非常接近）  
print("与b的误差:", np.linalg.norm(np.dot(A, x) - b))  
  
# 计算系数矩阵的秩（可选，以确认方程组是否有唯一解）  
print("系数矩阵的秩:", np.linalg.matrix_rank(A))  
  
# 由于秩等于未知数数量，且没有矛盾，我们可以确信有一个唯一解