Подробно расписана, Ваша система для
и
И это получается ровно та же система, которую я получал в своем методе (что и неудивительно).
Код:
function Main()
clear all; close all;
syms x y x0 y0 x1 y1 x2 y2 x3 y3 x4 y4 X1 Y1 X2 Y2 X3 Y3 X4 Y4
% через определители матриц
A = [2*x0, 1, 0, 0;
x0^2+y0^2, x0, y0, 1;
X1^2+Y1^2, X1, Y1, 1;
X2^2+Y2^2, X2, Y2, 1];
B = [2*y0, 0, 1, 0;
x0^2+y0^2, x0, y0, 1;
X3^2+Y3^2, X3, Y3, 1;
X4^2+Y4^2, X4, Y4, 1];
C = [2*x0, 1, 0, 0;
x0^2+y0^2, x0, y0, 1;
X3^2+Y3^2, X3, Y3, 1;
X4^2+Y4^2, X4, Y4, 1];
D = [2*y0, 0, 1, 0;
x0^2+y0^2, x0, y0, 1;
X1^2+Y1^2, X1, Y1, 1;
X2^2+Y2^2, X2, Y2, 1];
g = det(A)*det(B) - det(C)*det(D);
% напрямую через уравнения прямых
a(x) = x0 - x;
b(y) = y0 - y;
c(x, y) = (x0 - x)*x + (y0 - y)*y;
zp(x1, y1, x2, y2) = a(x1)*b(y2) - a(x2)*b(y1);
xp(x1, y1, x2, y2) = (b(y2)*c(x1, y1) - b(y1)*c(x2, y2));
yp(x1, y1, x2, y2) = (a(x1)*c(x2, y2) - a(x2)*c(x1, y1));
zp12 = zp(X1, Y1, X2, Y2);
zp34 = zp(X3, Y3, X4, Y4);
xp12 = xp(X1, Y1, X2, Y2);
xp34 = xp(X3, Y3, X4, Y4);
yp12 = yp(X1, Y1, X2, Y2);
yp34 = yp(X3, Y3, X4, Y4);
f = y0*(xp34*zp12 - xp12*zp34) - x0*(yp34*zp12 - yp12*zp34) - yp12*xp34 + xp12*yp34;
% сравнение результатов
factor(g)
factor(f)
simplify(f - g)
end