课程大实验:椭球面上两点间的最短距离

最新推荐文章于 2024-05-14 18:01:31 发布

一个毛毛虫

最新推荐文章于 2024-05-14 18:01:31 发布

阅读量627

点赞数 9

分类专栏：电子科技大学数学实验练习题文章标签： matlab 数学建模

本文链接：https://blog.csdn.net/ksnzuxbxjnd/article/details/136121015

版权

电子科技大学数学实验练习题专栏收录该内容

6 篇文章 7 订阅

订阅专栏

课程大实验:椭球面上两点间的最短距离

一、实验内容

已知椭球面方程: $\frac{x^2}{a^2}+\frac{y^2}{a^2}+\frac{z^2}{b^2}=1$ .设a=6000,b=5000.又已知该椭球面上两点 $P_1$ (2200,3600, $z_1$ ), $P_2$ (2900,3300, $z_2$ ).请设计算法估算 $P_1$ 和 $P_2$ 二点在椭球面上的最短距离.这里 $z_1,z_2$ 均大于0.

二、实验目的

1.掌握曲面的绘制方法；

2.基于贪心策略设计算法；

3.算法设计与分析的应用。

三、实验过程

1. 问题分析

1.1实验基础

首先考虑普通球面上两点之间的最短距离，根据数学知识不难得出，球面上任意两点之间的最短距离为这两点与球心三点形成的平面与球的截面圆上所对应的劣弧。例题如下：

已知一球面方程: $\frac{x^2}{a^2}+\frac{y^2}{a^2}+\frac{z^2}{a^2}=1$ .其中, .又已知球面上2点的 (300,320, $z_1$ ), (450,800, $z_2$ ).请编程计算 $P_1$ 和 $P_2$ 二点在球面上的最短距离.这里 $z_1,z_2$ 均大于0.

1.2算法分析

贪心算法是一种解决最优化问题的算法思想。在每一步选择中，贪心算法总是选择当前状态下最优的选择，而不考虑全局的最优解。它通过贪心的方式进行局部最优的选择，希望最终能够得到全局最优解。

贪心算法的基本思路是，每次选择当前状态下的最优解，并将其添加到解集中，然后更新问题的状态，继续进行下一步选择，直到达到问题的结束条件。

因此在应用贪心算法来计算近似解的基本思想是，将一个观测点从 $P_1$ 开始，逐步挪向 $P_2$ 在每一次挪动方向的选择上，我们都基于“贪心”思想，选择挪动后与之间的距离减少最多的方向，每一步挪动只考虑当前的最优解，而不考虑全局的最优解，期望通过每一步的最优解来得到全局的最优解。

2.模型假设

挪动步长足够小，观测点所在椭球面可近似为平面；
观测点每次移动的距离近似于对应椭球表面曲线弧长；
当观测点与 $P_2$ 于约定的固定值时 $\varepsilon$ ，就近似为到达 $P_2$ 。

3.变量与符号说明

变量	符号说明
M	观测点
$N_i$	观测点下一步所能到达的椭球面所有点的集合
$N_M$	$N_i$ 中使 $\vec{n}·\vec{MN_i}$ 值最大的点,即局部最优解
$d_i$	观测点每次移动的直线距离
D	观测点总共移动的距离
$\varepsilon$	模型假设2.3中的约定固定值
$\vec{n}$	$P_2$ 处椭球面向外的法向量

4.模型建立与算法设计

4.1模型建立

根据椭球面的参数方程，有：

$\begin{cases}x=acos\theta cos\alpha\\ y=acos\theta sin\alpha\\ z=bsin\theta \end{cases}$

利用MATLAB画出椭球面并标注出 $P_1和P_2$ 两点。
图1 椭球面

4.2算法设计

建立直角坐标系,计算出 $P_2$ 的外法向量 $\vec{n}$
设观测点为M(第一步时为 $P_1$ ), $N i$ 是在椭球面上以M为圆心,半径足够小的圆上的点的集合,计算出 $\vec{MN_i}$
计算 $\vec{n}·\vec{MN_i}$ 的值,并找到使该值最大的向量 $\vec{MN_M}$ ,该向量方向就是当前状态下的最优解.
计算当前M与 $N_M$ 的距离,该距离就是观测点此次移动的直线距离 $d_i$ ,用其近似代替曲线弧长
以 $N_M$ 为新的观测点A,重复2-4过程,知道 $d_i$ 小于固定值 $\varepsilon$
计算观测点总共移动的距离D,近似为两点之间的最短距离.

5.运行结果及结果分析

5.1 运行结果

绘制出椭球面以及面上两点之间的最短路径，如图所示

图2 最短路径
图3 挪动路径
图4 运行结果

5.2运行结果分析

通过实验结果可知,根据贪心算法,经过77次迭代后,求得 $P_1与P_2$ 之间的最短路径长度近似为783.457350.

相对来说 $P_1与P_2$ 间隔不远,可用两点之间直线距离来评估我们结果的合理性.

利用两点之间距离公式的 $P_1与P_2$ 之间距离为776.1416;我们所得计算结果略大于此值,从这个角度看,结果比较合理.

四、优缺点及改进方向

优点：算法简单，利用贪心算法，每一步仅需考虑当前最优解，而不用关注全局或过去的状态，运行结果也较为温和。

缺点：计算量较大，对于椭球面上两点间的最短距离，贪心算法并不是一个最优的解决方法。贪心算法通常在每一步选择当前最优解，然后希望通过这种局部最优解能够达到全局最优解。但在椭球面上计算最短距离涉及复杂的几何计算，简单地通过贪心策略来求解可能出错。

改进方向：Vincenty 公式是专门为计算椭球面上两点间的最短距离而设计的算法，它考虑了椭球体参数和大圆路径上的各种因素，能够得到较为精确的结果。因此，使用 Vincenty 公式或其他类似的算法是求解椭球面上两点间最短距离的可靠方法。

五、心得体会与总结

通过本次实验，将学习过的贪心算法应用到具体的实际问题中，数学建模和分析解决能力得到进一步提升。一方面，我们对MATLAB的语法及常用函数和功能特别是绘图功能有了更好的掌握；另一方面，我们对微积分和线性代数中的知识如空间立体、梯度等有了更深的理解。整体上提高了用知识解决实际问题的工程能力。

六、对本实验问题的设计提出改进意见

无

附件

附件1.主程序:main.m

a = 6000;
b = 5000;
ellipsoid_axes = [a a b];
P1 = [2200 3600];
P2 = [2900 3300];
goal = 5;
total_distance = 0;
total_steps = 1;
P1 = [P1 getEllipsoidZ(ellipsoid_axes, P1)];
P2 = [P2 getEllipsoidZ(ellipsoid_axes, P2)];
P2_normal = getEllipsoidNormal(ellipsoid_axes, P2);

cur_pos = P1;
nexts = [];

while (true)
    next_pos = getNextPos(ellipsoid_axes, cur_pos, P2_normal);
    nexts = [nexts; next_pos];
    walked_distance = getDistance(cur_pos, next_pos);
    total_distance = total_distance + walked_distance;
    dst_distance = getDistance(next_pos, P2);

    if dst_distance < goal
        break
    else
        cur_pos = next_pos;
        total_steps = total_steps+1;
    end

end

fprintf('total_distance=%f\n',total_distance);
fprintf('total_steps=%d\n',total_steps);
fprintf('P1=');
fprintf('%f ',P1);
fprintf('\nP2=');
fprintf('%f ',P2);

[theta, alpha] = meshgrid(linspace(0, pi / 2, 50), linspace(0, 2 * pi, 50));
z = b * sin(theta);
x = a * cos(theta) .* cos(alpha);
y = a * cos(theta) .* sin(alpha);
mesh(x, y, z)
hold on;

plot3(P1(1), P1(2), P1(3), 'r.', 'markersize', 6)
plot3(P2(1), P2(2), P2(3), 'r.', 'markersize', 6)
plot3(nexts(:, 1), nexts(:, 2), nexts(:, 3), 'r.', 'markersize', 6)

附件2.距离子函数getDistance.m

function distance = getDistance(pointA, pointB)
    % Calculate the distance between A and B.
    distance = sqrt((pointB(1) - pointA(1))^2 + (pointB(2) - pointA(2))^2 + (pointB(3) - pointA(3))^2);
end

附件3. 子函数getEllipsoidNormal

function vec = getEllipsoidNormal(ellipsoid_coord, point)
    % Get the normal vector of a point on an ellipsoid.
    a = ellipsoid_coord(1);
    b = ellipsoid_coord(2);
    c = ellipsoid_coord(3);
    x = point(1);
    y = point(2);
    z = point(3);

    vec = [2 * x / a^2, 2 * y / b^2, 2 * z / c^2];
    vec = getUnitVector([0 0 0], vec);
end

附件4.子函数getEllipsoidZ

function z = getEllipsoidZ(ellipsoid_coord, point)
    % Get Z coordinate of a point on an ellipsoid.
    a = ellipsoid_coord(1);
    b = ellipsoid_coord(2);
    c = ellipsoid_coord(3);
    x = point(1);
    y = point(2);

    z = c * sqrt(1 - x^2 / a^2 - y^2 / b^2);
end

附件5.子函数getNextPos

function next_pos = getNextPos(ellipsoid_axes, cur, normal_vector)
    % Get the rectangular coordinate the next_pos step should land on.
    max_product = 0;
    next_pos = [0 0 0];

    cur_polar = rectToPolar(ellipsoid_axes, cur);
    next_polars = getNextPolars(cur_polar);

    for next_polar = next_polars
        next_rect = polarToRect(ellipsoid_axes, next_polar);
        next_unit = getUnitVector(cur, next_rect);
        product = next_unit * normal_vector';

        if product > max_product
            max_product = product;
            next_pos = next_rect;
        end

    end

end

function rect_coord = polarToRect(ellipsoid_axes, polar_coord)
    % Convert polar coordinates to rectangular coordinates.
    a = ellipsoid_axes(1);
    b = ellipsoid_axes(2);
    c = ellipsoid_axes(3);
    theta = polar_coord(1);
    alpha = polar_coord(2);

    x = a * cos(theta) * cos(alpha);
    y = b * cos(theta) * sin(alpha);
    z = c * sin(theta);

    rect_coord = [x y z];
end

function polar_coord = rectToPolar(ellipsoid_axes, rect_coord)
    % Convert rectangular coordinates to polar cordinates.
    a = ellipsoid_axes(1);
    b = ellipsoid_axes(2);
    c = ellipsoid_axes(3);
    x = rect_coord(1);
    y = rect_coord(2);
    z = rect_coord(3);

    tan_alpha = (y / b) / (x / a);
    alpha = atan(tan_alpha);

    tan_theta = sqrt((z^2 / c^2) / (x^2 / a^2 + y^2 / b^2));
    theta = atan(tan_theta);

    polar_coord = [theta alpha];
end

function next_polars = getNextPolars(cur_polar)
    % Get all ends of vectors radiating from `cur_polar`, among which the next_pos step will select one.
    radius = 0.01;
    precision = 0.001;
    num = 2 * radius / precision + 1;

    theta = linspace(cur_polar(1) - radius, cur_polar(1) + radius, num);
    alpha = linspace(cur_polar(2) - radius, cur_polar(2) + radius, num);
    [theta, alpha] = meshgrid(theta, alpha);
    theta = reshape(theta, [1 num^2]);
    alpha = reshape(alpha, [1 num^2]);

    next_polars = [theta; alpha];
end

附件6.子函数getUnitVector

function vec = getUnitVector(pointA, pointB)
    % Get the unit vector pointing from A to B.
    dx = pointB(1) - pointA(1);
    dy = pointB(2) - pointA(2);
    dz = pointB(3) - pointA(3);

    vec = [dx, dy, dz];
    vec = vec / norm(vec);
end