Линия пересечения поверхностей. Метод Драгилева.

Форум пользователей пакета Maple

Модератор: Admin

алексей_алексей
Сообщения: 1776
Зарегистрирован: Вс май 01, 2005 9:02 pm

Линия пересечения поверхностей. Метод Драгилева.

Сообщение алексей_алексей » Сб мар 30, 2013 5:36 pm

Mетод Драгилева и вопрос с MaplePrimes:
http://www.mapleprimes.com/questions/14 ... t-Function

f=sin(theta)-2*sin(2*x-theta+(1/3)*Pi)/sin(x+(1/3)*Pi-theta)=0;

Имеем функцию двух переменных, равную 0. Такое выражение между переменными соответствует множеству точек на линии. Явно мы не можем выразить одну переменную через другую, и применяем метод Драгилева, чтобы найти все точки линии. Для начала попытаемся найти хотя бы одну точку на линии. Произведём удобную для глаз замену :

theta := x1; x := x2; sin(x1)-2*cos(-2*x2+x1+(1/6)*Pi)/cos(-x2+(1/6)*Pi+x1)=0;

В чем Maple помогает нам. Если, например, x1=0, то очевидное значение второй переменной Pi/3. Приняв обе переменные за функции от некой независимой переменной s , мы дифференцируем исходную функцию:

(df/dx1)*(dx1/ds) + (df/dx2)*(dx2/ds)=0; (1)

Что в конкретном случае имеет вид:

(cos(x1)+2*sin(-2*x2+x1+(1/6)*Pi)/cos(-x2+(1/6)*Pi+x1)-2*cos(-2*x2+x1+(1/6)*Pi)*sin(-x2+(1/6)*Pi+x1)/cos(-x2+(1/6)*Pi+x1)^2)*dx1/ds+(-4*sin(-2*x2+x1+(1/6)*Pi)/cos(-x2+(1/6)*Pi+x1)+2*cos(-2*x2+x1+(1/6)*Pi)*sin(-x2+(1/6)*Pi+x1)/cos(-x2+(1/6)*Pi+x1)^2)*dx2/ds = 0;

Относительно dx1/ds и dx2/ds (1) является линейной однородной алгебраической системой уравнений. Назначаем любую из переменных независимой, присваиваем ей любое значение, и получаем следующую систему автономных дифференциальных уравнений:

dx1/ds = - df/dx2 ;
dx2/ds = df/dx1 ;


Где s на самом деле будет длиной дуги при соответствующей нормировке. Вспоминая о найденной ранее точке на линии (0,Pi/3) , берём её в качестве начальной для решения дифференциального уравнения и строим решение в обе стороны от неё. Maple подсказывает нам пределы интервала интегрирования. Решаем диффуравнение, строим график, смотрим на знаменатель второго слагаемого данной функции, проверяем построение с помощью implicitplot, убеждаемся в ошибке Maple относительно прямолинейных участков и обнаруживаем периодичность полученного решения по обеим переменным. Период по каждой из них (2n - +1)Pi, что логично и имеет подтверждение на графике implicitplot.
Текст, вроде, находится по ссылке, и на всякий случай несколько точек одного из участков линии с невязкой (не забываем о периоде):

"theta =",-2.0925005007182698, "x=", 0.0011160467088066897
"f1=",-2.205349833306336e-6
"theta =",-0.990846823103417, "x=", 0.6448981494374012
"f1=", 1.1416547462950177e-7
"theta =", -0.5685327511784426, "x=", 0.8477725044734922
"f1=", -5.344985283262105e-7
"theta =", -0.2629667293645833, "x=", 0.9653037597258018
"f1=", -3.112680203587459e-7
"theta =", -0.012471214507812373, "x=", 1.0436479689520148
"f1=", -1.4131847879406134e-9
"theta =", 0.20758702295412537, "x=", 1.1029112642162227
"f1=", 2.273146131881454e-7
"theta =", 0.4100563779714732, "x=",1.1543541345922472
"f1=", 1.6856249668295575e-7
"theta =", 0.6029672301388762, "x=", 1.2053961186356148
"f1=", 1.302299429406517e-7
"theta =", 0.792226585268965, "x=", 1.2616203705829965
"f1=", 2.0308207171471082e-7
"theta =", 0.9828717674314443, "x=", 1.327725387189918
"f1=", 2.4272807730429946e-7
"theta =",1.1797129339720567, "x=", 1.408006884872145
"f1=", -9.972637131649975e-7
"theta =",1.3876441924765175, "x=",1.5065581954110674
"f1=", 6.29553041586739e-7
"theta =",1.611690876740141, "x=",1.6272048452050252
"f1=", -1.4371806885682403e-6
"theta =", 1.8567985568489118, "x=",1.7731447015165294
"f1=", -1.3781162930825985e-6
"theta =", 2.127373694060153, "x=",1.9462246402649659
"f1=", 6.63589013627508e-7
"theta =", 2.426861620626197, "x=",2.1459952911600784
"f1=", 2.5939935393015645e-6
"theta =", 2.7583692333313397, "x=",2.369098062557234
"f1=", 5.333817258423856e-6
"theta =", 3.1295744557588754, "x=",(2.610471426280491
"f1=", 4.528650775421569e-6
"theta =", 3.580135852720373, "x=", 2.872670272009993
"f1=", 7.649256314967712e-6
Последний раз редактировалось алексей_алексей Пт май 10, 2013 9:18 pm, всего редактировалось 2 раза.

алексей_алексей
Сообщения: 1776
Зарегистрирован: Вс май 01, 2005 9:02 pm

Re: Метод Драгилева

Сообщение алексей_алексей » Пт май 10, 2013 9:16 pm

Оказывается, по поиску “линия пересечения поверхностей” народ может набрести на несусветную муть, не смотря на наличие публикаций по методу Драгилева. Поэтому решил дать соответствующее название теме. Здесь постараюсь привести примеры, которые можно непосредственно называть нахождением линии пересечения поверхностей без таких слов, как: “анимация пространственных рычажных механизмов” или “расчёт траекторий всех его точек”, “решение системы нелинейных уравнений”, “неявные функции”. Что поделаешь, ищущие линию пересечения поверхностей, как правило, не в состоянии увидеть другой её ипостаси. Местность давно огорожена красными флажками пособий, книг и методичек от разного рода технических учёных, начавших немерено плодиться ещё в советское время и продолжающих это увлекательное занятие, похоже, как основное. Подумать только, через флажки не может перескочить даже Google. Ничего, будем расчищать дорогу прогрессу…

алексей_алексей
Сообщения: 1776
Зарегистрирован: Вс май 01, 2005 9:02 pm

Линия пересечения поверхностей. Метод Драгилева.

Сообщение алексей_алексей » Сб май 11, 2013 2:31 pm

Старый пример от Пахоменкова Ю. М. Пример из жизни. Мы договорились, что это не решение СНУ. Три поверхности в 4d образуют линию.
(4*(1-2*cos(x1)+2*cos(x2)-2*cos(x3)))/Pi-x4 = 0;
(4*(1-2*cos(5*x1)+2*cos(5*x2)-2*cos(5*x3)))/(5*Pi) = 0;
(4*(1-2*cos(7*x1)+2*cos(7*x2)-2*cos(7*x3)))/(7*Pi) = 0;

Посмотрим, каким образом эта линия может быть найдена. Придадим уравнениям алгебраический вид:

Код: Выделить всё

restart:
f1 := (4*(1-2*cos(x1)+2*cos(x2)-2*cos(x3)))/Pi-x4 = 0;
f2 := (4*(1-2*cos(5*x1)+2*cos(5*x2)-2*cos(5*x3)))/(5*Pi) = 0;
f3 := (4*(1-2*cos(7*x1)+2*cos(7*x2)-2*cos(7*x3)))/(7*Pi) = 0;
f1 := subs(cos(x1) = x1, f1):
f1 := subs(cos(x2) = x2, f1):
f1 := expand(Pi*subs(cos(x3) = x3, f1));
f2 := expand(f2):
f2 := subs(cos(x1) = x1, f2):
f2 := subs(cos(x2) = x2, f2):
f2 := expand(Pi*subs(cos(x3) = x3, f2));
f3 := expand(f3):
f3 := subs(cos(x1) = x1, f3):
f3 := subs(cos(x2) = x2, f3):
f3 := expand(Pi*subs(cos(x3) = x3, f3));

Чтобы воспользоваться процедурой Isolate

Код: Выделить всё

restart:
with(RootFinding):
f1 := 4-8*x1+8*x2-8*x3-evalf(Pi)*x4;
f2 := 4/5-(128/5)*x1^5+32*x1^3-8*x1+(128/5)*x2^5-32*x2^3+8*x2-(128/5)*x3^5+32*x3^3-8*x3;
f3 := 4/7-(512/7)*x1^7+128*x1^5-64*x1^3+8*x1+(512/7)*x2^7-128*x2^5+64*x2^3-8*x2-(512/7)*x3^7+128*x3^5-64*x3^3+8*x3;
f4 := x1+x2+x3+x4;
Isolate([f1, f2, f3, f4], [x1, x2, x3, x4]);

Добавив к существующим уравнениям уравнение, например, плоскости f4, мы попытаемся отыскать отдельные точки на линии...

Продолжение следует. Пока оно не последовало, можно предложить что-нибудь лучше.

алексей_алексей
Сообщения: 1776
Зарегистрирован: Вс май 01, 2005 9:02 pm

Re: Линия пересечения поверхностей. Метод Драгилева.

Сообщение алексей_алексей » Сб май 11, 2013 3:49 pm

Получив решения алгебраической системы, мы теперь можем применить метод Драгилева для вычисления линии пересечения поверхностей. Решения нам подходят не все, а те, где x1,x2,x3 принимают значения по модулю <= 1, дабы избежать непонятных и ненужных в данном случае комплексных вычислений. Эти решения берутся в качестве начальной точки задачи Коши, и от этих точек строятся участки линии 4d. Работа частично ручная, потому что не очень понятно, как автоматизировать процесс. Множество точек линии, скорее всего, не является связным, и мы строим только отдельные связные участки. Возможно попадание нескольких начальных точек на один участок. Проверку можно осуществлять графически, строя проекции линии на 3d.
Линия строится на основе исходных уравнений, то есть, до замены их на алгебраические. И помним, что x4 осталась без изменений, а остальные переменные поменяют свои значения соответствующим образом. Для x1,x2,x3 решение с периодом 2Pi в любом направлении и в любых сочетаниях, но для x4 периодичности нет.
Текст вычисления участка линии 4d взят для точки на линии x1 = .5, x2 = -.25, x3 = -.25, x4 = 0. И участок вычислен для диапазона параметра t от -0,5 до +0,5. Участок в моём понимании какой-то особенный, потому что x1 и x4 на нём остаются неизменными. На печать выводится значение параметра, значения переменных и все невязки.
Можно брать любые другие найденные начальные точки с учётом периодичности первых трёх переменных и при этом смотреть на невязки…

Код: Выделить всё

restart:
with(LinearAlgebra):
smin := -.5; smax := .5; h := 0.1e-1:
f1 := (4*(1-2*cos(x1)+2*cos(x2)-2*cos(x3)))/evalf(Pi)-x4:
f2 := (4*(1-2*cos(5*x1)+2*cos(5*x2)-2*cos(5*x3)))/(5*evalf(Pi)):
f3 := (4*(1-2*cos(7*x1)+2*cos(7*x2)-2*cos(7*x3)))/(7*evalf(Pi)):
x01 := arccos(.5);
x02 := arccos(-.25);
x03 := arccos(-.25);
x04 := -0.;
n := 3:
x := seq(eval(cat('x', i)), i = 1 .. n+1):
x0 := seq(eval(cat('x0', i)), i = 1 .. n+1):
F := [seq(unapply(eval(cat('f', i)), [x]), i = 1 .. n)]:
A := MTM:-jacobian(F(x), [x]):
A := MTM:-subs(A, [x], [x(s)]):
deqs := seq(diff(x[i](s), s) = simplify(Determinant(DeleteColumn(ColumnOperation(A, [i, n+1]), n+1))), i = 1 .. n), diff(x[n+1](s), s) = simplify(-Determinant(DeleteColumn(A, n+1))); ics := seq(x[i](0) = x0[i], i = 1 .. n+1):
sln := dsolve([deqs, ics], numeric, method = rkf45, abserr = (1/100000)*h, maxfun = 100000, range = smin .. smax); B := `~`[rhs](Matrix(`~`[sln]([seq(i, i = smin .. smax, h)]))[() .. (), 1 .. n+2]):
m := ArrayNumElems(B[() .. (), 1]): m;
k := 1; for j to n+1 do L[j] := B[() .. (), k+j] end do:
for i from 0 to m-1 do "t =", op(convert(map(rhs, sln(smin+h*i)[1 .. 1]), list)), "x =", op(convert(map(rhs, sln(smin+h*i)[2 .. n+2]), list)):
for j to n do print(f, j, "=", F[j](op(convert(map(rhs, sln(smin+h*i)[2 .. n+2]), list)))) end do end do

алексей_алексей
Сообщения: 1776
Зарегистрирован: Вс май 01, 2005 9:02 pm

Re: Линия пересечения поверхностей. Метод Драгилева.

Сообщение алексей_алексей » Вс май 12, 2013 1:47 pm

x1^2+x2^2-x3^2=0;
x1+2*x2+2*x3-1=0;

Кто помнит и знает, метод Драгилева способен вычислять линию пересечения поверхностей не только численно, но и в любимой многими замкнутой форме или в буквенном виде. И вот тому пример. Не знаю, как Ньютон со товарищи не путались в траекториях космических тел, отличая эллипс от параболы или гиперболы, но это, кажется, парабола. Да здесь это не суть важно, главное, что есть уравнение плоской кривой в пространстве…

Код: Выделить всё

restart:
with(LinearAlgebra):
with(plots):
f1 := x1^2+x2^2-x3^2;
f2 := x1+2*x2+2*x3-1;
x01 := -1;
x02 := 0.;
x03 := 1;
n := 2:
x := seq(eval(cat('x', i)), i = 1 .. n+1):
x0 := seq(eval(cat('x0', i)), i = 1 .. n+1):
F := [seq(unapply(eval(cat('f', i)), [x]), i = 1 .. n)]:
A := MTM:-jacobian(F(x), [x]):
A := MTM:-subs(A, [x], [x(s)]):
deqs := seq(diff(x[i](s), s) = simplify(Determinant(DeleteColumn(ColumnOperation(A, [i, n+1]), n+1))), i = 1 .. n), diff(x[n+1](s), s) = simplify(-Determinant(DeleteColumn(A, n+1))):
ics := seq(x[i](0) = x0[i], i = 1 .. n+1):
sln := dsolve([deqs, ics]):
"x1(s)=", rhs(sln[1]);
"x2(s)=", rhs(sln[2]);
"x3(s)=", rhs(sln[3]);
S1 := spacecurve({[rhs(sln[1]), rhs(sln[2]), rhs(sln[3])]}, s = -1.25 .. .45, axes = normal, color = red, thickness = 10, numpoints = 5000):
S2 := plot3d([(x1^2+x2^2)^.5, -(x1+2*x2-1)/(2.)], x1 = -2.5 .. 2.5, x2 = -2.5 .. 2.5, color = [blue, green], transparency = .6):
display(S1, S2);

алексей_алексей
Сообщения: 1776
Зарегистрирован: Вс май 01, 2005 9:02 pm

Re: Линия пересечения поверхностей. Метод Драгилева.

Сообщение алексей_алексей » Вс май 12, 2013 2:28 pm

Изображение
x1^2+x2^2-x3^2=0;
x1-2*x3+5=0;


x1(s)= -(1666666667/500000000)*cos(2*sqrt(3)*s)+1666666667/1000000000;
x2(s)= (1666666667/1000000000)*sin(2*sqrt(3)*s)*sqrt(3);
x3(s)= -(1666666667/1000000000)*cos(2*sqrt(3)*s)+1666666667/500000000;

Линия пересечения двух поверхностей в картинке и в буквенном виде.
Примерно так и круче лётают к нам кометы…

алексей_алексей
Сообщения: 1776
Зарегистрирован: Вс май 01, 2005 9:02 pm

Re: Линия пересечения поверхностей. Метод Драгилева.

Сообщение алексей_алексей » Пн май 13, 2013 2:39 pm

Наверное, это будет последний пример из негромоздких и визуализируемых. Тот же самый текст, но решение системы диффуравнений в буквах не получается, зато очень хорошо получается численно. (Пример из темы про анимацию механизмов, где эта же самая линия пересечения поверхностей служила направляющей для ползуна.) Всё стандартно: находим любым доступным способом одну точку на линии и решаем задачу Коши…

(x1-2)^2+(x2-2)^2+(x3+.5*sin(10*x1))^2-9=0;
x1^6+x2^6+x3^6-12=0;


Изображение

алексей_алексей
Сообщения: 1776
Зарегистрирован: Вс май 01, 2005 9:02 pm

Re: Линия пересечения поверхностей. Метод Драгилева.

Сообщение алексей_алексей » Ср май 15, 2013 7:43 pm

Код: Выделить всё

restart:
with(LinearAlgebra):
smin := 0: smax := 0.2e-4: h := 0.1e-5:
CD1 := -1.; CD2 := 2; CD3 := 1.; g1 := 5.; g2 := 3.; g3 := 2.;
L1 := 2.4; L2 := 7.2; L3 := 2.4; L4 := 1.; L5 := 5.6; L6 := 4.5;
f1 := (CD1-x4)^2+(CD2-x5)^2+(CD3-x6)^2-L1^2;
f2 := x1-5;
f3 := x4+.5*x5;
f4 := (g1-x1)^2+(g2-x2)^2+(g3-x3)^2-L3^2;
f5 := (x4-x1)^2+(x5-x2)^2+(x6-x3)^2-L2^2;
f6 := (x7-2.)^2+(x9-3.)^2-.25*(x8+1.5)^2;
f7 := (x4-x7)^2+(x5-x8)^2+(x6-x9)^2-L5^2;
f8 := (x1-x7)^2+(x2-x8)^2+(x3-x9)^2-L6^2;
x01 := 5.;
x02 := 4.326649916;
x03 := 0.;
x04 := -.6671120731;
x05 := 1.334224146;
x06 := 3.281650311;
x07 := 4.240771345;
x08 := 3.787034644;
x09 := 4.402543266;
n := 8;
x := seq(eval(cat('x', i)), i = 1 .. n+1):
x0 := seq(eval(cat('x0', i)), i = 1 .. n+1):
F := [seq(unapply(eval(cat('f', i)), [x]), i = 1 .. n)]:
A := MTM:-jacobian(F(x), [x]):
A := MTM:-subs(A, [x], [x(s)]):
deqs := seq(diff(x[i](s), s) = simplify(Determinant(DeleteColumn(ColumnOperation(A, [i, n+1]), n+1))), i = 1 .. n), diff(x[n+1](s), s) = simplify(-Determinant(DeleteColumn(A, n+1))):
ics := seq(x[i](0) = x0[i], i = 1 .. n+1):
sln := dsolve([deqs, ics], numeric, method = rkf45, abserr = (1/1000)*h, maxfun = 10000, range = smin .. smax):
B := `~`[rhs](Matrix(`~`[sln]([seq(i, i = smin .. smax, h)]))[() .. (), 1 .. n+2]):
m := ArrayNumElems(B[() .. (), 1]);
k := 1:
for j to n+1 do L[j] := B[() .. (), k+j] end do:
for i from 0 to m-1 do "t =", op(convert(map(rhs, sln(smin+h*i)[1 .. 1]), list)), "x =", op(convert(map(rhs, sln(smin+h*i)[2 .. n+2]), list)); for j to n do print(f, j, "=", F[j](op(convert(map(rhs, sln(smin+h*i)[2 .. n+2]), list)))) end do end do;


Кусочек линии из пространства 9d в количестве 21 точки (можно больше). Хотел, было, в 12d, но Maple, как уже сообщалось в одной из тем, чередует в алфавитном порядке переменные x1,…,x9,x10,x11,x12 так, что x10,x11,x12 становятся сразу вслед за x1, а x9 занимает законное последнее место. Справиться могу только руками, чем и занимался в механизмах, но здесь у нас цикл с координатами точек линии и невязками каждого из уравнений, а как справиться с такой чехардой в цикле, не знаю. Поэтому 9d. От предыдущих текстов отличается ничем, точнее, значением “n”.
Называется, напоминаю, линия пересечения поверхностей. Интересно, какие имеются альтернативы?...

Markiyan Hirnyk
Сообщения: 1366
Зарегистрирован: Вс дек 04, 2011 11:07 pm

Сила!

Сообщение Markiyan Hirnyk » Ср май 15, 2013 8:49 pm

Результат команды
plots:-intersectplot((x1-2)^2+(x2-2)^2+(x3+.5*sin(10*x1))^2-9 = 0, x1^6+x2^6+x3^6-12 = 0, x1 = -1.5 .. 1.6, x2 = -2 .. 2, x3 = -2 .. 2, numpoints = 10000000)
и в подметки не годится.

алексей_алексей
Сообщения: 1776
Зарегистрирован: Вс май 01, 2005 9:02 pm

Re: Сила!

Сообщение алексей_алексей » Ср май 15, 2013 9:24 pm

Markiyan Hirnyk писал(а):Результат команды
plots:-intersectplot((x1-2)^2+(x2-2)^2+(x3+.5*sin(10*x1))^2-9 = 0, x1^6+x2^6+x3^6-12 = 0, x1 = -1.5 .. 1.6, x2 = -2 .. 2, x3 = -2 .. 2, numpoints = 10000000)
и в подметки не годится.

Нет, в данном случае всё не совсем так. Эта процедура просто не реагирует на numpoints.
plots:-intersectplot((x1-2)^2+(x2-2)^2+(x3+.5*sin(10*x1))^2-9 = 0, x1^6+x2^6+x3^6-12 = 0, x1 = -1.5 .. 1.6, x2 = -2 .. 2, x3 = -2 .. 2, thickness = 6, grid = [50, 50, 50]);
Формально вид будет ничего так, если не говорить о скорости, точности и о возможности анимации. Вот примеры 12d (из механизмов) и здешний 9d уже вряд ли могут быть исполнены, тем более, если не будут алгебраическими…

Markiyan Hirnyk
Сообщения: 1366
Зарегистрирован: Вс дек 04, 2011 11:07 pm

Сообщение Markiyan Hirnyk » Ср май 15, 2013 9:37 pm

Спасибо. Ну Вы поднаторели в Maple.

алексей_алексей
Сообщения: 1776
Зарегистрирован: Вс май 01, 2005 9:02 pm

Сообщение алексей_алексей » Чт май 16, 2013 8:32 am

Да, уж...

Markiyan Hirnyk
Сообщения: 1366
Зарегистрирован: Вс дек 04, 2011 11:07 pm

Прямо-таки

Сообщение Markiyan Hirnyk » Пт май 17, 2013 9:36 pm

алексей_алексей писал(а):Да уж...

Прямо-таки.

алексей_алексей
Сообщения: 1776
Зарегистрирован: Вс май 01, 2005 9:02 pm

Re: Прямо-таки

Сообщение алексей_алексей » Сб май 18, 2013 9:09 am

Markiyan Hirnyk писал(а):
алексей_алексей писал(а):Да уж...

Прямо-таки.

Чтоб-таки, да, так нет. Зато линии пересечения, однако…

Markiyan Hirnyk
Сообщения: 1366
Зарегистрирован: Вс дек 04, 2011 11:07 pm

Re: Прямо-таки

Сообщение Markiyan Hirnyk » Сб май 18, 2013 1:52 pm

алексей_алексей писал(а): Чтоб-таки, да, так нет. Зато линии пересечения, однако…
Процитирую Б.Н. Ельцина:"Ну вы даете..."