uzluga.ru
добавить свой файл
1
Построение пространственных кривых, заданных неявно, с применением программы Mathematica


А.П.Мостовской


Пакет программ Mathematica, версий 4 и 5, компании Wolfram Research Inc. позволяет строить множества в пространстве или на плоскости, заданные разными способами. Так, на плоскости (2D - построения) с помощью встроенных или подключенных функций пакета можно построить множества, заданные явно (команда Plot) или неявно (ImplicitPlot), множества, заданные параметрическими уравнениями (ParametricPlot), в полярной системе координат и так далее. В пространстве (3D - построения) применяются аналоги таких команд - Plot3D, ParametricPlot3D, для построения множеств, заданных неявно одним уравнением, применяется команда ContourPlot3D. Есть и другие полезные программы.

Одним из пробелов пакета Mathematica в этом направлении является отсутствие встроенной или загружаемой команды, выполняющей построение кривой в 3-х мерном пространстве, заданной естественным образом, неявно - как решение системы двух уравнений. Например, нельзя построить в Mathematica линию Вивиани, заданную как пересечение сферы единичного радиуса и цилиндра, содержащего центр сферы и имеющего единичный диаметр:



В данной статье приводится программа в Mathematica 4(5), позволяющая строить кривые в пространстве, заданные неявно. При этом, кривую можно построить отдельно или вместе с каркасами порождающих ее поверхностей.

Идея построения кривой заключена в следующем: пространство разбивается параллельными плоскостями на одинаковые непересекающиеся кубы. Находятся многоугольники (полигоны) пересечения каждого куба и одной из поверхностей. Затем определяются отрезки-пересечения найденных многоугольников (составляющих первую поверхность) со второй поверхностью. Объединение найденных отрезков составляет искомую кривую.


1. Вспомогательные задачи


Рассмотрим несколько задач, необходимых в дальнейшем.

Длина вектора. Определим функцию norm[x], вычисляющую длину вектора, следующим образом:


norm[x_] := Sqrt[x.x]


Преобразование списка. Пусть ломаная задана своими вершинами gt. Составим список сторон данной ломаной в порядке их следования. Это можно сделать командой


neu[gt_] := Block[{vv = {}, ty = Length[gt] - 1},

Do[vv = Append[vv, {gt[[i]], gt[[i + 1]]}],

{i, 1, Length[gt] - 1}];

vv = Append[vv, {gt[[Length[gt]]], gt[[1]]}]]


Построение выпуклого многоугольника по его вершинам. Пусть даны координаты вершин плоского выпуклого многоугольника, например, e={{0.8, 0.5, 3.}, {-1, 0.1, 3.}, {0.7, -0.7, 3.}, {-0.1, 1, 3.}, {-0.5, -0.9, 3.}, {0.9, 0.4, 3.}}}.

Требуется построить многоугольник с вершинами в данных точках. Для построения многоугольника с вершинами в точках, координаты которых представлены в виде списка e, в программе Mathematica надо применить графический примитив Polygon в следующем виде:


Show[Graphics3D[Polygon[e]], Boxed -> False,

PlotRange -> {{-1, 1}, {-1, 1}, {2.8, 3.3}}]}


Результат действия команды приведен на рисунке. Следовательно, перед применением команды построения многоугольника, надо упорядочить вершины многоугольника в порядке их следования по границе многоугольника, если число вершин многоугольника Length[e] больше трех.

Для этого перейдем от точек к векторам. Вычтем из координат каждой точки списка e соответствующие координаты первой точки, то есть точки e[[1]].

Это можно сделать командой Map с применением "чистой" функции: (# - e[[1]]) & /@ e (замена каждого элемента списка q значением функции D[x] от элемента можно командой Map (другое ее обозначение /@) по следующему правилу: Map[D[#]&,q] (или D[#]&/@q)).

Получим список из векторов, в котором первый элемент будет нулевой вектор. Вырежем этот вектор из писка командой Drop. Получим список ww:


ww = Drop[(# - e[[1]])& /@ e, 1]


Рассмотрим теперь ориентированные углы между векторами списка ww. Для этого предположим, что векторов в списке по крайней мере два и они не коллинеарны. Возьмем два первых вектора списка ww[[1]] и ww[[2]] и построим нормальный вектор к плоскости многоугольника


nm = Cross[ww[[1]] - ww[[2]], ww[[1]] + ww[[2]]];


Под ориентированным углом между векторами, например, ww[[1]]

и ww[[4]] будем понимать следующее число


Sign[Det[{nm, ww[[1]], ww[[4]]}]]

ArcCos[ww[[1]].ww[[4]] /(norm[ww[[1]]]norm[ww[[4]]])]

Воспользуемся понятием "чистой" функции, командой Map и получим список u ориентированных углов между первым вектором ww[[1]] и всеми векторами списка ww:


u = (Sign[Det[{nm, ww[[1]],#}]]ArcCos[ww[[1]].# / (norm[ww[[1]]] norm[#])])& / @ ww


Полученный список ориентированных углов упорядочим по возрастанию:


u1 = Sort[u]


Пусть fg = Length[u1] - длина списка u1 и qq = {} - вспомогательный пустой список. Рассмотрим цикл


Do[po = Position[u, u1[[k]]][[1, 1]];

ve = Take[e, {po + 1}][[1]];

qq = Append[qq, ve], {k, 1, fg}];


Команда Position определяет позицию po вхождения элемента u1[[k]], k=1,2,...,fg, в список u, ve - вершина многоугольника из списка e соответствующая позиции po. Вершину ve прибавляем справа к списку qq командой Append. Наконец, прибавим к списку qq начальную точку e[[1]] и получим упорядоченный список состоящий из координат вершин многоугольника:


qq = Append[qq, e[[1]]]


Теперь команда


Show[Graphics3D[Polygon[qq]], Boxed -> False,

PlotRange -> {{-1, 1}, {-1, 1}, {2.8, 3.3}}]


нарисует выпуклый многоугольник.

Объединим все команды в единую функцию следующим образом:


er[e_]: = Block[{u1},

If[Length[e] > 3, ww = Drop[(# - e[[1]])&/@ e, 1];

nm = Cross[ww[[1]] - ww[[2]], ww[[1]] + ww[[2]]];

u = (Sign[Det[{nm, ww[[1]],#}]]ArcCos[ww[[1]].# / (norm[ww[[1]]] norm[#])])&/@ ww;

u1 = Sort[u];

fg = Length[u1];

qq = {};

Do[po = Position[u, u1[[k]]][[1, 1]];

ve = Take[e, {po + 1}][[1]];

qq = Append[qq, ve], {k, 1, fg}];

qq = Append[qq, e[[1]]];, qq = e]; qq]


2. Пересечение отрезка и поверхности


Пусть множество F задана уравнением вида f(x,y,z)=0,где f(x,y,z) - непрерывная функция, отрезок AB, , достаточно мал, точки A и B достаточно близки к поверхности. Если , то отрезок AB пересекает поверхность.

Действительно, координаты u=(x,y,z) радиус-вектора произвольной точки M отрезка можно выразить через координаты радиус-векторов концов отрезка:





где , а . Тогда функция

(1)

непрерывна и следовательно обращается в ноль при некотором ,

если на концах отрезка принимает значения разных знаков.

Предположим теперь, что поверхность F гладкая (функция f(x,y,z) дифференцируема на множестве F) и найдем отношение, в котором точка M пересечения поверхности и отрезка AB делит отрезок AB. Пусть есть проекция точки A на поверхность F. Применяя равенство (1) для точек A, и формулу Тейлора, можно записать, пренебрегая членами малости второго порядка

(2)

так как ; здесь - расстояние от точки A до поверхности F. Так как уравнение касательной плоскости к поверхности F в точке имеет вид



то расстояние от точки A до поверхности, равное расстоянию от точки A до касательной плоскости в точке , можно представить в виде




где



- нормальный вектор к поверхности F в точке . Отсюда и из (2) получаем, что f(A)==|n()|. Аналогично можно записать f(B)==|n()|.

Пренебрегая членами малости второго порядка, можно считать, что треугольники AM и BM подобны, а нормальные векторы в точках A и B равны. Отсюда получим, что



и следовательно координаты точки M можно представить в виде

(3)

В программе Mathematica координаты точки пересечения поверхности и отрезка можно вычислить по формуле (3) с помощью следующей функции:


xyzf[p_]:= Block[{a,b},a= Abs[f[p[[1]]]]; b = Abs[f[p[[2]]]];

If[a + b == 0, p, {p[[1]] b + p[[2]] a}/(a + b)]];


где p - список из координат концов отрезка, f[x,y,z] - функция, задающая поверхность. При этом будем считать, что пересечением будет весь отрезок, если концы отрезка принадлежат поверхности: a+b=|f(A)|+|f(B)|=0.


3. Построение кривой, заданной неявно


В этом параграфе составим сначала программу построения поверхности F, заданную неявно, то есть уравнением вида: f(x,y,z)=0. Эта программа позволит нам решить основную задачу – построение кривой, заданной неявно.

Рассмотрим куб с длиной ребра h и со сторонами параллельными осям координат. Если вершина куба {i, j, k} имеет минимальные координаты по отношению к соответствующим координатам остальных вершин куба, то координаты всех остальных вершин куба можно получить

прибавляя к соответствующим координатам длину ребра h. Поэтому следующий список


cb:={{{i,j,k}, {i+h,j,k}}, {{i+h,j,k}, {i+h,j+h,k}}, {{i+h,j+h,k}, {i,j+h,k}}, {{i,j + h,k}, {i,j,k}}, {{i,j,k+h}, {i+h,j,k+h}},{{i+ h,j,k+h}, {i+h,j+h,k+h}}, {{i+h,j+h,k+h}, {i,j+h,k+h}}, {{i, j+h, k+h}, {i,j,k+h}}, {{i,j,k}, {i,j,k+h}}, {{i+h,j,k}, {i+h, j, k+h}}, {{i+h,j+h,k}, {i+h,j+h,k+h}}, {{i,j+h,k}, {i,j+h,k+h}}};


порождает список, состоящий из координат концов ребер куба, если задана минимальная вершина куба {i, j, k} и длина его ребра h.

Для простоты будем рисовать только часть поверхности, заключенную в куб с центром в начале координат, и гранями, параллельными координатным плоскостям. Пусть длина ребра такого куба равна 2m. Разобьем такой куб плоскостями, параллельными координатным плоскостям, на меньшие непересекающиеся кубы с длиной ребра h, и составим список, элементами которого будут списки сторон (заданных своими концевыми точками) каждого куба разбиения. Это можно сделать следующей командой:


p = Flatten[Table[cb, {i, -m, m - h, h}, {j, -m, m - h, h}, {k, -m, m - h, h}], 2]

Выберем командой Select из списка p только те элементы, которые связаны с кубами, пересекающими поверхность F. Для этого надо рассмотреть элементы списка p. Каждый такой элемент содержит список концевых точек ребер соответствующего куба.

Список q, образованный следующей командой, содержит стороны только тех кубов, которые пересекают поверхность.


q = Select[p, MemberQ[Map[f[#1[[1]]]f[#1[[2]]] <=

0&,#], True] &}]

В команде MemberQ[Map[f[#1[[1]]]f[#1[[2]]]<=0&,#],True]& значениями переменной # являются по очереди все элементы списка p.

Команда Map[f[#1[[1]]]f[#1[[2]]] <= 0& воздействует на каждую сторону элемента списка p - проверяет, верно ли утверждение f[#1[[1]]]f[#1[[2]]]<=0, где "чистая" переменные #1 будут заменены всякий раз координатами концевых точек одной стороны. Если неравенство истинно, то соответствующий список из сторон присоединяется к списку q.

Из каждого такого списка команда Select выбирает стороны, удовлетворяющие условию f[#1[[1]]]f[#1[[2]]] <= 0,


evrr = Map[Select[#, f[#1[[1]]]f[#1[[2]]] ≤ 0 &] &, q];


где вместо решетки #1 подставляется список из концов каждой стороны. Здесь, #1[[1]] - это координаты конца одной стороны, а #1[[2]] - координаты второго конца этой же стороны. Список evrr будет состоять из тех ребер списка q, которые пересекают поверхность F.

Каждый отрезок из списка evrr заменяем точкой пересечения отрезка и поверхности, используя функцию xyzf, построенную в п. 2. Появившиеся в результате этой процедуры лишние скобки убираем командами Flatten, после чего разбиваем полученный список на тройки координат точек пересечения командой Partition. В результпте получим новый список rre:


rre = Map[Partition[Flatten[Map[xyzf, #], 2], 3] &, evrr];


Уберем из каждого элемента списка rre дублирующие координаты:


wer1 = Map[Union[#] &, rre];


Список wer1 состоит из полигонов, образованных пересечением маленьких кубиков и поверхности F. Теперь можно построить каркас поверхности F. Для этого упорядочим точки каждого элемента списка wer1 командой er из п. 2:


wer = Map[er[N[#]]&,wer1]


Добавим к каждому элементу этого списка его начальную точку и, тем самым, замкнем каждый полигон:


wer2 = Map[Append[#, First[#]] &, wer]


Теперь можно построить каркас поверхности F:


d1 = Graphics3D[Line /@ wer2]


Вернемся к списку wer. Каждый элемент этого списка есть список вершин одного из полигонов. Составим из каждого такого списка вершин список сторон полигона:


rre2=Map[neu[#]&,wer];


Из списка rre2 выберем те полигоны, которые пересекаю вторую поверхность G, заданную, например, уравнением g(x,y,z)=0:


q=Select[rre2,MemberQ[Map[g[#1[[1]]]g[#1[[2]]]0&,#],True]&];


Из каждого элемента списка q возьмем только те стороны, который пересекают поверхность G:


evrr2 = Map[Select[#, g[#1[[1]]]g[#1[[2]]] ≤ 0 &] &, q];


Каждый отрезок списка evrr2 заменяем его точкой пересечения с поверхностью G:


rre3 = Map[Partition[Flatten[Map[xyzg, #], 2], 3] &, evrr2];


Теперь можно построить кривую пересечения поверхностей F и G:


d3 = Graphics3D[{Hue[0.7], Thickness[0.01], Line /@ rre3}];

Show[d3]


Сделаем дополнение к написанной программе. Приведем программу позволяющую нарисовать кривую и каркасы поверхностей ее образующих.

Построим каркас поверхности G по аналогии с построением каркаса поверхности F:


q=Select[p,MemberQ[Map[g[#1[[1]]]g[#1[[2]]]≤0&,#],True] &];

evrr = Map[Select[#, g[#1[[1]]]g[#1[[2]]] ≤ 0 &] &, q];

rre = Map[Partition[Flatten[Map[xyzg, #], 2], 3] &, evrr];

werg = Map[er[#] &, rre];

werg = Map[Append[#, First[#]] &, werg];

d2 = Graphics3D[Line /@ werg];


Здесь функция xyzg аналогична функции xyzf.

Нарисуем теперь кривую и каркасы поверхностей:


Show[d1, d2, d3]


4. Пример. Построение линии Вивиани.


Приведем полный список команд построения кривой.

Сначала надо выполнить следующие вспомогательные команды:


xyzf[p_]:=

Block[{a,b},a=Abs[N[f[p[[1]]]]];b=Abs[N[f[p[[2]]]]];

If[a + b == 0, {p[[1]], p[[2]]},{p[[1]] b + p[[2]] a}/

(a + b)]];

xyzg[p_]:=

Block[{a, b},a = Abs[N[g[p[[1]]]]];b=Abs[N[g[p[[2]]]]];

If[a + b == 0, {p[[1]], p[[2]]}, {p[[1]] b + p[[2]] a}/

(a + b)]];

norm[x_] := Sqrt[x.x] // N

er[e_] :=

Block[{u1}, If[Length[e] > 3, ww = Drop[(# - e[[1]]) & /@ e, 1];

nm = Cross[ww[[1]], ww[[2]]];

u = (Sign[Det[{nm, ww[[1]], #}]]ArcCos[

ww[[1]].#/(norm[ww[[1]] ]norm[#])]) & /@ ww;

u1 = Sort[u];

fg = Length[u1];

qq = {};

Do[po = Position[u, u1[[k]]][[1, 1]];

ve = Take[e, {po + 1}][[1]];

qq = Append[qq, ve], {k, 1, fg}];

qq = Append[qq, e[[1]]];, qq = e]; qq]

cb :={{{i, j, k}, {i + h, j, k}}, {{i + h, j, k}, {i +

h, j + h, k}}, {{i + h, j + h, k}, {i, j + h, k}}, {{i, j+h,

k}, {i, j, k}}, {{i,j, k + h}, {i + h, j, k + h}}, {{i + h, j, k + h}, {i + h, j + h, k + h}}, {{i + h, j + h, k + h}, {i, j + h, k + h}}, {{i, j + h, k + h}, {i, j, k + h}}, {{i, j, k}, {i, j, k + h}}, {{i + h, j, k}, {i + h, j, k + h}}, {{i + h, j + h, k}, {i + h, j + h, k + h}}, {{i, j + h, k}, {i, j + h, k + h}}};

neu[gt_] :=

Block[{vv = {}, ty = Length[gt] -

1}, Do[vv = Append[vv, {gt[[i]],

gt[[i + 1]]}], {i, 1, Length[gt] - 1}];

vv = Append[vv, {gt[[Length[gt]]], gt[[1]]}]]

f[{x1_, y1_, z1_}] := ff /. {x -> x1, y -> y1, z -> z1};

g[{x1_, y1_, z1_}] := gg /. {x -> x1, y -> y1, z -> z1};


Теперь можно ввести функции:


ff =(x-1/2)^2+y^2-1/4;

gg=x^2+y^2+z^2-1;


С точки зрения программы Mathematica это будут выражения. Две последние строки в предыдущем блоке превращают эти выражения в функции.

Далее вводим команды построения кривой и кривой и поверхностей:

m = 1;

h = .2001;

p = Flatten[Table[cb, {i, -m, m, h}, {j, -m, m, h}, {k, -m, m, h}], 2];

q = Select[p, MemberQ[Map[f[#1[[1]]]f[#1[[2]]] ≤

0 &, {{#[[1, 1]], #[[6, 2]]}, {#[[7, 2]], #[[1, 2]]}, {#[[9, 2]], #[[2,

2]]}, {#[[3, 2]], #[[6, 1]]}}], True] &];

evrr = Map[Select[#, f[#1[[1]]]f[#1[[2]]] ≤ 0 &] &, q];

rre = Map[Partition[Flatten[Map[xyzf, #], 2], 3] &, evrr];

wer1 = Map[Union[#] &, rre];

wer = Map[er[#] &, wer1];

rre2 = Map[neu[#] &, wer];

q = Select[rre2, MemberQ[Map[g[#1[[1]]]g[#1[[2]]] ≤ 0 &, #], True] &];

evrr2 = Map[Select[#, g[#1[[1]]]g[#1[[2]]] ≤ 0 &] &, q];

rre3 = Map[Partition[Flatten[Map[xyzg, #], 2], 3] &, evrr2];

d3 = Graphics3D[{Hue[0.7], Thickness[0.01], Line /@ rre3}];

Show[d3]

q = Select[p, MemberQ[Map[g[#1[[1]]]g[#1[[2]]] ≤ 0 &, #], True] &];

evrr = Map[Select[#, g[#1[[1]]]g[#1[[2]]] ≤ 0 &] &, q];

rre = Map[Partition[Flatten[Map[xyzg, #], 2], 3] &, evrr];

werg = Map[er[#] &, rre];

werg = Map[Append[#, First[#]] &, werg];

d2 = Graphics3D[Line /@ werg];

Show[d1, d2, d3]


В результате получим:




Здесь выбрано m=1, шаг (и размеры маленьких кубов) h=.2001.