сфера nd

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

Модератор: Admin

uni
Сообщения: 1817
Зарегистрирован: Сб ноя 13, 2004 3:06 pm
Откуда: п.г.т. Излучинск
Контактная информация:

Сообщение uni » Вс окт 02, 2011 11:01 pm

Если всё дело в формуле гиперсферы, то её точки, на мой взгляд, можно вычислить при помощи трёх вложенных циклов не ущемляя общности задачи для произвольного n при фиксированном количестве шагов по параметрам.

Оценим количество свободных параметров в формуле: это n - размерность сферы (количество координат), a - это коэффициенты при параметрическом описании сферы и последний - N - это количество шагов по каждому параметру. Т.е. мы перебираем все коэффициенты по количеству шагов и вычисляем все координаты точек по порядку. Получится прямоугольная таблица.

Ничего особо сложного не вижу, FOR_DO можно отложить для специальных случаев. Вам в коде ещё чуть не хватило опыта для допиливания его до нужной кондиции.

Должно быть три цикла while do. В самом верхнем цикле мы перебираем шаги по параметрам a(i), в следующем цикле мы перебираем сами коэффициенты a, а в последнем перебираем уже координаты точек u. Итак, три цикла, три индекса по параметрам (шаги, коэффициенты, точки).

Получится, что в теле самого вложенного цикла мы пробегаем все координаты для фиксированного коэффициента и фиксированного шага для этого коэффициента. На втором уровне вложенности мы пробегаем все коэффициенты для всех точек при фиксированном шаге. А на первом уровне мы пробегаем все шаги для всех коэффициентов для всех точек.

Циклы, как мне думается, можно менять местами. Попробую реализовать эту идею.

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

Сообщение алексей_алексей » Вс окт 02, 2011 11:33 pm

Слава, знаю, только, что точек должно быть N^(n-1) (что само по себе может оказаться проблематичным для памяти, а потому, наверное, надо разбивать поверхность на части.) Но в остальном голова уже не варит, и как при трёх циклах получится столько точек пока тоже не соображу… Хорошо, конечно, получись оно без премудростей…

uni
Сообщения: 1817
Зарегистрирован: Сб ноя 13, 2004 3:06 pm
Откуда: п.г.т. Излучинск
Контактная информация:

Сообщение uni » Пн окт 03, 2011 2:39 am

Вроде бы это должно выглядеть так.

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

restart:

# Подключаемые библиотеки
with( plots ):

# Параметры:
R := 1.0:   # - радиус сферы

N := 3:      # - размерность сферы
M := 50:    # - количество точек вдоль одной координаты (сферической)
K := N - 1: # - количество параметров при параметрическом задании формулы

# Индексы
m := 1:     # - индекс по шагам
k := 1:     # - индекс по коэффициентам a[i]
n := 1:     # - индекс по текущим координатам точек

# Шаг 1. Зададим начальные значения для коэффициентов a[i]:
for k from 1 to K do a[k] := 1: end do:

# Шаг 2. Вычисляем точки гиперсферы

# Инициализируем счётчики
m := 1:     # - индекс по шагам
k := 1:     # - индекс по коэффициентам a[i]
n := 1:     # - индекс по текущим координатам точек

# Постоянный множитель
L := 2 * Pi / M:

# Цикл по шагам и коэффициентам a[i] одновременно
while m <= M ^ K do

    # Цикл по координатам точек
    while n <= N do

        # Первое уравнение не содержит функции cos(), поэтому обрабатываем особо
        if ( n = 1 ) then
           
            u[ m, 1 ] := R * ( product( sin( L * a[jj] ), jj = 1 .. K ) ):
       
        else
       
            u[ m, n ] := R * cos( L * a[ n - 1 ] ) * ( product( sin( L * a[jj] ) , jj = n .. K ) ):
       
        end if:
       
        n := n + 1:
   
    end do:
   
    # Готовим счётчик к следующей итерации
    n := 1:
   
    # Увеличиваем шаг
    m := m + 1:
   
    # Главный трюк
    # Здесь коэффициенты a[i] имеют целые значения, которые обновляются как в зацепном
    # счётчике. Можно представить коэффициенты a[i] в виде вектора целых значений, где
    # самый младший "разряд" перебирает значения от 1 до M, а остальные увеличиваются
    # в зависимости от переполнения предыдущего коэффициента.
    for k from 1 to K do
   
        if ( k = 1 ) then
           
            a[1] := m mod M + 1:
       
        else
           
            a[k] := floor( m / ( k * M ) ) + 1:
           
        end if:
   
    end do:

end do:

# Подготавливаем расчётную таблицу к отображению
L := seq( [ u[ii, 1], u[ii, 2], u[ii, 3] ], ii = 1 .. m - 1 ):
 
# Рисуем картинку
pointplot3d( [L], axes = normal, color = red );

Дайте мне медаль :)

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

Сообщение алексей_алексей » Пн окт 03, 2011 3:53 pm

Орден, Слава, только орден или звезду Героя. Сообщение видел сутра, но немного позапускать смог только что. Не все операторы мне ещё знакомы – буду смотреть. А работает в разы быстрее моего текста – наверное, из-за процедуры? Вечером постараюсь заняться вплотную…

Kitonum
Сообщения: 2084
Зарегистрирован: Ср дек 31, 2008 1:55 pm
Откуда: г. Пенза

Сообщение Kitonum » Пн окт 03, 2011 7:56 pm

алексей_алексей писал(а):...только орден или звезду Героя...

Позволю себе немного подпортить атмосферу всеобщего ликования! В коде уважаемого uni заметил ряд шероховатостей:

1) Размерность сферы определяется числом независимых параметров, её определяющих. Для обычной сферы в трёхмерном пространстве N=2 а не 3.

2) В списке точек, лежащих на сфере, отсутствует "северный полюс", т.е. точка (0,0,R), а другие точки повторяются, особенно часто "южный полюс" (0,0,-R). В самом деле
nops([L]);
100

nops({L});
42

Конечно, эти недочёты не принципиальны и их легко исправить. Более интересно, что сам код уважаемого uni мохно заметно подсократить. Если использовать рекурсию, то можно вообще обойтись без циклических конструкций. Легко заметить, что N-мерная сфера получается из (N-1)-мерной сферы умножением каждой координаты на синус соответствующего параметра и добавлением координаты, равной произведению радиуса на косинус того же параметра. На этом сообрахении основана следующая процедура с именем P, которая вычисляет точки на N-мерной сфере для любого N>=1. Смысл параметра M несколько другой, чем в коде uni. На M частей делится промежуток от 0 до Pi. Так удобнее, чтобы исключить повторы.

Код процедуры:

P:=proc(R, N, M)
P(R,1,M):=[seq([R*cos(Pi*k/M),R*sin(Pi*k/M)],k=0..2*M-1)];
if N>1 then P(R,N,M):=[[seq(0,j=1..N),1],seq(seq([seq(P (R,N-1,M)[k][n]*sin(Pi*m/M),n=1..N),R*cos(Pi*m/M)],k=1..nops(P(R,N-1,M))),m=1..M-1),[seq(0,j=1..N),-1]];
fi;
end proc:


Примеры использования:

Список точек:

P(1,2,10);

Рисунок:

plots[pointplot3d](P(1,2,10),symbol=box,symbolsize=4, color = red, axes = normal, view = [-1.9 .. 1.9, -1.9 .. 1.9, -1.4 .. 1.4],scaling=constrained);

Анимация рисунка:

Seq:=seq(plots[pointplot3d](P(1,2,10)[1..k],symbol=box,symbolsize=4, color = red, axes = normal, view = [-1.9 .. 1.9, -1.9 .. 1.9, -1.4 .. 1.4]),k=1..nops(P(1,2,10))):
plots[display](Seq,insequence=true,scaling=constrained);

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

Сообщение алексей_алексей » Пн окт 03, 2011 9:01 pm

Kitonum писал(а):...

У Славы отдельной строкой выделено число параметров, а именно: K:=N-1. То есть, размерность пространства минус 1. Но дело не в этом. Всё обнаружится в процессе отладки и доводки, если у участвующих хватит терпения. Больше всего важна скорость счёта, ну, понятно, как там растёт время в зависимости от размерности. Вряд ли удастся непосредственно на сфере находить точку пересечения с линиями невязок – точки надо будет уточнять отдельным способом (в голове держу пакет Моисеева)… Сравните, если Вам не трудно, время счёта своей процедуры и, несомненно (он же был первым), орденоносца Вячеслава Николаевича…
Для сокращения времени счёта думается на этапе вычисления точек сферы ввести ограничение до 2-3-х значащих цифр. Чичас буду это пробовать…

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

Сообщение алексей_алексей » Пн окт 03, 2011 9:12 pm

О, Digits очень сильно влияет на время счёта…

Kitonum
Сообщения: 2084
Зарегистрирован: Ср дек 31, 2008 1:55 pm
Откуда: г. Пенза

Сообщение Kitonum » Пн окт 03, 2011 9:43 pm

алексей_алексей писал(а):… Сравните, если Вам не трудно, время счёта своей процедуры и, несомненно (он же был первым), орденоносца Вячеслава Николаевича…
Для сокращения времени счёта думается на этапе вычисления точек сферы ввести ограничение до 2-3-х значащих цифр...

Время счёта окончательного списка у uni для R=1, N=4, M=40:
st:= time(): evalf[3]([seq( [seq( u[ii, k], k=1..N)], ii = 1 .. m - 1 )]): time() - st;
2.870

Время счёта того же примера у меня:
st:= time(): evalf[3](P(1,3,20)): time() - st;
0.827

Софт: Maple 13 (классический интерфейс)
Компьютер: Athlon, 2 Ггц, Ram 2 Ггб.

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

Сообщение алексей_алексей » Пн окт 03, 2011 10:00 pm

Kitonum писал(а):Время счёта окончательного списка у uni для R=1, N=4, M=40:
st:= time(): evalf[3]([seq( [seq( u[ii, k], k=1..N)], ii = 1 .. m - 1 )]): time() - st;
2.870

Время счёта того же примера у меня:
st:= time(): evalf[3](P(1,3,20)): time() - st;
0.827

Софт: Maple 13 (классический интерфейс)
Компьютер: Athlon, 2 Ггц, Ram 2 Ггб.


Хотел поработать с Вашей процедурой, но выдаёт ошибку Warning, `P` is implicitly declared local to procedure `P` Да, а в классике работает… Непонятно.
А так разница во времени впечатляет…

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

Сообщение алексей_алексей » Пн окт 03, 2011 10:20 pm

Kitonum писал(а):Время счёта того же примера у меня...

Да у Вас просто мгновенно считает… Как бы текст в привычную среду поместить?

uni
Сообщения: 1817
Зарегистрирован: Сб ноя 13, 2004 3:06 pm
Откуда: п.г.т. Излучинск
Контактная информация:

Сообщение uni » Пн окт 03, 2011 10:43 pm

Да, эта штука работает и быстро. Вот чем мне нравится Maple, когда можно сократить обычный код до такого уровня.

Буду потом в каждом листинге указывать свою версию пакета (у меня Maple 10).

Синусы, косинусы и повторные вычисления всегда напрашиваются на быструю реализацию. Это второй этап после написания работающего кода. Спасибо за замечания, я их учту на будущее.

Половину (лучшую) ордена отдаю уважаемому Kitonum, т.к. лучше я вряд ли напишу. Спецом по Maple не являюсь.

Kitonum
Сообщения: 2084
Зарегистрирован: Ср дек 31, 2008 1:55 pm
Откуда: г. Пенза

Сообщение Kitonum » Пн окт 03, 2011 11:43 pm

алексей_алексей писал(а):...выдаёт ошибку…

Мне тоже непонятно, почему один и тот же код, скопированный в классику, работает, а в стандарт - нет. Почему-то в стандарте некорректно работает механизм вызова процедуры внутри самой себя (рекурсия). Что-то создатели "нахимичили", либо я что-то не понимаю!

Следующий вариант кода работает в стандарте:

P:= proc (R, N, M)
local A, t;
A[1] := [seq([R*cos(Pi*k/M), R*sin(Pi*k/M)], k = 0 .. 2*M-1)];
if 1 < N then for t from 2 to N do A[t] := [[seq(0, j = 1 .. t), 1], seq(seq([seq(A[t-1][k][n]*sin(Pi*m/M), n = 1 .. t), R*cos(Pi*m/M)], k = 1 .. nops(A[t-1])), m = 1 .. M-1), [seq(0, j = 1 .. t), -1]] end do end if;
A[N];
end proc:


Время счёта примерно то же.

Kitonum
Сообщения: 2084
Зарегистрирован: Ср дек 31, 2008 1:55 pm
Откуда: г. Пенза

Сообщение Kitonum » Вт окт 04, 2011 12:24 am

Kitonum писал(а):... Что-то создатели "нахимичили"...

Зря я погрешил на разработчиков! Предыдущий вариант кода (с рекурсией) тоже работает в стандарте, если убрать пробел между P и следующей скобкой ( в выражении P (R,N-1,M) . Непонятно как он там оказался!

Подправленный вариант кода (при инициализации выбрать remember table assignment) :

P:=proc(R, N, M)
P(R,1,M):=[seq([R*cos(Pi*k/M),R*sin(Pi*k/M)],k=0..2*M-1)];
if N>1 then
P(R,N,M):=[[seq(0,j=1..N),1],seq(seq([seq(P(R,N-1,M)[k][n]*sin(Pi*m/M),n=1..N),R*cos(Pi*m/M)],k=1..nops(P(R,N-1,M))),m=1..M-1),[seq(0,j=1..N),-1]]; fi;
end proc:

uni
Сообщения: 1817
Зарегистрирован: Сб ноя 13, 2004 3:06 pm
Откуда: п.г.т. Излучинск
Контактная информация:

Сообщение uni » Вт окт 04, 2011 12:28 am

Это, кстати, была бы неплохая задачка для изучающих Maple и программирование в нём. Три разных подхода к решению одной задачи (без FOR_DO).

Ещё замечу. При цифровой обработке сигналов функции sin() (и cos(), что то же самое) не вычисляются. Там обычно используют заранее посчитанную с нужной точностью таблицу синуса, а косинус получается сдвигом на Pi/2, т.е. индекс в таблице берётся со смещением в четверть.

Я не знаю как вычисляет эти функции Maple, но, думаю, что он делает это "честно". Можно попробовать ускорить алгоритм не вычисляя функции, а обращаясь к заранее подготовленной таблице нужного размера. Такая вот идея.

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

Сообщение алексей_алексей » Вт окт 04, 2011 9:01 pm

Никогда не пользовался рекуррентными (или как правильно, рекурсией?) кодами, и немного посмотрел, как оно пишется. (Не понял, зачем мы перед оператором конца процедуры пишем A[N]?) Наверное, исходя из универсальности, лучше использовать текст для стандартного интерфейса, хотя первый вариант понятней. И по ощущениям время счёта в классике меньше, не смотря на одинаковые показатели счётчика.
uni писал(а): Ещё замечу. При цифровой обработке сигналов функции sin() (и cos(), что то же самое) не вычисляются. Там обычно используют заранее посчитанную с нужной точностью таблицу синуса, а косинус получается сдвигом на Pi/2, т.е. индекс в таблице берётся со смещением в четверть.

Я не знаю как вычисляет эти функции Maple, но, думаю, что он делает это "честно". Можно попробовать ускорить алгоритм не вычисляя функции, а обращаясь к заранее подготовленной таблице нужного размера. Такая вот идея.

А как это реализовать?

Господа орденоносцы, будем продолжать или остановимся на абстрактной сфере? Путь – полное численное решение полиномиальных систем со СНУ на горизонте…