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

[Padawan]

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

Короче, смотрю я - полная каша

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

[paha]

Padawan в сообщении #304462 писал(а):

В итоге получается система дифференциальных уравнений с частными производными, типа уравнения теплопроводности

постановки корректной задачи для которой мы с нетерпением ожидаем

[ewert]

ДУ: $y'_t=\dfrac{y''_{xx}}{1+{y'_x}^2}$.

Граничные и начальные условия были приведены выше.

Я по-прежнему считаю, что степень суммы в знаменателе должна быть именно первая. Вполне возможно, что я и не прав (лень думать). Но даже если и так -- качественно свойства этой задачи от той степени не зависят.

[paha]

ewert в сообщении #304512 писал(а):

ДУ: $y'_t=\dfrac{y''_{xx}}{1+{y'_x}^2}$.

Граничные и начальные условия были приведены выше.

Уравнение увидел, оно имеет вид $F(y'_t,y'_x,y''_{xx})=0$. Осталось два вопроса

paha в сообщении #304266 писал(а):

1) Правда ли, что при каждом фиксированном $t$ проэволюционировавшая к этому моменту кривая $y(x,t)$ является графиком функции?

2) что такое <пересчитанная на вертикаль кривизна>?

Точно правильное уравнение написал Padawan вот тут (ниже $t$ -- параметр вдоль исходной кривой, а $\tau$ -- "время сдвига"):

Padawan в сообщении #301327 писал(а):

$$\frac{\partial\mathbf{r}}{\partial\tau}=\frac{\mathbf{r}''}{|\mathbf{r}'|^2}-\mathbf{r}'\frac{(\mathbf{r}',\mathbf{r}'')}{|\mathbf{r}'|^4}$$

Или, если проще, то $r'_{\tau}=kn$, где кривизна вычисляется по параметру $t$ и нормаль к этой же кривой.
Переписываю уравнение в координатах:

$$ x'_\tau=-k\frac{y'_t}{\sqrt{x'_t^2+y'_t^2}}=-\frac{(x'_ty''_{tt}-y'_tx''_{tt})y'_t}{(x'_t^2+y'_t^2)^2} $$

$$ y'_\tau=k\frac{x'_t}{\sqrt{x'_t^2+y'_t^2}}=\frac{(x'_ty''_{tt}-y'_tx''_{tt})x'_t}{(x'_t^2+y'_t^2)^2} $$

Знаки и степени тут правильные.

[Yu_K]

Маленькая простая моделька. Деформация замкнутой ломаной - для каждых трех соседних точек ломаной $(A_k,A_{k+1},A_{k+2})$ строится центр описанной окружности треугольника с вершинами $(A_k,A_{k+1},A_{k+2})$ , и вычисляется радиус описанной окружности - R и затем средняя точка $A_{k+1}$, смещается в сторону центра на расстояние пропорциональное 1/R, с одним и тем же коэффициентом пропорциональности. Грубая модель конечно. Мультик такой вот получился.

http://www.youtube.com/watch?v=5nV-CmrNLdg

Тест с окружностью проходит корректно - окружность схлопывается. Но а если взять, например, замкнутую кривую в форме "восьмерки" - то там возникают проблемы с бесконечной кривизной. Задача получается корректна только для достаточно хороших начальных кривых.

То paha

А как связать $\tau$ и $t$? Я тоже такую систему ДУ написал.

[paha]

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

А как связать $\tau$ и $t$? Я тоже такую систему ДУ написал.

Их вязать не надо, это независимые переменные, а не преступники:)

Мы же должны однопараметрическое семейство кривых $r(t,\tau)$ получить.

[AKM]

Лемниската и эллипс.

1. %!PS-Adobe-2.0 
2. %%BoundingBox: 0 0 595 760 
3.  
4. % ----- This part of code is borrowed from PS-library, presented by "Алексей К./dxdy.ru" ----- 
5. % ----- Errors, iterations and main program (at end) are due to AKM/dxdy.ru -----  
6.  
7. /Red {1 0 0 setrgbcolor} bind def 
8. /Red2 {1 0 1 setrgbcolor} bind def 
9. /Green{0 .6 0 setrgbcolor} bind def 
10. %/Green{0.64 0 0.95 0.40 setcmykcolor} bind def 
11. /Blue {0 0 1 setrgbcolor} bind def 
12. /Brown{1 .8 0 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.  
35. /Rpoint {% Abs.Radius (curr.point)  
36. gsave setlinewidth currentpoint newpath 
37. [] 0 setdash 1 setlinecap moveto 0 0 rlineto stroke 
38. grestore 
39. } bind def 
40.  
41. /forfor {% x1 x2 N forfor : define float limits for "for" loop 
42. 0.01 add 1 index 3 index sub exch div exch 
43. } bind def 
44.  
45. /Pstack {% 
46. count [exch dup (items on stack:) 3 -1 roll % S0 S1 ... S_n-1 [N (txt) N 
47. dup 3 add exch % S0 S1 ... S_n-1 [ i=n+3 n 
48. {dup index dup [ eq {pop (mark)} if exch} repeat 
49. pop] == flush 
50. } bind def 
51.  
52. /Pdict {% dict 
53. count 0 eq {currentdict} {dup type /dicttype ne {currentdict} if}ifelse 
54. {[3 1 roll]==} forall 
55. } bind def 
56.  
57. /Args {% /Title: N  
58. [ 3 1 roll % a1 ... aN [ /T N 
59. dup 2 add exch % a1 ... aN [ /T N+2 N 
60. {dup index dup [ eq {pop (mark)} if exch} repeat pop % 2009! 
61. ] == flush 
62. } bind def 
63. % 0 1 2 3 /Test 3 Args 
64.  
65. /XYdraw {% draw array [x y x y ...] or --mark-- x y x y ... 
66. dup type /arraytype eq  
67. {aload length} 
68. {counttomark dup 2 add -1 roll pop} ifelse  
69. 2 idiv 1 sub 3 1 roll moveto {lineto} repeat 
70. } bind def 
71.  
72.  
73. % -------------- Iteration for one time step (AKM): 
74.  
75. /Tstep {% [x1 y1 x2 y2 ... xN yN] 
76. /I 0 def 
77. /DT DT0 def 
78. /DRmax 0 def 
79. [ exch aload length dup /NN ED 
80. -2 roll 4 copy NN 4 add 4 roll % PN P1 P2 ... PN P1 
81.  
82. NN 2 idiv {% % x1 y1 x2 y2 x3 y3 
83. %/P1P2P3 6 Args  
84. 6 copy XYsub neg exch neg atan % x1 y1 x2 y2 x3 y3 x1 y1 mu2 
85. 3 1 roll % x1 y1 x2 y2 x3 y3 mu2 x1 y1 
86. 6 index 6 index XYsub neg exch neg atan % x1 y1 x2 y2 x3 y3 mu2 mu1 
87. %/mu2mu1 2 Args 
88. 2 copy sub sin 5 1 roll add 4 1 roll % x1 y1 x2 y2 mu1+mu2 sin(rho) x3 y3  
89. 7 index 7 index XYsub ToPolar % x1 y1 x2 y2 mu1+mu2 sin(rho) h angref  
90. % gsave Blue newpath 1 setlinewidth 6 copy 4 2 roll pop pop 4 2 roll moveto Polar rl^ stroke grestore  
91. %/ha 2 Args 
92. 3 1 roll div 2 mul 3 1 roll 1 mul sub % x1 y1 x2 y2 k tau 
93. % gsave Green newpath 1 setlinewidth 4 copy 4 2 roll moveto exch pop 20 exch Polar rl^ stroke grestore  
94. exch DT mul exch 90 add Polar % x1 y1 x2 y2 dx2 dy2 
95. 2 copy dup mul exch dup mul exch add sqrt 
96. dup DRmax lt {pop}{/DRmax ED} ifelse 
97. % /DT DT DRmax div 3 mul def 
98. gsave newpath Black .2 setlinewidth  
99. 4 copy 4 2 roll moveto rlineto stroke  
100. grestore 
101.  
102. 4 copy XYadd counttomark 2 roll % [x2' y2' ... x1 y1 x2 y2 dx2 dy2 
103. % 3 index 3 index moveto rlineto 
104. pop pop 
105. /I I 1 add def 
106. } repeat 
107. pop pop pop pop]  
108.  
109. %[/DRmax DRmax] == dup length 2 div == 
110. } def 
111.  
112.  
113. 300 200 translate 
114.  
115. % ---------------------------------------------------- Ellipse picture 
116.  
117. /a 250 def % big half-axis 
118. /b 120 def % small half-axis  
119. /DT0 6 def % time step for one iteration 
120. /Npt 36 def % number of point on ellipse 
121. /N 25 def % numer of curves to draw 
122. /M 100 def % number of iterations between two drawings 
123.  
124. /Xtest a def 
125. /dt 360 Npt div def 
126. /dt2 dt 2. div def 
127. /Data0 [ 
128. 0 dt 360 dt2 sub {% 
129. dup cos a mul exch sin b mul 
130. } for 
131. ] def 
132.  
133.  
134. .6 setlinewidth 
135. Data0 XYdraw closepath stroke 
136. Red 
137.  
138.  
139. Data0 N {% 
140. M {Tstep} repeat 
141. dup  
142. dup dup 0 get exch 1 get [ 3 1 roll (Dx=) Xtest 3 index sub] == 
143. dup 0 get /Xtest ED 
144. XYdraw closepath stroke 
145. currentrgbcolor 3 1 roll setrgbcolor 
146. } repeat pop 
147. Black 
148.  
149.  
150.  
151.  
152. % ----------------------------------------------------------- Lemniscate picture 
153.  
154. 0 300 translate 
155.  
156. /a 250 def 
157. /b 250 def 
158. /DT0 10 def 
159. /Npt 30 def 
160. /N 13 def 
161. /M 50 def 
162.  
163. /Npt Npt 2 div round 2 mul cvi def 
164. /dt 180 Npt 2 sub div def 
165. /dt2 dt 2. div def 
166.  
167. /Data0 [ 
168. 0 0 
169. -45 dt2 add dt 45 {% 
170. dup 2 mul cos % fi r^2 
171. sqrt exch Polar b mul exch a mul exch 
172. } for 
173. 0 0 
174. -45 dt2 add dt 45 {% 
175. dup 2 mul cos % fi r^2 
176. sqrt exch Polar b mul exch a mul neg exch 
177. } for 
178. ] def 
179.  
180.  
181. 0.6 setlinewidth 
182. Data0 XYdraw closepath stroke 
183. %Data0 aload length 2 idiv {moveto 3 Rpoint} repeat 
184.  
185. Red 
186.  
187. Data0 N 
188. {% 
189. M {Tstep} repeat 
190. dup XYdraw closepath stroke 
191. currentrgbcolor 3 1 roll setrgbcolor 
192. } repeat pop 
193.  
194. flush showpage 
195. % (eof) 

[ewert]

paha в сообщении #304564 писал(а):

Точно правильное уравнение написал Padawan вот тут

В правильности сомневаться трудно, однако для численного моделирования выгоднее переписать это уравнение так:

$${\partial\vec r\over\partial\tau}={1\over|\vec{r_t}'|}\cdot\left({\vec{r_t}'\over|\vec{r_t}'|}\right)'_t.$$

(Оно, кстати, и геометрически очевиднее.) А для краевой задачи разумнее перейти от параметрического задания к явной зависимости от иксов: $y=y(x,\tau)$. Тогда $\displaystyle{\partial\vec r\over\partial\tau}=\vec n\cdot{y''_{xx}\over\left(1+{y'_x}^2\right)^{3/2}}$. Здесь $|\vec n|=1$, $|y'_{\tau}|=\dfrac{|\vec{r_{\tau}}'|}{\cos\alpha}$ и $\dfrac{1}{\cos\alpha}=\sqrt{1+\tg^2\alpha}=\sqrt{1+{y'_x}^2}$. Откуда

$$y'_{\tau}={y''_{xx}\over1+{y'_x}^2}.$$

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

Задача получается корректна только для достаточно хороших начальных кривых.

По-видимому, для любых кривых без самопересечений. А с самопересечениями, конечно, будут особенности.

[Yu_K]

ewert в сообщении #304783 писал(а):

По-видимому, для любых кривых без самопересечений. А с самопересечениями, конечно, будут особенности.

У лемнискаты точка ее самопересечения остается на месте и там нет проблем.
Проблемы возникают в процессе деформации кривой с быстрыми переходами от малых радиусов к большим - точки в окрестности малых радиусов (там где большая кривизна) начинают слишком быстро "вырываться вперед" и некоторая часть кривой начинает отрываться от части движущейся более медленно, и там начинают появляться "клювики".

[ewert]

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

точки в окрестности малых радиусов (там где большая кривизна) начинают слишком быстро "вырываться вперед"

Куда вырываться?... кривизны же выравниваются. Т.е. если в какой-то точке радиус был мал -- именно за счёт её ускоренного движения по сравнению с соседними радиус в этой точке начнёт увеличиваться.

Собственно, радиус в начале может быть и просто нулевым. Задайте петлю в форме капельки (т.е. с уголком). Этот уголок мгновенно закруглится и пойдёт дальше сглаживаться, сглаживаться.

[paha]

ewert в сообщении #304783 писал(а):

$|y'_{\tau}|=\dfrac{|\vec{r_{\tau}}'|}{\cos\alpha}$

Я наверное туплю спросонок, но разве катет больше гипотенузы?

[ewert]

paha в сообщении #304832 писал(а):

, но разве катет больше гипотенузы?

Нет, но зато гипотенуза больше катета. А перемещение от точки на кривой к соседней кривой по вертикали -- это именно гипотенуза (относительно перемещения к той же кривой по нормали).

[paha]

Итак, уравнение

ewert в сообщении #304783 писал(а):

$$y'_{t}={y''_{xx}\over1+{y'_x}^2}.$$

и

ewert в сообщении #304145 писал(а):

задача:

$\displaystyle\begin{cases}{\partial y\over\partial t}=\text{<пересчитанная на вертикаль кривизна>}; \\ y\big|_{x=0}=0,\quad y\big|_{x=1}=0; \\ y\big|_{t=0}=w(x).\end{cases}$

<пересчитанная на вертикаль кривизна> - это и есть правая часть нашего уравнения, т.е. Вы поставили такую задачу?

$\displaystyle\begin{cases}y'_{t}={y''_{xx}\over1+{y'_x}^2}; \\ y\big|_{x=0}=0,\quad y\big|_{x=1}=0; \\ y\big|_{t=0}=w(x).\end{cases}$

[ewert]

Да.

[vvvv]

Кривизна кривой величина скалярная, неясно , что значит
пересчитать на вертикаль и почему именно на вертикаль?