海洋工程建模基础

本文最后更新于 2025年1月11日 下午

作业一

作业一
作业一

我写的作业，第一题答案算错了，没改回来：作业一

作业二

1. 假设6个矩阵A1,A2,A3,A4,A5,A6,的阶数分别为6x7,7x3,3x1,1x2,2x4,4x5,求它们相乘的所需元素乘法的最小个数以及对应的最优乘法次序。

1
2
3

2.

\[ m,n为非负整数，p(m,n)=\begin{cases} 0 & \text{若 } m=0, n \neq 0 \\ 1 & \text{若 } m \neq 0, n=0 \\ \frac{p(m,n-1)}{2}+\frac{p(m-1,n+1)}{2} & \text{若 } m \geq 1,n \geq 1 \end{cases} \]

  1. 分别在 \[ n=1,2....,50 \] 时,(用动态规划方法)求 \(m=f(n)\) 满足 \(|p(m,n)-0.5|\)最小;
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
import numpy as np

input_n=50

def calculate_optimal_mAndp(n_array):
    max_n = max(n_array)
    m = 1
    # 初始化 p 矩阵，大小为 (m+1) x (max_n + m + 1)
    p = np.zeros((m + 1, max_n + m + 1))

    # 设置初始条件
    p[0, 1:max_n + 1] = 0  # 当 m=0, n != 0
    p[1, 0] = 1            # 当 m != 0, n=0
    for d in range(1, max_n + 1):
        p[1, d] = (p[1, d - 1] + p[0, d + 1]) / 2  # 填充 p(1, n)

    # 动态规划填充 p 矩阵直到满足条件
    while np.max(p[:, max_n]) < 0.5:
        m += 1
        # 补全矩阵 p 的第1到(m-1)+1行
        p = np.pad(p, ((0, 1), (0, 1)), 'constant')  # 扩展 p 矩阵
        p[0, m + max_n] = 0
        for l in range(1, m):
            p[l, m + max_n - l] = (p[l, m + max_n - l - 1] + 
                                   p[l - 1, m + max_n - l + 1]) / 2

        # 补全矩阵 p 的第 m + 1 行
        p[m, 0] = 1
        for d in range(1, max_n + 1):
            p[m, d] = (p[m, d - 1] + p[m - 1, d + 1]) / 2

    # 计算答案并填充结果
    results = np.zeros((3, len(n_array)))
    for i, n in enumerate(n_array):
        answer = 1
        min_diff = float('inf')  # 初始化最小差值为无穷大
        best_m = 1  # 跟踪最优的 m 值

        # 查找满足 p(m,n) >= 0.5 的 m
        while p[answer, n] < 0.5:
            answer += 1

        # 遍历当前 m 值及其前一个 m 值，寻找最优 m
        for current_m in (answer - 1, answer):
            if current_m >= 1:
                diff = abs(p[current_m, n] - 0.5)
                if diff < min_diff:
                    min_diff = diff
                    best_m = current_m

        # 更新结果
        results[0, i] = n
        results[1, i] = best_m
        results[2, i] = p[best_m, n]

    return p, results

# 调用
n_array = np.array(range(1, input_n+1))  # 1到50
p, results = calculate_optimal_mAndp(n_array)

# 输出结果
for n, m, p in zip(results[0], results[1], results[2]):
    print(f"n={int(n)}, Optimal m={int(m)}, p(m,n)={p:.8f}")

n=1, Optimal m=1, p(m,n)=0.50000000
n=2, Optimal m=4, p(m,n)=0.50781250
n=3, Optimal m=9, p(m,n)=0.50344467
n=4, Optimal m=16, p(m,n)=0.49955983
n=5, Optimal m=25, p(m,n)=0.49661744
n=6, Optimal m=37, p(m,n)=0.49989688
n=7, Optimal m=51, p(m,n)=0.50076639
n=8, Optimal m=67, p(m,n)=0.50062926
n=9, Optimal m=85, p(m,n)=0.50005461
n=10, Optimal m=105, p(m,n)=0.49929721
n=11, Optimal m=128, p(m,n)=0.50009713
n=12, Optimal m=153, p(m,n)=0.50038148
n=13, Optimal m=180, p(m,n)=0.50035225
n=14, Optimal m=209, p(m,n)=0.50013030
n=15, Optimal m=240, p(m,n)=0.49979040
n=16, Optimal m=274, p(m,n)=0.50014290
n=17, Optimal m=310, p(m,n)=0.50028242
n=18, Optimal m=348, p(m,n)=0.50027180
n=19, Optimal m=388, p(m,n)=0.50015490
n=20, Optimal m=430, p(m,n)=0.49996279
n=21, Optimal m=475, p(m,n)=0.50016021
n=22, Optimal m=521, p(m,n)=0.49983998
n=23, Optimal m=570, p(m,n)=0.49986929
n=24, Optimal m=622, p(m,n)=0.50016584
n=25, Optimal m=675, p(m,n)=0.50004254
n=26, Optimal m=730, p(m,n)=0.49988010
n=27, Optimal m=788, p(m,n)=0.49995576
n=28, Optimal m=848, p(m,n)=0.49997180
n=29, Optimal m=910, p(m,n)=0.49993981
n=30, Optimal m=975, p(m,n)=0.50008585
n=31, Optimal m=1041, p(m,n)=0.49997036
n=32, Optimal m=1110, p(m,n)=0.50002161
n=33, Optimal m=1181, p(m,n)=0.50003133
n=34, Optimal m=1254, p(m,n)=0.50000646
n=35, Optimal m=1329, p(m,n)=0.49995279
n=36, Optimal m=1407, p(m,n)=0.50002567
n=37, Optimal m=1487, p(m,n)=0.50006262
n=38, Optimal m=1568, p(m,n)=0.49993388
n=39, Optimal m=1653, p(m,n)=0.50004914
n=40, Optimal m=1739, p(m,n)=0.50000707
n=41, Optimal m=1827, p(m,n)=0.49994599
n=42, Optimal m=1918, p(m,n)=0.49997933
n=43, Optimal m=2011, p(m,n)=0.49998875
n=44, Optimal m=2106, p(m,n)=0.49997739
n=45, Optimal m=2204, p(m,n)=0.50004427
n=46, Optimal m=2303, p(m,n)=0.49999500
n=47, Optimal m=2405, p(m,n)=0.50002063
n=48, Optimal m=2509, p(m,n)=0.50002738
n=49, Optimal m=2615, p(m,n)=0.50001748
n=50, Optimal m=2723, p(m,n)=0.49999288
  1. 用蒙特卡罗方法计算\(p(f(50),50)\)，将结果与与(1)中的对应结果比较,若差别明显，则分析原因，采用适当方法改进(1)中的算法，并重新计算.

蒙特卡罗方法计算\(p(m,n)\)的流程:

(1).重复以下过程直到m=0或者n=0

对0-1随机取值，若0则m减1;若1则m加1,n减1 或者取(0,1)间随机数，分为>0.5和<0.5两种情况

(2).重复过程(1)足够多次(比如10000),其中最后n=0的次数比例即为\(p(m,n)\)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
import random
# 蒙特卡罗模拟方法来估计 p(m, n)
def monte_carlo_simulation(m, n, trials):
    success_count = 0
    
    for _ in range(trials):
        current_m, current_n = m, n
        
        while current_n > 0:
            # 模拟随机过程，选择下一个状态
            if random.random() < 0.5:
                current_n -= 1  # 相当于选择了 p(m, n-1)
            else:
                current_m -= 1  # 相当于选择了 p(m-1, n+1)
                current_n += 1
            
            # 如果 m 变成 0，则概率是 0
            if current_m == 0:
                break
        
        # 如果成功达到 n=0，说明当前模拟成功
        if current_n == 0 and current_m > 0:
            success_count += 1
    
    # 返回成功次数占总试验次数的比例作为概率估计
    return success_count / trials


# 根据(1)的结果，假设 f(50) = m_optimal
m_optimal = results[1][-1]         # (1)中 f(50) 对应的 m 值
trials = 10000  # 设置模拟的次数

# 进行蒙特卡罗模拟
monte_carlo_result = monte_carlo_simulation(m_optimal, input_n, trials)

# 打印蒙特卡罗结果
print(f"蒙特卡罗方法估计的 p({m_optimal}, {n}) = {monte_carlo_result}")
蒙特卡罗方法估计的 p(2723.0, 50.0) = 0.5009
  1. 画出\((n,f(n))\)的散列图,分别用2次和3次多项式\(g(x)\)和\(h(x)\)去拟合
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
import numpy as np
import matplotlib.pyplot as plt

# 提取 n 和 f(n)
n_values = results[0]
f_values = results[1]

# 绘制散点图
plt.scatter(n_values, f_values, label='数据点', color='blue')

# 多项式拟合
coeffs_2 = np.polyfit(n_values, f_values, 2)  # 2次多项式
coeffs_3 = np.polyfit(n_values, f_values, 3)  # 3次多项式

# 打印拟合的多项式系数
print(f"2次多项式拟合系数: {coeffs_2}")
print(f"3次多项式拟合系数: {coeffs_3}")

# 创建多项式函数
poly_2 = np.poly1d(coeffs_2)
poly_3 = np.poly1d(coeffs_3)

# 生成用于绘图的 n 值
n_fit = np.linspace(1, 50, 200)
f_fit_2 = poly_2(n_fit)
f_fit_3 = poly_3(n_fit)

# 绘制拟合曲线
plt.plot(n_fit, f_fit_2, label='2次多项式拟合曲线', color='orange',linewidth=2)
plt.plot(n_fit, f_fit_3, label='3次多项式拟合曲线', color='green',linewidth=2, 
                                                    linestyle='--')

# 图形美化
plt.title('散点图和拟合曲线')
plt.xlabel('n')
plt.ylabel('f(n)')
plt.legend()
plt.grid()
plt.show()

2次多项式拟合系数: [ 1.0990881  -0.50199234  0.51367347]
3次多项式拟合系数: [ 2.18411794e-05  1.09741725e+00 -4.67568453e-01  3.60178029e-01]
4

(4)用蒙特卡罗方法计算\(p(g(500),500)\)和\(p(h(500),500)\),通过比较以确定\(g\)和\(h\)哪个较优(比较是否过拟合)

1
2
3
4
5
6
7
8
9
10
11
12
13
#计算g(500)和h(500)
g_500=coeffs_2[0]*500**2+coeffs_2[1]*500+coeffs_2[2]
h_500=coeffs_3[0]*500**3+coeffs_3[1]*500**2+coeffs_3[2]*500+coeffs_3[3]
print("g(500)=",int(g_500))
print("h(500)=",int(h_500))
trials=10000
#用蒙特卡洛方法计算p(g(500),500)和pp(g(500),500)
p_g_500=monte_carlo_simulation(int(g_500),500,trials)
p_h_500=monte_carlo_simulation(int(h_500),500,trials)

print("P(g(500),500)=",p_g_500)
print("P(h(500),500)=",p_h_500)

g(500)= 274521
h(500)= 276851
P(g(500),500)= 0.5019
P(h(500),500)= 0.502

通过蒙特卡罗方法计算\(p(g(500),500)\)和\(p(h(500),500)\)，发现\(|p(h(500),500)-0.5| > |p(g(500),500)-0.5|\), 所以可以认为3次拟合曲线\(h\)过拟合

3.

\[ E_1 = \begin{pmatrix} 1 & -2 \\ 0 & 2 \end{pmatrix}, E_2 = \begin{pmatrix} 1 & 0 \\ 0 & 2 \end{pmatrix}, R= \left\{ \begin{pmatrix} x \\ y \end{pmatrix} \mid max\{|x|,|y|\} = 1 \ \right\} , S= \left\{ \begin{pmatrix} x \\ y \end{pmatrix} \mid x^2 + y^2 = 1 \ \right\} \]

  1. 在同一个平面坐标系中分别画出R, E1R, E2R, E1E2R;
5
6 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
import numpy as np
import matplotlib.pyplot as plt

# 设置matplotlib的字体
plt.rcParams['font.sans-serif'] = ['SimHei']  # 指定默认字体为黑体
plt.rcParams['axes.unicode_minus'] = False  

# 定义矩阵
E1 = np.array([[1, -2], [0, 2]])
E2 = np.array([[1, 0], [0, 2]])

# 定义集合 R 的边界
# R 是一个正方形
R_x = np.array([-1, 1, 1, -1, -1])
R_y = np.array([-1, -1, 1, 1, -1])
R = np.array([R_x, R_y]).T

# 应用 E1 和 E2 变换
E1R = np.dot(R, E1.T)
E2R = np.dot(R, E2.T)

# 计算 E1E2 变换
E1E2 = np.dot(E1, E2)
E1E2R = np.dot(R, E1E2.T)

# 绘制图形
plt.figure(figsize=(8, 8))

# 绘制 R
plt.plot(R[:, 0], R[:, 1], 'b-', label='R ')

# 绘制 E1R
plt.plot(E1R[:, 0], E1R[:, 1], 'r-', label='E1R')

# 绘制 E2R
plt.plot(E2R[:, 0], E2R[:, 1], 'g-', label='E2R')

# 绘制 E1E2R
plt.plot(E1E2R[:, 0], E1E2R[:, 1], 'm-', label='E1E2R')

# 设置坐标轴和图例
plt.axhline(0, color='black', lw=0.5)
plt.axvline(0, color='black', lw=0.5)
plt.xlim(-5, 5)
plt.ylim(-5, 5)
plt.gca().set_aspect('equal', adjustable='box')
plt.title('在同一坐标系下对 R 进行 E1, E2, E1E2 变换')
plt.legend()
plt.grid()
plt.show()

7
  1. 在同一个平面坐标系中分别画出S, E1S, E2S, E1E2S;
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
import numpy as np
import matplotlib.pyplot as plt

# 设置matplotlib的字体
plt.rcParams['font.sans-serif'] = ['SimHei']  # 指定默认字体为黑体
plt.rcParams['axes.unicode_minus'] = False  # 解决保存图像时负号'-'显示为方块的

# 定义矩阵
E1 = np.array([[1, -2], [0, 2]])
E2 = np.array([[1, 0], [0, 2]])

# 定义 S 的边界，单位圆
theta = np.linspace(0, 2 * np.pi, 100)
x = np.cos(theta)
y = np.sin(theta)
S = np.array([x, y]).T

# 应用 E1
E1S = np.dot(S, E1.T)

# 应用 E2
E2S = np.dot(S, E2.T)

# 应用 E1E2
E1E2 = np.dot(E1, E2)
E1E2S = np.dot(S, E1E2.T)

# 绘制图形
plt.figure(figsize=(8, 8))

# 绘制 S
plt.plot(S[:, 0], S[:, 1], 'b-', label='S ')

# 绘制 E1S
plt.plot(E1S[:, 0], E1S[:, 1], 'r-', label='E1S')

# 绘制 E2S
plt.plot(E2S[:, 0], E2S[:, 1], 'g-', label='E2S')

# 绘制 E1E2S
plt.plot(E1E2S[:, 0], E1E2S[:, 1], 'm-', label='E1E2S')

# 设置坐标轴和图例
plt.axhline(0, color='black', lw=0.5)
plt.axvline(0, color='black', lw=0.5)

# 自动调整坐标轴范围
plt.xlim(-5, 5)
plt.ylim(-5, 5)
plt.gca().set_aspect('equal', adjustable='box')
plt.title('在同一坐标系下对 S 进行 E1, E2, E1E2 变换')
plt.legend()
plt.grid()
plt.show()

8
  1. 在同一个平面坐标系中分别画出R, E1E2R, S, E1E2S
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
import numpy as np
import matplotlib.pyplot as plt

# 定义矩阵
E1 = np.array([[1, -2], [0, 2]])
E2 = np.array([[1, 0], [0, 2]])

# 定义集合 R 的边界（正方形）
R_x = np.array([-1, 1, 1, -1, -1])
R_y = np.array([-1, -1, 1, 1, -1])
R = np.array([R_x, R_y]).T

# 定义集合 S 的边界（单位圆）
theta = np.linspace(0, 2 * np.pi, 100)
S_x = np.cos(theta)
S_y = np.sin(theta)
S = np.array([S_x, S_y]).T

# 计算 E1E2 变换
E1E2 = np.dot(E1, E2)

# 变换 R 和 S
E1E2R = np.dot(R, E1E2.T)
E1E2S = np.dot(S, E1E2.T)

# 绘制图形
plt.figure(figsize=(8, 8))

# 绘制 R
plt.plot(R[:, 0], R[:, 1], 'b-', label='R ')

# 绘制 E1E2R
plt.plot(E1E2R[:, 0], E1E2R[:, 1], 'r-', label='E1E2R')

# 绘制 S
plt.plot(S[:, 0], S[:, 1], 'g-', label='S ')

# 绘制 E1E2S
plt.plot(E1E2S[:, 0], E1E2S[:, 1], 'm-', label='E1E2S')

# 设置坐标轴和图例
plt.axhline(0, color='black', lw=0.5)
plt.axvline(0, color='black', lw=0.5)
plt.xlim(-5, 5)
plt.ylim(-5, 5)
plt.gca().set_aspect('equal', adjustable='box')
plt.title('在同一个平面坐标系中分别画出R, E1E2R, S, E1E2S')
plt.legend()
plt.grid()
plt.show()

9
  1. 画出一个\[S'= \left\{ \begin{pmatrix} x \\ y \end{pmatrix} \mid x^2 -xy + \frac{y^2}{2} = 1 \ \right\} \] 以及一个包含于S'的面积最大的六边形
10 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
import numpy as np
import matplotlib.pyplot as plt

# 定义单位圆
theta = np.linspace(0, 2 * np.pi, 100)
circle_x = np.cos(theta)
circle_y = np.sin(theta)

# 定义最大面积的六边形
def hexagon(radius=1):
    angles = np.linspace(0, 2 * np.pi, 7)[:-1]  # 6个顶点
    return radius * np.array([np.cos(angles), np.sin(angles)]).T

# 线性变换矩阵 E3
E3 = np.array([[1, -1/2], [0, 1/2]])

# 生成六边形顶点
hexagon_vertices = hexagon()

# 应用线性变换
circle_transformed = np.dot(np.vstack((circle_x, circle_y)).T, E3.T)
hexagon_transformed = np.dot(hexagon_vertices, E3.T)

# 绘制图形
plt.figure(figsize=(8, 8))

# 绘制单位圆
plt.subplot(1, 2, 1)
plt.plot(circle_x, circle_y, 'b-', label=' $S$')
plt.plot(np.append(hexagon_vertices[:, 0], hexagon_vertices[0, 0]),  # 确保六边形封闭
         np.append(hexagon_vertices[:, 1], hexagon_vertices[0, 1]), 'g--', label='六边形')
plt.axhline(0, color='black', lw=0.5)
plt.axvline(0, color='black', lw=0.5)
plt.xlim(-1.5, 1.5)
plt.ylim(-1.5, 1.5)
plt.gca().set_aspect('equal', adjustable='box')
plt.title('单位圆和内接六边形')
plt.legend()
plt.grid()

# 绘制经过变换后的图形
plt.subplot(1, 2, 2)
plt.plot(circle_transformed[:, 0], circle_transformed[:, 1], 'b-', label='变换后的$S$')
plt.plot(np.append(hexagon_transformed[:, 0], hexagon_transformed[0, 0]),  # 确保变换后的六边形封闭
         np.append(hexagon_transformed[:, 1], hexagon_transformed[0, 1]), 'g--', label='变换后的六边形')
plt.axhline(0, color='black', lw=0.5)
plt.axvline(0, color='black', lw=0.5)
plt.xlim(-1.5, 1.5)
plt.ylim(-1.5, 1.5)
plt.gca().set_aspect('equal', adjustable='box')
plt.title('变换后的圆和内接六边形')
plt.legend()
plt.grid()

plt.tight_layout()
plt.show()

11
大三上课程
#海洋工程建模基础
海洋工程建模基础
https://qingyun0118.github.io/2024/10/06/海洋工程建模基础/
作者
Wei Jicheng
发布于
2024年10月6日
许可协议