《Python数学建模算法与应用》第二章习题解答

代码仅供参考,如有谬误或不足之处请留言。

1.同一图形界面画多条曲线并标注

例1:

示例代码:

import math
import pylab as plt
import numpy as np

plt.rc('text', usetex=True)  # 调用字库
x = np.linspace(-10, 10, 100)
y1 = np.cosh(x)
y2 = np.sinh(x)
y3 = math.e**x/2
plt.plot(x, y1, label='$\\mathrm{cosh}(x)$')
plt.plot(x, y2, label='$\\mathrm{sinh}(x)$')
plt.plot(x, y3, label='$\\frac{1}{2} \cdot e^x$')
plt.legend()
plt.xlabel('$x$')
plt.ylabel('$y$')
plt.show()

运行结果: 

例2: 

示例代码:

import numpy as np
import matplotlib.pyplot as plt

k_values = [1, 2, 3, 4, 5, 6]  # k的取值
x = np.linspace(-10, 10, 100)  # x的范围

for k in k_values:
    y = k * x ** 2 + 2 * k  # 计算y值
    label = f'k={k}'  # 设置标注
    plt.plot(x, y, label=label)  # 绘制曲线

plt.xlabel('x')  # 添加x轴标签
plt.ylabel('y')  # 添加y轴标签
plt.legend()  # 添加图例
plt.grid(True)  # 添加网格线
plt.show()  # 显示图形

 运行结果:

2. 积分函数图像

示例代码:

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import quad

def fun(t, x):
    return np.exp(-t) * (t ** (x - 1))


x = np.linspace(0, 10, 100)  # x 的范围
y = [quad(fun, 0, np.inf, args=i)[0] for i in x]  # 计算积分

plt.plot(x, y)
plt.xlabel('x')
plt.ylabel('$ y = \int_0^{\infty} e^{-t} \cdot t^{x-1} dt $')
plt.grid(True)
plt.show()

 运行结果:

3.分别在子图中绘制曲线 

 

示例代码:

import numpy as np
import matplotlib.pyplot as plt

plt.rc('font', family='SimHei')  # 用来正常显示中文标签
plt.rc('axes', unicode_minus=False) # 用来正常显示负号
k_values = [1, 2, 3, 4, 5, 6]  # k的取值
x = np.linspace(-10, 10, 100)  # x的范围

fig, axs = plt.subplots(2, 3, figsize=(10, 6))  # 创建2行3列的子窗口

for i, k in enumerate(k_values):
    y = k * x ** 2 + 2 * k  # 计算y值
    row = i // 3  # 计算子窗口所在的行数
    col = i % 3  # 计算子窗口所在的列数
    ax = axs[row, col]  # 获取当前子窗口
    label = f'k={k}'  # 设置标注
    ax.plot(x, y, label=label)  # 绘制曲线
    ax.set_xlabel('x')  # 添加x轴标签
    ax.set_ylabel('y')  # 添加y轴标签
    ax.set_title(f'图 {i+1}')  # 添加子窗口标题
    ax.legend()  # 添加图例
    ax.grid(True)  # 添加网格线

plt.tight_layout()  # 自动调整子窗口布局
plt.show()  # 显示图形

运行结果:

4.绘制二次曲面 

(1)

示例代码:

import numpy as np
import matplotlib.pyplot as plt

a = 2
b = np.sqrt(10)
c = np.sqrt(8)

phi = np.arange(0, 2*np.pi+0.1, 0.1)
theta = np.arange(-1, 1.1, 0.1)[:, np.newaxis]

x = a * np.cosh(theta) * np.cos(phi)
y = b * np.cosh(theta) * np.sin(phi)
z = c * np.sinh(theta)

fig = plt.figure()
ax = fig.add_subplot(projection='3d')
ax.plot_surface(x, y, z,  cmap='viridis')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
ax.set_box_aspect([1, 1, 1])
ax.set_title('$\\frac{x^2}{4}+\\frac{y^2}{10}-\\frac{z^2}{8}=1$')
plt.show()

运行结果:

(2)

示例代码:

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# 创建一个三维坐标系
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

# 生成x,y的网格
x = np.linspace(-50, 50, 1000)
y = np.linspace(-50, 50, 1000)
X, Y = np.meshgrid(x, y)

# 计算z的值
Z = (X**2)/4+(Y**2)/6

# 绘制二次曲面
ax.plot_surface(X, Y, Z, cmap='viridis')

# 添加坐标轴标签
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')

# 添加标题
ax.set_title('$(\\frac{x^2}{4}+\\frac{y^2}{6})=z$')

# 设置z轴的范围
ax.set_zlim(0, 1000)

# 显示图形
plt.show()

运行结果:

 5.读取Excel中的数据并绘制网格图和等高线图

        由于笔者没有附件资源,无法获悉Excel表格中的内容,该题目以自制的Excel表格为例。

        [1:,0]为x坐标0到100,间隔为1;[0,1:]为y坐标0到100,间隔为1;[1:,1:]为z坐标,是0到1000的随机数。由于高度数据是随机生成的,所以绘制的结果不符合常理,只需要掌握绘制方法,在实际应用过程中替换数据即可。

生成Excel文件的代码:

import pandas as pd
import numpy as np

# 生成x、y和z数据
x = np.arange(0, 101, 1)
y = np.arange(0, 101, 1)
z = np.random.randint(0, 1001, size=(101, 101))

# 创建DataFrame对象
df = pd.DataFrame(data=z, index=x, columns=y)

# 写入Excel文件
df.to_excel('data.xlsx')

绘制等高线: 

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import RectBivariateSpline


df = pd.read_excel('data.xlsx', header=None)

x = df.iloc[1:, 0].values
y = df.iloc[0, 1:].values
z = df.iloc[1:, 1:].values

# 绘制等高线图
plt.contour(x, y, z)
plt.xlabel('x')
plt.ylabel('y')
plt.title('Contour Plot')
plt.colorbar()  # 添加颜色条
# 在已知点处添加标注
point1 = (30, 0)
point2 = (43, 30)
plt.annotate('(30,0)', point1, textcoords="offset points", xytext=(0,10), ha='center')
plt.annotate('(43,30)', point2, textcoords="offset points", xytext=(0,10), ha='center')
# 创建插值函数
interp_func = RectBivariateSpline(x, y, z)

# 定义积分区间
x_min, x_max = min(x), max(x)
y_min, y_max = min(y), max(y)

# 计算积分区间的网格点数
grid_size = 100  # 可根据需要调整

# 生成网格点
x_grid = np.linspace(x_min, x_max, grid_size)
y_grid = np.linspace(y_min, y_max, grid_size)

# 计算网格点上的插值值
z_grid = interp_func(x_grid, y_grid)

# 计算曲面面积
area = np.trapz(np.trapz(z_grid, x_grid), y_grid)

print('区域面积:', area)

plt.show()

运行结果:

区域面积: 2491714.2083680276

绘制网格图:

import pandas as pd
import matplotlib.pyplot as plt
import numpy as np

df = pd.read_excel('data2.xlsx', header=None)

x = df.iloc[1:, 0].values
y = df.iloc[0, 1:].values
z = df.iloc[1:, 1:].values

# 创建网格
X, Y = np.meshgrid(y, x)

# 创建画布和3D坐标轴
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

# 绘制三维网格图
ax.plot_surface(X, Y, z)

# 设置坐标轴标签
ax.set_xlabel('Y')
ax.set_ylabel('X')
ax.set_zlabel('Z')

# 显示三维网格图
plt.show()

运行结果: 

 6、求线性方程组的唯一解、最小二乘解、最小范数解

例1:

示例代码:

仅给出第一问的求解过程,第二问将代码中的矩阵替换即可。

import numpy as np

# 定义系数矩阵 A 和常数向量 b
A = np.array([[4, 2, -1], [3, -1, 2], [11, 3, 0]])
b = np.array([2, 10, 8])

# 解线性方程组 Ax = b
x = np.linalg.solve(A, b)

# 判断解的情况
if np.linalg.matrix_rank(A) == np.linalg.matrix_rank(np.column_stack((A, b))):
    if np.linalg.matrix_rank(A) == A.shape[1]:
        print("线性方程组有唯一解")
        x_unique = np.linalg.solve(A, b)
        print("唯一解 x =", x_unique)
    else:
        print("线性方程组有无穷多解")
        x_least_squares = np.linalg.lstsq(A, b, rcond=None)[0]
        print("最小二乘解 x =", x_least_squares)
else:
    print("线性方程组无解")
    x_least_squares = np.linalg.lstsq(A, b, rcond=None)[0]
    print("最小二乘解 x =", x_least_squares)
    x_min_norm = np.linalg.pinv(A).dot(b)
    print("最小范数解 x =", x_min_norm)

 运行结果:

线性方程组无解
最小二乘解 x = [ 1.21304348 -1.44782609  1.95652174]
最小范数解 x = [ 1.21304348 -1.44782609  1.95652174]

例2:

import numpy as np

# 生成系数矩阵A
A = np.zeros((1000, 1000))
np.fill_diagonal(A, 4)
np.fill_diagonal(A[:, 1:], 1)
np.fill_diagonal(A[1:, :], 1)
# 生成常数向量b
b = np.arange(1, 1001)
# 判断解的情况
if np.linalg.matrix_rank(A) == np.linalg.matrix_rank(np.column_stack((A, b))):
    if np.linalg.matrix_rank(A) == A.shape[1]:
        print("线性方程组有唯一解")
        x_unique = np.linalg.solve(A, b)
        print("唯一解 x =", x_unique)
    else:
        print("线性方程组有无穷多解")
        x_least_squares = np.linalg.lstsq(A, b, rcond=None)[0]
        print("最小二乘解 x =", x_least_squares)
else:
    print("线性方程组无解")
    x_least_squares = np.linalg.lstsq(A, b, rcond=None)[0]
    print("最小二乘解 x =", x_least_squares)
    x_min_norm = np.linalg.pinv(A).dot(b)
    print("最小范数解 x =", x_min_norm)


有唯一解,运行结果太长,这里不做展示。

7.求非线性方程组的符号解和数值解

例1:

 示例代码:

from sympy import symbols, Eq, solve

# 定义未知数
x, y = symbols('x y')

# 定义非线性方程组
equations = (Eq(x**2 - y - x - 3, 0), Eq(x + 3*y - 2, 0))

# 求解符号解
symbolic_solution = solve(equations, (x, y))
print("符号解:", symbolic_solution)

# 数值解需要使用数值方法求解
from scipy.optimize import fsolve

# 定义非线性方程组的函数
def equations(variables):
    x, y = variables
    return [x**2 - y - x - 3, x + 3*y - 2]

# 初始猜测值
initial_guess = [1, 1]

# 求解数值解
numeric_solution = fsolve(equations, initial_guess)
print("数值解:", numeric_solution)

运行结果:

符号解: [(1/3 - sqrt(34)/3, 5/9 + sqrt(34)/9), (1/3 + sqrt(34)/3, 5/9 - sqrt(34)/9)]
数值解: [ 2.27698396 -0.09232799]

例2:

 示例代码:

from sympy import symbols, pi, integrate, sqrt

# 定义符号变量
y = symbols('y')

# 定义曲线方程
curve1 = y - sqrt(4*y - y**2)
curve2 = sqrt(4 - y)

# 计算容器的体积
volume = pi * integrate((curve1**2 - curve2**2), (y, 1, 3))
print("容器的体积:", volume)

# 定义重力加速度和水的密度
g = 9.8
p = 1000

# 计算至少需要做多少功
work = p * g * volume
print("至少需要做多少功:", work)

运行结果:

容器的体积: pi*(-8*pi/3 - 4*sqrt(3) + 12)
至少需要做多少功: 9800.0*pi*(-8*pi/3 - 4*sqrt(3) + 12)

例3: 

 

示例代码:

import numpy as np
from scipy.optimize import fsolve

def f(x):
    return (abs(x + 1) - abs(x - 1)) / 2 + np.sin(x)

def g(x):
    return (abs(x + 3) - abs(x - 3)) / 2 + np.cos(x)

def equations(variables):
    x1, x2, y1, y2 = variables[:4]
    eq1 = 2 * x1 - (3 * f(y1) + 4 * g(y2) - 1)
    eq2 = 3 * x2 - (2 * f(y1) + 6 * g(y2) - 2)
    eq3 = y1 - (f(x1) + 3 * g(x2) - 3)
    eq4 = 5 * y2 - (4 * f(x1) + 6 * g(x2) - 1)

    return [eq1, eq2, eq3, eq4]

initial_guess = [0, 0, 0, 0]  # 初始猜测值

numeric_solution = fsolve(equations, initial_guess)
print("数值解:", numeric_solution)

 运行结果:

数值解: [-0.85877135  0.3022151  -0.84512781  0.01562419]

例4:

示例代码: 

import numpy as np
from scipy.optimize import least_squares

def f(x):
    return (np.abs(x + 1) - np.abs(x - 1)) / 2 + np.sin(x)

def g(x):
    return (np.abs(x + 3) - np.abs(x - 3)) / 2 + np.cos(x)

def equations(variables):
    x1, x2, y1, y2 = variables[:4]
    eq1 = 2 * x1 - (3 * f(y1) + 4 * g(y2) - 1)
    eq2 = 3 * x2 - (2 * f(y1) + 6 * g(y2) - 2)
    eq3 = y1 - (f(x1) + 3 * g(x2) - 3)
    eq4 = 5 * y2 - (4 * f(x1) + 6 * g(x2) - 1)
    eq5 = x1 + y1 - (f(y2) + g(x2) - 2)
    eq6 = x2 - 3 * y2 - (2 * f(x1) - 10 * g(y1) - 5)

    return [eq1, eq2, eq3, eq4, eq5, eq6]

initial_guess = [0, 0, 0, 0]  # 初始猜测值

result = least_squares(equations, initial_guess)
numeric_solution = result.x
print("最小二乘解:", numeric_solution)

运行结果:

最小二乘解: [-0.69842439  0.10960544 -1.15185365  0.07400321]

8.求矩阵特征值和特征向量的数值解和符号解

 示例代码:

import numpy as np

# 定义矩阵
matrix = np.array([[-1, 1, 0], [-4, 3, 0], [1, 0, 2]])

# 求解特征值和特征向量的数值解
eigenvalues, eigenvectors = np.linalg.eig(matrix)

print("特征值:", eigenvalues)
print("特征向量:\n", eigenvectors)

# 求解特征值和特征向量的符号解
symbolic_eigenvalues = np.linalg.eigvals(matrix)
symbolic_eigenvectors = np.linalg.eig(matrix)[1]

print("特征值的符号解:", symbolic_eigenvalues)
print("特征向量的符号解:\n", symbolic_eigenvectors)

运行结果:

特征值: [2. 1. 1.]
特征向量:
 [[ 0.          0.40824829  0.40824829]
 [ 0.          0.81649658  0.81649658]
 [ 1.         -0.40824829 -0.40824829]]
特征值的符号解: [2. 1. 1.]
特征向量的符号解:
 [[ 0.          0.40824829  0.40824829]
 [ 0.          0.81649658  0.81649658]
 [ 1.         -0.40824829 -0.40824829]]

  • 17
    点赞
  • 17
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值