习题2

学号后四位3026 1班
习题

2.1
import numpy as np
import matplotlib.pyplot as plt

生成x数据点
x = np.linspace(-5, 5, 400)

计算y值
y_cosh = np.cosh(x)
y_sinh = np.sinh(x)
y_exp = 0.5 * np.exp(x)

创建图形和坐标轴
plt.figure(figsize=(10, 6))

绘制cosh(x)
plt.plot(x, y_cosh, label=r'$y = \cosh(x)$', color='blue')

绘制sinh(x)
plt.plot(x, y_sinh, label=r'$y = \sinh(x)$', color='green')

绘制1/2 * e^x
plt.plot(x, y_exp, label=r'$y = \frac{1}{2}e^x$', color='red')

添加图例
plt.legend(loc='upper left')

添加网格
plt.grid(True)

设置x轴和y轴标签
plt.xlabel('x')
plt.ylabel('y')

设置标题
plt.title('Graphs of cosh(x), sinh(x), and 1/2*e^x')

显示图形
plt.show()

2.2
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import gamma

定义x的范围
x = np.linspace(0.01, 10, 400) # Gamma函数在x=0处未定义，所以我们从一个小的正数开始

计算Gamma函数的值
y = gamma(x)

绘制图形
plt.figure(figsize=(10, 6))
plt.plot(x, y, label=r'$\Gamma(x)$')
plt.title('Graph of $\Gamma(x)$ for $x \in [0.01, 10]$')
plt.xlabel('x')
plt.ylabel(r'$\Gamma(x)$')
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) # 从-10到10，共400个点

创建一个图形和坐标轴
plt.figure(figsize=(10, 6))

循环绘制6条曲线
for k in range(1, 7):
y = k * (x ** 2) + 2 * k # 计算y值
plt.plot(x, y, label=f'k={k}') # 绘制曲线并添加标签

添加图例
plt.legend()

设置坐标轴标签
plt.xlabel('x')
plt.ylabel('y')

设置标题
plt.title('Plot of y_k(x^2) + 2k for k = 1 to 6')

显示图形
plt.grid(True) # 显示网格
plt.show()

2.4
import matplotlib.pyplot as plt
import numpy as np

创建一个2行3列的子图
fig, axs = plt.subplots(2, 3, figsize=(12, 8))

定义x的范围
x = np.linspace(-10, 10, 400)

遍历k的值，绘制曲线
for k, ax in enumerate(axs.flat, start=1):
y = k * (x ** 2) + 2 * k
ax.plot(x, y, label=f'k={k}')
ax.set_title(f'y = {k}x^2 + 2{k}')
ax.legend()
ax.grid(True)

调整子图之间的间距
plt.tight_layout()

显示图形
plt.show()

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

定义x和y的范围
x = np.linspace(-5, 5, 400)
y = np.linspace(-5, 5, 400)
X, Y = np.meshgrid(x, y)

根据双曲面方程计算z的值
注意：这里我们需要计算z的两个解（因为双曲面有两个分支），但这里我们只取其中一个分支
Z = np.sqrt((8 / 4) * X ** 2 + (8 / 10) * Y ** 2 - 8)

创建一个图形和一个3D坐标轴
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

绘制双曲面的线框图
ax.plot_wireframe(X, Y, Z, color='b')

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

显示图形
plt.show()

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

定义x和y的范围
x = np.linspace(-5, 5, 400)
y = np.linspace(-5, 5, 400)
X, Y = np.meshgrid(x, y)

根据椭圆抛物面方程计算z的值
Z = (X ** 2) / 4 + (Y ** 2) / 6

创建一个图形和一个3D坐标轴
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

绘制椭圆抛物面
surf = ax.plot_surface(X, Y, Z, cmap='viridis', edgecolor='none')

添加颜色条
fig.colorbar(surf, ax=ax, shrink=0.5, aspect=5)

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

设置坐标轴的显示范围（可选，根据需要调整）
ax.set_xlim([-5, 5])
ax.set_ylim([-5, 5])
ax.set_zlim([0, 15]) # 根据需要调整z轴的范围

显示图形
plt.show()

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

假设elevations是已经加载的高程数据，其形状为(582, 4356)
这里我们需要将x和y坐标转换为km单位，假设数据已经是km单位或者进行了相应的缩放
x = np.linspace(0, 58.2, 4356)
y = np.linspace(0, 43.56, 582)
X, Y = np.meshgrid(x, y)

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
surf = ax.plot_surface(X, Y, elevations, cmap='viridis', edgecolor='none')

ax.set_xlabel('X (km)')
ax.set_ylabel('Y (km)')
ax.set_zlabel('Elevation (km)')

plt.show()

fig, ax = plt.subplots()
CS = ax.contour(X, Y, elevations, colors='k')
ax.clabel(CS, inline=1, fontsize=10)

标记A和B点
ax.plot(30, 0, 'ro', markersize=5, label='A(30,0)')
ax.plot(43, 30, 'ro', markersize=5, label='B(43,30)')
ax.legend()

ax.set_xlabel('X (km)')
ax.set_ylabel('Y (km)')

plt.show()

假设dx和dy是x和y方向上的网格间距
dx = x[1] - x[0]
dy = y[1] - y[0]

计算每个网格单元的平均高程
avg_elevations = 0.5 * (elevations[:-1, :-1] + elevations[1:, 1:])

近似计算表面积（这里只计算了顶面，忽略了底面）
注意：这只是一个非常粗略的近似
approx_surface_area = np.sum(avg_elevations * dx * dy)

print(f"Approximate surface area: {approx_surface_area} km^2")

2.7
import numpy as np

定义系数矩阵A和常数项向量b
A = np.array([[4, 2, -1],
[3, -1, 2],
[0, 0, 0]]) # 第三行不影响前两个未知数的解
b = np.array([2, 10, 0]) # 第三个方程的常数项设为0

仅考虑前两个方程来求解x1和x2
A_reduced = A[:2, :2] # 取前两行和前两列
b_reduced = b[:2]

使用numpy.linalg.solve求解前两个方程
x_reduced = np.linalg.solve(A_reduced, b_reduced)

x3的值是任意的，因为它没有被方程约束
x3 = np.nan # 表示x3是未知的

打印解
print("解为：")
print("x1 =", x_reduced[0])
print("x2 =", x_reduced[1])
print("x3 是任意的（在这个方程组中没有足够的方程来确定它）")

2.7-2
import numpy as np

定义系数矩阵A和常数项向量b
A = np.array([[4, 2, -1],
[3, -1, 2],
[0, 0, 0]]) # 第三行不影响前两个未知数的解
b = np.array([2, 10, 0]) # 第三个方程的常数项设为0

仅考虑前两个方程来求解x1和x2
A_reduced = A[:2, :2] # 取前两行和前两列
b_reduced = b[:2]

使用numpy.linalg.solve求解前两个方程
x_reduced = np.linalg.solve(A_reduced, b_reduced)

x3的值是任意的，因为它没有被方程约束
x3 = np.nan # 表示x3是未知的

打印解
print("解为：")
print("x1 =", x_reduced[0])
print("x2 =", x_reduced[1])
print("x3 是任意的（在这个方程组中没有足够的方程来确定它）")

2.8
import numpy as np

构建系数矩阵 A
n = 1000 # 方程组的变量数
A = np.zeros((n, n))
for i in range(n):
A[i, i] = 4 # 主对角线上的元素
if i < n - 1:
A[i, i + 1] = 1 # 次对角线上的元素
if i > 0:
A[i, i - 1] = 1 # 上次对角线上的元素

构建常数项向量 b
b = np.arange(1, n + 1)

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

打印前几个解以验证
print(x[:10]) # 打印前10个解

2.9
import sympy as sp

定义符号变量
x, y = sp.symbols('x y')

定义方程组
equation1 = sp.Eq(x ** 2 - y - x, 3)
equation2 = sp.Eq(x + 3 * y, 2)

解方程组
solutions = sp.solve((equation1, equation2), (x, y), dict=True)

print(solutions)

2.10
import math

圆锥体底面半径和高度
R_cone = 2 # 圆锥体底面半径
h_cone = 2 # 圆锥体高度（从y=1到y=3，但考虑到整个旋转体，我们取2）

水的密度和重力加速度
rho = 1000 # kg/m^3
g = 9.81 # m/s^2

近似计算将容器内水移出所需做的功
我们将容器视为一个高度为2、底面半径为2的圆锥体
并使用简单的近似方法
简化计算，将圆锥体分为多个薄层
num_layers = 10
layer_height = h_cone / num_layers # 注意这里使用h_cone作为总高度
total_work = 0

for i in range(num_layers):
y_mid = 1 + (i + 0.5) * layer_height # 每一层的中间高度
radius = R_cone * (y_mid - 1) / h_cone # 每一层的半径（基于圆锥体的线性变化）
area = math.pi * radius ** 2 # 每一层的面积
volume = area * layer_height # 每一层的体积
mass = rho * volume # 每一层的水的质量
work_for_layer = mass * g * (3 - y_mid) # 注意这里的高度差可能需要根据实际情况调整
total_work += work_for_layer

注意：上面的计算中，3 - y_mid 可能是不准确的，因为容器的顶部不是在y=3处
如果容器的顶部确实在y=3（且容器是封闭的），则应该使用其他方法来计算重力势能的变化
但由于这是一个近似的圆锥体，我们可以暂时接受这个计算
print(f"将容器内水移出所需做的功（近似值）: {total_work} 焦耳")

2.11
import numpy as np

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)

假设的初始值
x1, x2, y1, y2 = 0.5, 1.5, 0.2, 0.8

计算函数值
fy1 = f(y1)
gy2 = g(y2)
fx1 = f(x1)
gx2 = g(x2)

构建方程组
A = np.array([[2, 0, -3 * fy1, -4 * gy2],
[0, 3, -2 * fy1, -6 * gy2],
[fx1 + 3 * gx2 - 3, 0, 1, 0],
[0, 4 * fx1 + 6 * gx2 - 1, 0, 5]])
b = np.array([1, 2, y1, 5 * y2])

求解方程组
try:
solution = np.linalg.solve(A, b)
print("解为:", solution)
except np.linalg.LinAlgError:
print("方程组无解或系数矩阵奇异")

2.12
import sympy as sp

定义符号变量
x, y, z = sp.symbols('x y z')

定义矩阵
A_sym = sp.Matrix([[-1, 1, 0],
[-4, 3, 0],
[1, 0, 2]])

计算特征多项式
char_poly = A_sym.charpoly(x)

打印特征多项式
print("特征多项式:")
print(char_poly)

求解特征多项式
eigenvalues_sym = sp.solve(char_poly, x)

打印符号特征值
print("符号特征值:")
print(eigenvalues_sym)

对于每个符号特征值，计算对应的特征向量
eigenvectors_sym = []
for eigenval in eigenvalues_sym:
I = sp.eye(3) # 3x3单位矩阵
B = A_sym - eigenval * I
null_space = B.nullspace() # 计算B的零空间（即特征向量）
eigenvectors_sym.append(null_space[0][:]) # 提取第一个（或唯一的）零空间向量

打印符号特征向量
print("符号特征向量:")
for vec in eigenvectors_sym:
print(vec)

2.13
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 residuals(vars):
x1, x2, y1, y2 = vars
r1 = 2 * x1 - (3 * f(y1) + 4 * g(y2) - 1)
r2 = 3 * x2 - (2 * f(y1) + 6 * g(y2) - 2)
r3 = y1 - (f(x1) + 3 * g(x2) - 3)
r4 = 5 * y2 - (4 * f(x1) + 6 * g(x2) - 1)
r5 = x1 + y1 - (f(y2) + g(x2) - 2)
r6 = x2 - 3 * y2 - (2 * f(x1) - 10 * g(y1) - 5)
return np.array([r1, r2, r3, r4, r5, r6])

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

使用least_squares求解最小二乘解
result = least_squares(residuals, initial_guess)

输出结果
print("最小二乘解:", result.x)