|
ru.algorithms
- RU.ALGORITHMS ----------------------------------------------------------------
From : Nick Gorev 2:5020/400 25 Jan 2002 07:06:42
To : Nick Gorev
Subject : Реализация с помощью метода Штурма
--------------------------------------------------------------------------------
Hello, All!
Меня попросили подробнее рассказать решение.
Hачну с постановки задачи:
Имеется тор, стандартно вложенный в R^3:
Окружность радиуса r c центром в точке (0,p,0),
лежащая в плоскости XOZ вращается вдоль окружности
с центром в точке (0,0,0) радиуса p, лежащей
в плоскости XOY. Получается тор (поверхность бублика).
Алгебраическое уравнение этой поверхности:
(x^2 + y^2)^2 + 2*(z^2 - p^2 - r^2)*(x^2 + y^2) +
(z^2 + p^2 - r^2)^2 = 0
Имеется отрезок с концами (x0,y0,z0) и (x1,y1,z1).
Требуется определить, пересекаются ли отрезок с тором.
Параметрическое уравнение отрезка
x(t) = x0 + t * vx
y(t) = y0 + t * vy
z(t) = z0 + t * vz
где
vx = x1 - x0
vy = y1 - y0
vz = z1 - z0
параметр t пробегает отрезок [0..1].
Подставляем выражения для x(t), y(t), z(t) в уравнение
тора и получаем алгебраическое уравнение 4-й степени
относительно переменной t.
f(t) = a0 + a1*t + a2*t^2 + a3*t^3 + a4*t^4 = 0
где
a0 = (x0^2 + y0^2 + z0^2 - 2*r^2)*(x0^2 + y0^2 + z0^2) +
2*p^2*(z0^2 - y0^2 - x0^2) + (p^2 + r^2)^2
a1 = 4*p^2*(z0*vz - x0*vx - y0*vy) +
4*(x0^2 + y0^2 + z0^2 - r^2)*(x0*vx + y0*vy + z0*vz)
a2 = 2*(x0^2 + y0^2 + z0^2 - r^2)*(vx^2 + vy^2 + vz^2) +
4*(x0*vx + y0*vy + z0*vz)*(x0*vx + y0*vy + z0*vz) +
2*p^2*(vz^2 - vx^2 - vy^2)
a3 = 4*(vx^2 + vy^2 + vz^2)*(x0*vx + y0*vy + z0*vz)
a4 = (vx^2 + vy^2 + vz^2)^2
Случай когда отрезок имеет нулевую длину можно сразу
исключить. Поэтому можно считать, что a4=/=0.
Теперь надо узнать, есть ли у этого уравнения корни
на отрезке [0..1].
Для этого используется метод Штурма.
Hапомню его.
Строится последовательность многочленов вида:
f_0(t) = f(t) = g_0(t) * f_1(t) - f_2(t),
f_1(t) = f'(t) = g_1(t) * f_2(t) - f_3(t),
f_2(t) = g_2(t) * f_3(t) - f_4(t),
...,
где при i>1 каждый многочлен f_i есть взятый
с коэффициентом -1 остаток, получаемый при делении
f_(i-2) на f_(i-1).
Далее, N(t) есть число перемен знака в последовательности
значений этих многочленов (члены, обращающиеся в нуль, не
учитываются). Тогда число N(a) - N(b) есть число корней
между a и b, причем каждый кратный корень считается только
один раз.
В нашем случае имеем:
f_1(t) = b0 + b1*t + b2*t^2 + b3*t^3 =
f'(t) = a1 + 2*a2*t + 3*a3*t^2 + 4*a4*t^3
Остальные многочлены последовательности находим
делением многочленов. Это легко сделать вручную,
а ленивым (вроде меня) на помощь приходит функция
rem пакета Maple.
f_2(t) = c0 + c1*t + c2*t^2
f_3(t) = d0 + d1*t
f_4(t) = e0
Hу я уж не буду выписывать здесь значения
коэффициентов c0, c1,... (см. ниже листинг программы),
замечу только, что остатки от деления многочленов
могут получаться с нулевыми старшими коэффициентами.
В этом случае последовательность получается короче,
и это тоже надо учитывать.
Реализация:
/*
Функция, определяющая пересекается ли отрезок с тором
x0,y0,z0,x1,y1,z1 - координаты концов отрезка
r - радиус образующей тора
p - радиус направляющей тора
(образующая - это маленькая окружность, вращая ее вдоль
большой направляющей окружности, мы и получаем тор)
Функция возвращает ноль, если пересечения нет
*/
#include <math.h>
#define eps 1.e-20 // типа ноль
int CrossTor(double x0, double y0, double z0,
double x1, double y1, double z1,
double r, double p) {
double vx, vy, vz, lv2, xv, lx2, a0, a1, a2, a3, a4, b0,
b1, b2, b3, c0, c1, c2, d0, d1, e0,
v0[5], v1[5], w0[5], w1[5];
int s1, s0, I, N0, N1;
vx = x1 - x0;
vy = y1 - y0;
vz = z1 - z0;
lv2 = vx * vx + vy * vy + vz * vz;
if(lv2 < eps) return 0;
xv = x0 * vx + y0 * vy + z0 * vz;
lx2 = x0 * x0 + y0 * y0 + z0 * z0;
for(I = 0; I < 5; I++) {
v0[I] = 0.;
v1[I] = 0.;
}
// Hулевой многочлен последовательности Штурма - 4-й степени
a4 = lv2 * lv2;
a3 = 4. * lv2 * xv;
a2 = 2. * (lx2 - r * r) * lv2 - 2. * p * p *
(lv2 - 2. * vz * vz) + 4. * xv * xv;
a1 = -4. * p * p * (xv - 2. * z0 * vz) + 4. *
(lx2 - r * r) * xv;
a0 = (lx2 - 2. * r * r) * lx2 + (r * r - p * p) *
(r * r - p * p) - 2. * p * p * (lx2 - 2. * z0 * z0);
v0[0] = a0;
v1[0] = a0 + a1 + a2 + a3 + a4;
// Первый многочлен последовательности - производная нулевого -
// 3-й степени
b0 = a1;
b1 = 2. * a2;
b2 = 3. * a3;
b3 = 4. * a4;
v0[1] = b0;
v1[1] = b0 + b1 + b2 + b3;
// Второй многочлен - взятый со знаком минус остаток
// от деления нулевого многочлена на первый
c0 = a1 * a3 / (a4 * 16.) - a0;
c1 = a2 * a3 / (a4 * 8.) - 0.75 * a1;
c2 = 0.1875 * a3 * a3 / a4 - 0.5 * a2;
v0[2] = c0;
v1[2] = c0 + c1 + c2;
if(fabs(c2) <= eps) {
if(fabs(c1) > eps) {
d0 = c0 * b1 / c1 - c0 * c0 * b2 / (c1 * c1)
+ c0 * c0 * c0 * b3 / (c1 * c1 * c1) - b0;
v0[3] = d0;
v1[3] = d0;
}
} else {
d0 = c0 * b2 / c2 - b0 - c0 * b3 * c1 / (c2 * c2);
d1 = b3 * c0 / c2 - b1 + c1 * b2 / c2 -
b3 * c1 * c1 / (c2 * c2);
v0[3] = d0;
v1[3] = d0 + d1;
if(fabs(d1) > eps) {
e0 = d0 * c1 / d1 - c0 - d0 * d0 * c2 / (d1 * d1);
v0[4] = e0;
v1[4] = e0;
}
}
// Избавляемся от нулей в последовательности Штурма
s0 = 0;
s1 = 0;
for(I = 0; I < 5; I++) {
if(v0[I] != 0.) {w0[s0] = v0[I]; s0++;}
if(v1[I] != 0.) {w1[s1] = v1[I]; s1++;}
}
// Вычисление числа перемен знака
N0 = 0;
N1 = 0;
for(I = 1; I <= s0; I++)
if(w0[I] * w0[I - 1] < 0.) N0++;
for(I = 1; I <= s1; I++)
if(w1[I] * w1[I - 1] < 0.) N1++;
// Возвращаем число корней
return N0 - N1;
}
With best regards, Nick Gorev. E-mail: NickGorev@mtu-net.ru
--- ifmail v.2.15dev5
* Origin: MTU-Intel ISP (2:5020/400)
Вернуться к списку тем, сортированных по: возрастание даты уменьшение даты тема автор
|
