Numerical Methods II - Part 2

时隔一年,老师添加了一点新题,所以又来更新一下。

Numerical Integral, Matrices of Geometrical Transforms

近似计算定积分

Write an M-file for using quadrature formulas.
The name of the function be: numint

  • Input arguments: integrand (as a string),
    the endpoints of the interval (a, b),
    number of divisors (n),
    type of the quadrature (rectangle, trapezoid, simpson)

  • Output argument: the result of the integral.

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
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
% 7.1 近似计算定积分
% 积分(字符串)
% 区间 (a, b) 的端点
% 除数 (n)
% 正交类型(矩形、梯形、辛普森形)
function integralResult = numint(integrand, b, a, n, method)
arguments
% 检查integrand是否为字符串
integrand (1,1) string

% 检查a和b是否为数值,且a小于b
b (1,1) double {mustBeNumeric}
a (1,1) double {mustBeNumeric, mustBeLessThan(a, b)}

% 检查n是否为正整数
n (1,1) double {mustBeInteger, mustBePositive}

% 检查method是否为有效的字符串选项
method (1,1) string {mustBeMember(method, ["rectangle", "trapezoid", "simpson"])}
end

% 将字符串形式的被积函数转换为函数句柄
f = str2func(integrand);

% 根据所选的方法计算积分
switch method
case 'rectangle'
integralResult = rectangleMethod(f, a, b, n);
case 'trapezoid'
integralResult = trapezoidMethod(f, a, b, n);
case 'simpson'
integralResult = simpsonMethod(f, a, b, n);
otherwise
error('Unknown method. Please choose rectangle, trapezoid, or simpson.');
end

end

% 首先计算每个小区间的宽度 h,然后遍历每个小区间
% 计算每个区间中点的函数值,并将所有这些值累加起来
% 最后,将累加的和乘以小区间的宽度,得到积分的近似值
% 这个方法在处理高度振荡或非常不规则的函数时可能不够精确
% 在这种情况下,可能需要更多的区间(即更大的 n)
function result = rectangleMethod(f, a, b, n)
% 计算每个小区间的宽度
h = (b - a) / n;

% 初始化积分结果
integralSum = 0;

% 对每个小区间进行迭代
for i = 1:n
% 计算当前小区间的中点
x_mid = a + (i - 0.5) * h;

% 计算函数在中点的值并累加
integralSum = integralSum + f(x_mid);
end

% 计算最终的积分近似值
result = h * integralSum;
end

% 首先计算每个小区间的宽度 h
% 然后,我们初始化积分和为区间两端点处函数值的平均值
% (这是因为在梯形法中,区间两端的点只被计算一次)
% 接着,我们遍历每个小区间的内部点,将这些点处的函数值累加到积分和中
% 最后,将累加的和乘以小区间的宽度,得到积分的近似值
% 这种方法比矩形法更精确,尤其是在函数比较平滑的情况下
% 然而,对于高度振荡或非常不规则的函数,它仍然可能不够精确
function result = trapezoidMethod(f, a, b, n)
% 计算每个小区间的宽度
h = (b - a) / n;

% 初始化积分结果
integralSum = 0.5 * (f(a) + f(b));

% 对每个小区间的内部点进行迭代
for i = 1:(n-1)
x = a + i * h;
integralSum = integralSum + f(x);
end

% 计算最终的积分近似值
result = h * integralSum;
end

% 我们首先确保 n 是一个正的偶数。然后,我们计算每个小区间的宽度 h
% 接着,我们初始化积分和为区间两端点处函数值的和
% 在遍历每个小区间的内部点时,我们根据这些点是奇数位置还是偶数位置
% 分别乘以 4 或 2(这是辛普森法的特点)
% 最后,将累加的和乘以 h/3,得到积分的近似值
function result = simpsonMethod(f, a, b, n)
% 确保n为偶数
if mod(n, 2) ~= 0
error('n must be a positive even integer.');
end

% 计算每个小区间的宽度
h = (b - a) / n;

% 初始化积分结果
integralSum = f(a) + f(b);

% 对每个小区间的内部点进行迭代
for i = 1:n-1
x = a + i * h;
if mod(i, 2) == 0
integralSum = integralSum + 2 * f(x);
else
integralSum = integralSum + 4 * f(x);
end
end

% 计算最终的积分近似值
result = (h / 3) * integralSum;
end

原点 affine 变换

Write an M-file, what gives the matrix of any affin transform with fixed point in Origin
The name of the file: affin1

  • Input arguments: The images of points (0, 1) and (1, 0).
  • Output argument: The matrix of the transform.
  • If someone calls the function without input arguments give them an opportunity for the graphical input.
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
% 7.2 计算具有原点固定点的 affine 变换矩阵
% 输入参数:
% imageOf01 - 点 (0, 1) 变换后的图像
% imageOf10 - 点 (1, 0) 变换后的图像
% 输出参数:
% A - 变换的矩阵
function A = affin1(imageOf01, imageOf10)
if nargin == 0
% 如果没有输入参数,通过图形界面获取输入
disp('请图形化输入点 (0, 1) 和 (1, 0) 变换后的图像');
[x1, y1] = ginput(1); % 获取点 (0, 1) 变换后的图像
[x2, y2] = ginput(1); % 获取点 (1, 0) 变换后的图像
imageOf01 = [x1; y1];
imageOf10 = [x2; y2];
end

% 验证输入参数
if ~isequal(size(imageOf01), [2, 1]) || ~isequal(size(imageOf10), [2, 1])
error('输入参数的大小必须是 2x1');
end

% 计算 affine 变换矩阵
A = [imageOf10, imageOf01];

% 输出变换矩阵
disp('Affine 变换矩阵为:');
disp(A);
end
1
2
% 输入示例
A = affin1([2; 3], [4; 1])

三角变换

Write an M-file, what gives the matrix of any affin transform.
The name of the m-file: affin2

  • Input arguments: a triangle (with the coordinates of the edges) and the image of the triangle (also with the edges)
  • Output argument: The matrix of the transform.
  • Give an opportunity for graphical input (as in the previous exercise)
  • Draw a figure in both cases.
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
% 7.3 计算基于一个三角形及其映射后的三角形的仿射变换矩阵
%
% 使用方法:
% A = affin2(triangle, image_triangle)
% 其中 triangle 和 image_triangle 是包含三角形顶点坐标的 3x2 矩阵
% 如果没有给出参数,则使用图形输入
function A = affin2(triangle, image_triangle)
if nargin == 0
% 没有输入参数,使用图形输入
disp('点击以定义原始三角形的顶点:');
triangle = ginput(3);
disp('点击以定义映射后三角形的顶点:');
image_triangle = ginput(3);
end

if size(triangle, 1) ~= 3 || size(image_triangle, 1) ~= 3
error('两个三角形都必须有 3 个顶点。');
end

% 构造变换矩阵
% 添加一行 1 以处理仿射变换
T = [triangle, ones(3, 1)];
T_image = [image_triangle, ones(3, 1)];

% 求解变换矩阵
A = T_image' / T';

% 绘制原始三角形和变换后的三角形
figure;
plot([triangle(:,1); triangle(1,1)], [triangle(:,2); triangle(1,2)], 'b-o');
hold on;
plot([image_triangle(:,1); image_triangle(1,1)], [image_triangle(:,2); image_triangle(1,2)], 'r-o');
legend('原始三角形', '变换后的三角形');
title('三角形的仿射变换');
hold off;
end
1
2
3
4
% 输入示例
triangle = [0, 0; 1, 0; 0, 1];
image_triangle = [2, 3; 4, 1; 5, 5];
A = affin2(triangle, image_triangle)

等边三角形

Give 3 points in the plane. Two of them, P (2; 3) and Q (4; 2), lie on different sides of an equilateral triangle.
The third point, S (3; 3), is the centroid of the triangle.
Use MATLAB to find the vertices of a triangle (by adding the coordinates of the vertices) and make an illustration for the exercise.

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
% 定义点 P, Q 和 S
P = [2, 3];
Q = [4, 2];
S = [3, 3];

% 计算三角形的边长
side_length = sqrt((Q(1) - P(1))^2 + (Q(2) - P(2))^2);

% 计算外接圆的半径
R = side_length / sqrt(3);

% 计算三角形的顶点
% 我们已经有两个顶点(P 和 Q),我们需要找到第三个顶点(R)

% 计算 PQ 的中点
mid_PQ = (P + Q) / 2;

% 计算中点到质心的向量
vec_mid_to_centroid = S - mid_PQ;

% 计算 PQ 的垂直向量
perpendicular = [vec_mid_to_centroid(2), -vec_mid_to_centroid(1)];

% 归一化垂直向量
perpendicular = perpendicular / norm(perpendicular);

% 计算第三个顶点
R_vertex = S + R * perpendicular;

% 检查 R_vertex 是否与 S 在 PQ 的同一侧
% 如果不是,反转方向
if dot(R_vertex - mid_PQ, vec_mid_to_centroid) < 0
R_vertex = S - R * perpendicular;
end

% 绘制三角形
figure;
hold on;
grid on;
axis equal;

% 绘制点
plot(P(1), P(2), 'ro');
plot(Q(1), Q(2), 'ro');
plot(S(1), S(2), 'bo');
plot(R_vertex(1), R_vertex(2), 'go');

% 绘制三角形的边
plot([P(1), Q(1)], [P(2), Q(2)], 'r');
plot([P(1), R_vertex(1)], [P(2), R_vertex(2)], 'r');
plot([Q(1), R_vertex(1)], [Q(2), R_vertex(2)], 'r');

% 添加注释
text(P(1), P(2), ' P');
text(Q(1), Q(2), ' Q');
text(S(1), S(2), ' S (质心)');
text(R_vertex(1), R_vertex(2), ' R');

title('具有顶点 P, Q 和 R 的等边三角形');
xlabel('X 轴');
ylabel('Y 轴');
hold off;
Author

Aloento

Posted on

2023-12-03

Updated on

2024-03-01

Licensed under

CC BY-NC-SA 4.0

Your browser is out-of-date!

Update your browser to view this website correctly.&npsb;Update my browser now

×