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

[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]

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

Используется синтаксис PostScript

%!PS-Adobe-2.0
%%BoundingBox: 0 0 595 760

% ----- This part of code is borrowed from PS-library, presented by "Алексей К./dxdy.ru" -----
% ----- Errors, iterations and main program (at end) are due to AKM/dxdy.ru ----- 

/Red  {1  0 0 setrgbcolor} bind def
/Red2 {1  0 1 setrgbcolor} bind def
/Green{0  .6 0 setrgbcolor} bind def
%/Green{0.64 0 0.95 0.40 setcmykcolor} bind def
/Blue {0  0 1 setrgbcolor} bind def
/Brown{1 .8 0 setrgbcolor} bind def
/Black{0  0 0 setrgbcolor} bind def

/ED {exch def} bind def

/XYadd {% x1 y1 x2 y2  --> x1+x2 y1+y2
  3 -1 roll add 3 1 roll add exch
} bind def

/XYsub {% x1 y1 x2 y2  --> x1-x2 y1-y2
  3 -1 roll sub neg 3 1 roll sub exch
} bind def

/Polar {%   r phi --> x=r*cos  y=r*sin
  2 copy cos mul 3 1 roll sin mul
} bind def 

/ToPolar {% x y --> r  phi
   exch 2 copy dup mul exch dup mul add sqrt 3 1 roll
   2 index 0. eq {pop pop 0.} {atan} ifelse
} bind def


/Rpoint {%  Abs.Radius   (curr.point) 
    gsave setlinewidth currentpoint newpath
          [] 0 setdash 1 setlinecap moveto 0 0 rlineto stroke
    grestore
} bind def

/forfor {%  x1 x2 N  forfor :  define float limits for "for" loop
    0.01 add 1 index 3 index sub exch div exch
} bind def

/Pstack {%
  count [exch dup (items on stack:) 3 -1 roll % S0 S1 ... S_n-1 [N (txt) N
      dup 3 add exch                          % S0 S1 ... S_n-1 [ i=n+3 n
      {dup index dup [ eq {pop (mark)} if exch} repeat
      pop] == flush
} bind def

/Pdict {%  dict
  count 0 eq {currentdict} {dup type /dicttype ne {currentdict} if}ifelse
  {[3 1 roll]==} forall
} bind def

/Args {%  /Title: N 
    [ 3 1 roll        % a1 ... aN [ /T N
      dup 2 add exch  % a1 ... aN [ /T N+2 N
      {dup index dup [ eq {pop (mark)} if exch} repeat pop % 2009!
    ] == flush
} bind def
% 0 1 2 3   /Test 3 Args

/XYdraw {% draw array [x y x y ...]  or --mark-- x y x y ...
   dup type /arraytype eq 
   {aload length}
   {counttomark dup 2 add -1 roll pop} ifelse  
   2 idiv 1 sub 3 1 roll moveto {lineto} repeat
} bind def


% -------------- Iteration for one time step (AKM):

/Tstep  {% [x1 y1 x2 y2  ... xN yN]
/I 0 def
/DT DT0 def
/DRmax 0 def
   [ exch aload length dup /NN ED
     -2 roll 4 copy NN 4 add 4 roll               % PN P1 P2 ... PN P1

     NN 2 idiv {%                                 % x1 y1 x2 y2 x3 y3
%/P1P2P3 6 Args 
        6 copy XYsub neg exch neg atan            % x1 y1 x2 y2 x3 y3 x1 y1 mu2
        3 1 roll                                  % x1 y1 x2 y2 x3 y3 mu2 x1 y1
        6 index 6 index XYsub neg exch neg atan   % x1 y1 x2 y2 x3 y3 mu2 mu1
%/mu2mu1 2 Args
        2 copy sub sin 5 1 roll add 4 1 roll      % x1 y1 x2 y2 mu1+mu2 sin(rho) x3 y3 
        7 index 7 index XYsub ToPolar             % x1 y1 x2 y2 mu1+mu2 sin(rho) h  angref 
%        gsave Blue newpath 1 setlinewidth 6 copy 4 2 roll pop pop 4 2 roll moveto  Polar rl^ stroke grestore  
%/ha 2 Args
        3 1 roll div 2 mul  3 1 roll 1 mul sub   % x1 y1 x2 y2  k tau
%        gsave Green newpath 1 setlinewidth 4 copy 4 2 roll moveto exch pop 20 exch Polar rl^ stroke grestore  
        exch DT mul exch 90 add Polar                  % x1 y1 x2 y2 dx2 dy2
       2 copy dup mul exch dup mul exch add sqrt
       dup DRmax lt {pop}{/DRmax ED} ifelse
%       /DT DT DRmax div 3 mul def
       gsave newpath Black .2 setlinewidth 
            4 copy 4 2 roll moveto rlineto stroke 
       grestore

        4 copy XYadd counttomark 2 roll           % [x2' y2' ... x1 y1 x2 y2 dx2 dy2
%        3 index 3 index moveto rlineto
        pop pop
/I I 1 add def
     } repeat
     pop pop pop pop] 

%[/DRmax DRmax] ==  dup length 2 div ==
} def


300 200 translate

%  ----------------------------------------------------   Ellipse picture

/a  250 def            % big   half-axis
/b  120 def            % small half-axis  
/DT0 6 def             % time step for one iteration
/Npt 36 def            % number of point on ellipse
/N 25 def              % numer of curves to draw
/M 100 def             % number of iterations between two drawings

/Xtest a def
/dt 360 Npt div def
/dt2 dt 2. div def
/Data0 [
   0 dt 360 dt2 sub {%
    dup cos a mul exch sin b mul
  } for
] def


.6 setlinewidth
Data0 XYdraw closepath stroke
Red


Data0 N {%
   M {Tstep} repeat
   dup  
   dup dup 0 get exch 1 get [ 3 1 roll (Dx=) Xtest 3 index sub] ==
   dup 0 get /Xtest ED
   XYdraw closepath stroke
   currentrgbcolor 3 1 roll setrgbcolor
} repeat pop
Black




%  ----------------------------------------------------------- Lemniscate picture

0 300 translate

/a  250 def
/b  250 def
/DT0 10 def
/Npt 30 def
/N 13 def
/M 50 def

/Npt Npt 2 div round 2 mul cvi def
/dt 180 Npt 2 sub div def
/dt2 dt 2. div def

/Data0 [
    0 0
   -45 dt2 add dt 45 {%
     dup 2 mul cos   % fi r^2
        sqrt exch Polar b mul exch a mul exch
  } for
    0 0
   -45 dt2 add dt 45 {%
     dup 2 mul cos   % fi r^2
        sqrt exch Polar b mul exch a mul neg exch
  } for
] def


0.6 setlinewidth
Data0 XYdraw closepath stroke
%Data0 aload length 2 idiv {moveto 3 Rpoint} repeat

Red

Data0 N
 {%
   M {Tstep} repeat
   dup  XYdraw closepath stroke
   currentrgbcolor 3 1 roll setrgbcolor
} repeat pop

flush showpage
% (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]

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