问MATLAB:计算机科学的数学模型

EN

// The five points, assumed all distinct.
A = [1 0 0]';
B = [0 1 0]';
C = [0 0  1]';

D1 = [0 0 0]';
D2 = [0.5 0.5 0.5]';

// Not necessary to form the matrix M. Used for convenient plotting.
M = [A B C];

hold on
grid on

// Plot all the points and the line.
for i = 1:3
    plot3(M(1,i), M(2,i), M(3,i), '.', 'MarkerSize', 24)
end
plot3(D1(1), D1(2), D1(3), 'r.', 'MarkerSize', 24)
plot3(D2(1), D2(2), D2(3), 'r.', 'MarkerSize', 24)
plot3([D1(1) D2(1)], [D1(2) D2(2)], [D1(3) D2(3)], 'k', 'LineWidth', 2)
view(120, 20)

// Cross product of B-A and C-A
cp = cross(B-A, C-A);

// Are A, B, and C collinear?
if norm(cp) == 0
    disp('The points are collinear')
else
    disp('The points are not collinear')
    // In this case, find the plane defined by the three points.
    syms x  y z
    plane_eq = [x y z]*cp;
    offset = subs(plane_eq, {x, y, z}, {A(1), A(2), A(3)});
    disp(['The equation of the plane is ' char(plane_eq) ' = ' num2str(offset)])
    u = get(gca, 'XLim');
    v = get(gca, 'YLim');
    h = ezmesh(char(solve(plane_eq - offset, z)), [u v]);
    set(h, 'EdgeColor', [0 0 0], 'FaceColor', 'none')
    if norm((D1-A)'*cp) ~= 0
        disp('D1 does not lie in the plane.')
    end

    if norm((D2-A)'*cp) ~= 0
        disp('D2 does not lie in the plane.')
    end

    M = [D1-D2 B-A C-A];
    b = D1 - A;
    t = inv(M)*b;
    ip = D1 + t(1)*(D2 - D1);
    disp(['The intersection point is [' num2str(ip') '].'])
    plot3(ip(1), ip(2), ip(3), 'g.', 'MarkerSize', 24)

end