Моделирование преобразования кривой...

[Padawan]

alexey007
Если для поверхностей, то лучше брать среднюю кривизну - полусумму главных кривизн. Так как она имеет смысл величины давления, вызываемого поверхностным натяжением, если мне не изменяет память. То есть Вы опишите динамику мыльной пленки.

[alexey007]

vvvv в сообщении #306534 писал(а):

А что, с плоской кривой уже полная ясность?

Да, все ясно, уравнения получены для плоской кривой, нужен численный метод чтобы решить и все, что еще тут не ясно может быть?

[vvvv]

Вот, вот - решить уравнения и посмотреть на кривую

[alexey007]

vvvv в сообщении #306556 писал(а):

Вот, вот - решить уравнения и посмотреть на кривую

Вот щас занимаюсь разбором численной схемы которая хорошо бы решала эту задачу. vvvv, эволюцию какой кривой вы бы хотели увидеть?

[AKM]

Повторюсь --- желающие взглянуть на кривульки, и владеющие Маплом, могут себя удовлетворить, положив что-то вроде Digits:=40 и тыщу точек на кривой. Не знаю, есть ли такая фича в Матлабе, например. Я лично подожду, пока кто-нибудь нарисует или код выложит --- мне лень изучать там рисовалки, и код переписывать.
Если не выложат, и кончатся детективы --- сосчитать точки в Мапле и экспортнуть в свой PS сумею.

Грубых ошибок у меня, похоже, нет.

AKM в сообщении #303590 писал(а):

Замкнутая окружность у меня схлопывается/расширяется по закону

$$R(t)=\sqrt{R_0^2\pm 2t}.$$

Радиус 256 должен уменьшиться вдвое через $t=3\cdot 32\cdot 256$. Берём a=b=256, DT0=3, M=256, N=32, заранее пририсовываем итоговую окружность радиуса 128, и к ней схлопываемся:

%  ----------------------------------------------------   Ellipse picture
/a  256 def            % big   half-axis
/b  256 def            % small half-axis  
/DT0 3 def             % time step for one iteration
/Npt 50 def            % number of point on ellipse
/N 32 def              % numer of curves to draw
/M 256 def             % number of iterations between two drawings
% So, total number of iterations is M*N

 0 0 128 0 360 arc stroke %   <-----------------------  ОКРУЖНОСТЬ радиуса 128

Этот самопроверочный тест, возможно, сгодится и автору темы.

[Yu_K]

АКМ Эволюцию, например, такой полярной кривой $r(t)=2+cos(12t)$ - хотелось бы увидеть. Просто картинку без кода покажите.

В дискретном аналоге исходной задачи нет шага $t$ по параметру кривой - а я то думал, что тут для устойчивости счета понадобится условие похожее на условие Куранта для гиперболических уравнений $u*\tau/h<1$.

[alexey007]

Yu_K, пишите пожалуйста крививые в парметрическом виде в полярных координатах $p\in[0,2\pi]$ $x(p,0)=f(p), y(p,0)=g(p)$, у меня в коде так сделано, а переделывать, на это времени нет

[AKM]

Yu_K в сообщении #306572 писал(а):

АКМ Эволюцию, например, такой полярной кривой $r(t)=2+cos(12t)$ - хотелось бы увидеть. Просто картинку без кода покажите.

Картинку без кода не смогу --- это дольше (код-то всё равно надо написать) и, оказывается, неинтересно. Интересно вьючить в GSview --- как бы анимация получается.

Замысел Ваш понятен. Мне тоже казалось, что всё должно схлопываться в окружность, и здесь это здорово получается. Только я взял не 12, а 6.

1. %!PS-Adobe-2.0 
2. %%BoundingBox: 0 0 595 760 
3.  
4. % ==================== !!! ОЙ КАК СМЕШНО ПОЛУЧАЕТСЯ !!! ===================== 
5.  
6. % ----- This part of code is borrowed from PS-library, presented by "Алексей К./dxdy.ru" ----- 
7. % ----- Errors, iterations and main program (at end) are due to AKM/dxdy.ru -----  
8.  
9. /Red {1 0 0 setrgbcolor} bind def 
10. /Red2 {1 0 1 setrgbcolor} bind def 
11. /Green{0 .6 0 setrgbcolor} bind def 
12. %/Green{0.64 0 0.95 0.40 setcmykcolor} bind def 
13. /Blue {0 0 1 setrgbcolor} bind def 
14. /Brown{1 .8 0 setrgbcolor} bind def 
15. /Black{0 0 0 setrgbcolor} bind def 
16.  
17. /ED {exch def} bind def 
18.  
19. /XYadd {% x1 y1 x2 y2 --> x1+x2 y1+y2 
20. 3 -1 roll add 3 1 roll add exch 
21. } bind def 
22.  
23. /XYsub {% x1 y1 x2 y2 --> x1-x2 y1-y2 
24. 3 -1 roll sub neg 3 1 roll sub exch 
25. } bind def 
26.  
27. /Polar {% r phi --> x=r*cos y=r*sin 
28. 2 copy cos mul 3 1 roll sin mul 
29. } bind def  
30.  
31. /ToPolar {% x y --> r phi 
32. exch 2 copy dup mul exch dup mul add sqrt 3 1 roll 
33. 2 index 0. eq {pop pop 0.} {atan} ifelse 
34. } bind def 
35.  
36.  
37. /Rpoint {% Abs.Radius (curr.point)  
38. gsave setlinewidth currentpoint newpath 
39. [] 0 setdash 1 setlinecap moveto 0 0 rlineto stroke 
40. grestore 
41. } bind def 
42.  
43. /forfor {% x1 x2 N forfor : define float limits for "for" loop 
44. 0.01 add 1 index 3 index sub exch div exch 
45. } bind def 
46.  
47. /Pstack {% 
48. count [exch dup (items on stack:) 3 -1 roll % S0 S1 ... S_n-1 [N (txt) N 
49. dup 3 add exch % S0 S1 ... S_n-1 [ i=n+3 n 
50. {dup index dup [ eq {pop (mark)} if exch} repeat 
51. pop] == flush 
52. } bind def 
53.  
54. /Pdict {% dict 
55. count 0 eq {currentdict} {dup type /dicttype ne {currentdict} if}ifelse 
56. {[3 1 roll]==} forall 
57. } bind def 
58.  
59. /Args {% /Title: N  
60. [ 3 1 roll % a1 ... aN [ /T N 
61. dup 2 add exch % a1 ... aN [ /T N+2 N 
62. {dup index dup [ eq {pop (mark)} if exch} repeat pop % 2009! 
63. ] == flush 
64. } bind def 
65. % 0 1 2 3 /Test 3 Args 
66.  
67. /XYdraw {% draw array [x y x y ...] or --mark-- x y x y ... 
68. dup type /arraytype eq  
69. {aload length} 
70. {counttomark dup 2 add -1 roll pop} ifelse  
71. 2 idiv 1 sub 3 1 roll moveto {lineto} repeat 
72. } bind def 
73.  
74.  
75. % -------------- Iteration for one time step (AKM): 
76.  
77. /Tstep {% [x1 y1 x2 y2 ... xN yN] 
78. /I 0 def 
79. /DT DT0 def 
80. /DRmax 0 def 
81. [ exch aload length dup /NN ED 
82. -2 roll 4 copy NN 4 add 4 roll % PN P1 P2 ... PN P1 
83.  
84. NN 2 idiv {% % x1 y1 x2 y2 x3 y3 
85. %/P1P2P3 6 Args  
86. 6 copy XYsub neg exch neg atan % x1 y1 x2 y2 x3 y3 x1 y1 mu2 
87. 3 1 roll % x1 y1 x2 y2 x3 y3 mu2 x1 y1 
88. 6 index 6 index XYsub neg exch neg atan % x1 y1 x2 y2 x3 y3 mu2 mu1 
89. %/mu2mu1 2 Args 
90. 2 copy sub sin 5 1 roll add 4 1 roll % x1 y1 x2 y2 mu1+mu2 sin(rho) x3 y3  
91. 7 index 7 index XYsub ToPolar % x1 y1 x2 y2 mu1+mu2 sin(rho) h angref  
92. % gsave Blue newpath 1 setlinewidth 6 copy 4 2 roll pop pop 4 2 roll moveto Polar rl^ stroke grestore  
93. %/ha 2 Args 
94. 3 1 roll div 2 mul 3 1 roll 1 mul sub % x1 y1 x2 y2 k tau 
95. % gsave Green newpath 1 setlinewidth 4 copy 4 2 roll moveto exch pop 20 exch Polar rl^ stroke grestore  
96. exch DT mul exch 90 add Polar % x1 y1 x2 y2 dx2 dy2 
97. 2 copy dup mul exch dup mul exch add sqrt 
98. dup DRmax lt {pop}{/DRmax ED} ifelse 
99. % /DT DT DRmax div 3 mul def 
100. gsave newpath Black .2 setlinewidth  
101. 4 copy 4 2 roll moveto rlineto stroke  
102. grestore 
103.  
104. 4 copy XYadd counttomark 2 roll % [x2' y2' ... x1 y1 x2 y2 dx2 dy2 
105. % 3 index 3 index moveto rlineto 
106. pop pop 
107. /I I 1 add def 
108. } repeat 
109. pop pop pop pop]  
110.  
111. %[/DRmax DRmax] == dup length 2 div == 
112. } def 
113.  
114.  
115. 305 320 translate 
116.  
117. % ---------------------------------------------------- Yu_K picture 
118.  
119. % R(t) = a + b*cos(wt)  
120.  
121. /a 200 def  
122. /b 100 def  
123. /w 6 def 
124. /DT0 3 def % time step for one iteration 
125. /Npt 120 def % number of points on ellipse 
126. /N 20 def % numer of curves to draw 
127. /M 400 def % number of iterations between two drawings 
128. % So, total number of iterations is M*N 
129.  
130. % Нет, точки перегиба лежат не на этой окружности; её надо аккуратно сосчитать... 
131. gsave .2 setlinewidth [3 3] 0 setdash 0 0 a 0 360 arc stroke grestore 
132.  
133.  
134. /Xtest a def 
135. /dt 360 Npt div def 
136. /dt2 dt 2. div def 
137. /Data0 [ 
138. 0 dt 360 dt2 sub {% 
139. dup w mul cos b mul a add exch Polar 
140. } for 
141. ] def 
142.  
143. 1.2 setlinewidth 
144. Data0 XYdraw closepath stroke 
145.  
146. Red .6 setlinewidth 
147.  
148. Data0 N {% 
149. M {Tstep} repeat 
150. dup  
151. dup dup 0 get exch 1 get [ 3 1 roll (Dx=) Xtest 3 index sub] == 
152. dup 0 get /Xtest ED 
153. XYdraw closepath stroke 
154. currentrgbcolor 3 1 roll setrgbcolor 
155. } repeat pop 
156. flush showpage 
157. % (eof) 
158.  
159.  

-- Пн апр 05, 2010 17:16:53 --

Красное цветовыделение --- свойство тэга syntax, а вовсе не злоупотребление модераторскими правами.

(Кому ещё какую кривульку забацать?)

Я стараюсь взять себя в руки, и на 2-5 дней удалиться от темы. С форума вряд ли выгонят, а с работы точно!
alexey007, Вас я в этом винить не буду, это я сам такой безответственный. Как бы в игрушки весёленькие играю вместо того, чтобы...

[Padawan]

AKM
Красота!

[vvvv]

АКМ, но мочему сперва кривая рисуется во впадинах наружу,а в выступах вовнутрь. Должно, по-моему, или вовнутрь или наружу.
А, ясно - к центру кривизны

[AKM]

Padawan, я рад, что Вам понравилось. Сам оторваться не могу.
Но всё равно не уверен, что при более аккуратном моделировании результат будет таким же.

vvvv в сообщении #305703 писал(а):

Файл АКМ у меня не открывается -(неправильный формат или повреждена) ?

Да, действительно, в Опере появляется коррумпированный zip и не открывается, а IE нормально открывает eps.
Но его можно вытащить и чисто кнопкой "ВЫДЕЛИТЬ".

-- Пн апр 05, 2010 19:03:56 --

Надо бы координаты точек перегиба сосчитать (это же так легко --- запузырить соотв. темку "Срочно! Помогите найти точку перегиба!") и там точек наставить. Пунктирная окружность по замыслу была ГМТ точек перегиба, но на самом деле это где $\cos t=0$, а перегибы наверняка оттуда переехали...

[vvvv]

АКМ, у меня оба файла открылись и я картинки посмотрел.Действительно, смотрится красиво

[Yu_K]

У меня здесь нет возможности посмотреть. vvvv, можете показать картинки.

И вот еще одна кривая, которая мне интересна $x(t,0)=sin(t)-2sin(2t), y(t,0)=cos(t)+2cos(2t)$ отсюда
http://www.mathcurve.com/courbes2d/hypotrochoid/hypotrochoid.shtml.
В Маткаде не получалось ничего хорошего сделать при грубой явной аппроксимации. Можно конечно пробовать неявную схему. АКМ а Вы какой алгоритм используете при построении?

[ewert]

(Оффтоп)

AKM в сообщении #306535 писал(а):

Или когда детективы в магазине кончились.

Ещё было бы забавно убедить ewert'a, что оно в квадрате, или наоборот.

детективы там давно уж закончились; давно уж все напечатаны и распроданы -- во всяком случае, стОящие (по смыслу, а не по ценникам).

И не убедите. Я неверующий.

А насчёт мусора -- скорее всего, так. По мере сжатия предельной окружности, наверное, ужесточаются требования устойчивости к временному шагу. Поскольку радиус окружности схлопывается с ускорением (там типа $r(t)\sim\sqrt{C-t}$). Поэтому с ростом времени следует уменьшать временной шаг во избежание неустойчивости. Впрочем, серьёзно я это не продумывал.

[AKM]

Yu_K в сообщении #306628 писал(а):

АКМ а Вы какой алгоритм используете при построении?

Тупо провожу ориентированную окружность через соседние 3 точки, использую её кривизну и нормаль в средней точке. Нормаль будет кажется

$$ \arg(z_2-z_1)+\arg(z_3-z_2)-\arg(z_3-z_1)+\frac\pi2 $$

(т.е. касательная + 90). Там же в коде всё вроде написано!

-- Пн апр 05, 2010 22:00:22 --

Берём

$$x(t)=r(t)\cos t,\quad y(t)=r(t)\sin t, \quad\text{где}\quad r(t)=a+b\cos(4 t),\quad b>a.$$

Т.е. ф-ция $r(t)$ отрицательной бывает.

1. %!PS-Adobe-2.0 
2. %%BoundingBox: 0 0 595 760 
3.  
4. % Вот сволочь, --- чего вытворяет! Через ноль на скорости проскочила...  -------- 
5.  
6. % ----- This part of code is borrowed from PS-library, presented by "Алексей К./dxdy.ru" ----- 
7. % ----- Errors, iterations and main program (at end) are due to AKM/dxdy.ru -----  
8.  
9. /Red  {1  0 0 setrgbcolor} bind def 
10. /Red2 {1  0 1 setrgbcolor} bind def 
11. /Green{0  .6 0 setrgbcolor} bind def 
12. /Blue {0  0 1 setrgbcolor} bind def 
13. /Black{0  0 0 setrgbcolor} bind def 
14.  
15. /ED {exch def} bind def 
16.  
17. /XYadd {% x1 y1 x2 y2  --> x1+x2 y1+y2 
18.   3 -1 roll add 3 1 roll add exch 
19. } bind def 
20.  
21. /XYsub {% x1 y1 x2 y2  --> x1-x2 y1-y2 
22.   3 -1 roll sub neg 3 1 roll sub exch 
23. } bind def 
24.  
25. /Polar {%   r phi --> x=r*cos  y=r*sin 
26.   2 copy cos mul 3 1 roll sin mul 
27. } bind def  
28.  
29. /ToPolar {% x y --> r  phi 
30.    exch 2 copy dup mul exch dup mul add sqrt 3 1 roll 
31.    2 index 0. eq {pop pop 0.} {atan} ifelse 
32. } bind def 
33.  
34. /XYdraw {% draw array [x y x y ...]  or --mark-- x y x y ... 
35.    dup type /arraytype eq  
36.    {aload length} 
37.    {counttomark dup 2 add -1 roll pop} ifelse   
38.    2 idiv 1 sub 3 1 roll moveto {lineto} repeat 
39. } bind def 
40.  
41. % -------------- Iteration for one time step (AKM): 
42.  
43. /Tstep  {% [x1 y1 x2 y2  ... xN yN] 
44. /I 0 def 
45. /DT DT0 def 
46. /DRmax 0 def 
47.    [ exch aload length dup /NN ED 
48.      -2 roll 4 copy NN 4 add 4 roll               % PN P1 P2 ... PN P1 
49.  
50.      NN 2 idiv {%                                 % x1 y1 x2 y2 x3 y3 
51.         6 copy XYsub neg exch neg atan            % x1 y1 x2 y2 x3 y3 x1 y1 mu2 
52.         3 1 roll                                  % x1 y1 x2 y2 x3 y3 mu2 x1 y1 
53.         6 index 6 index XYsub neg exch neg atan   % x1 y1 x2 y2 x3 y3 mu2 mu1 
54.         2 copy sub sin 5 1 roll add 4 1 roll      % x1 y1 x2 y2 mu1+mu2 sin(rho) x3 y3  
55.         7 index 7 index XYsub ToPolar             % x1 y1 x2 y2 mu1+mu2 sin(rho) h  angref  
56.         3 1 roll div 2 mul  3 1 roll 1 mul sub   % x1 y1 x2 y2  k tau 
57.         exch DT mul exch 90 add Polar                  % x1 y1 x2 y2 dx2 dy2 
58.        2 copy dup mul exch dup mul exch add sqrt 
59.        dup DRmax lt {pop}{/DRmax ED} ifelse 
60.        gsave newpath Black .2 setlinewidth  
61.             4 copy 4 2 roll moveto rlineto stroke  
62.        grestore 
63.  
64.         4 copy XYadd counttomark 2 roll           % [x2' y2' ... x1 y1 x2 y2 dx2 dy2 
65.         pop pop 
66. /I I 1 add def 
67.      } repeat 
68.      pop pop pop pop]  
69. } def 
70.  
71.  
72. 320 320 translate  
73.  
74. %  ----------------------------------------------------   Yu_K picture 
75.  
76. %              R(t) = a + b*cos(wt)     
77.  
78. /a  100 def            %  
79. /b  320 def            %    
80. /w  4 def 
81. /DT0 2 def             % time step for one iteration 
82. /Npt 128 def           % number of points on ellipse 
83. /N 25 def              % numer of curves to draw 
84. /M 300 def             % number of iterations between two drawings 
85.                        % So, total number of iterations is M*N 
86.  
87. 45 rotate 
88.  
89.  
90. /Xtest a def 
91. /dt 360 Npt div def 
92. /dt2 dt 2. div def 
93. /Data0 [ 
94.    0 dt 360 dt2 sub {% 
95.      dup w mul cos b mul a add exch Polar 
96.   } for 
97. ] def 
98.  
99. 1 setlinewidth 
100. Data0 XYdraw closepath stroke 
101. Red 
102.  
103. .6 setlinewidth 
104. Data0 N {% 
105.    M {Tstep} repeat 
106.    dup   
107.    dup dup 0 get exch 1 get [ 3 1 roll (Dx=) Xtest 3 index sub] == 
108.    dup 0 get /Xtest ED 
109.    XYdraw closepath stroke 
110.    currentrgbcolor 3 1 roll setrgbcolor 
111. } repeat pop 
112. Black 
113.  
114. flush showpage 
115. % (eof)