Описание бильярда

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

Модератор: Admin

Markiyan Hirnyk
Сообщения: 1366
Зарегистрирован: Вс дек 04, 2011 11:07 pm

Описание бильярда

Сообщение Markiyan Hirnyk » Пн фев 24, 2014 12:24 pm

Возможно ли описание двухмерного бильярда в единичном квадрате в терминах events и if (см. http://ru.wikipedia.org/wiki/%D0%91%D0% ... 1%80%D0%B4)? Например, задана динамическая система
{x''(t)=0,y''(t)=0, x(0)=0.5, x'(0)=0.2, y(0)=0.5, y'(0)= sqrt(3)/5},
при достижении границы {x=0}, {x=1}, {y=0}, {y=1} траектория {x(t),y(t)} отражается по закону Гюйгенса, т. е. угол падения равен углу отражения, без потери энергии.
Последний раз редактировалось Markiyan Hirnyk Пн фев 24, 2014 5:23 pm, всего редактировалось 1 раз.

Markiyan Hirnyk
Сообщения: 1366
Зарегистрирован: Вс дек 04, 2011 11:07 pm

Ответ

Сообщение Markiyan Hirnyk » Пн фев 24, 2014 5:12 pm

Ответ на MaplePrimes http://www.mapleprimes.com/questions/20 ... eply=reply
Тот же вопрос для единичного круга остается открытым.

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

Re: Ответ

Сообщение алексей_алексей » Пн фев 24, 2014 6:05 pm

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

restart; with(plots):
 sys := {diff(x(t), t, t) = 0, diff(y(t), t, t) = 0, x(0) = 0., y(0) = -0., (D(x))(0) = .2, (D(y))(0) = (1/5)*sqrt(3)}; res := dsolve(sys, numeric, events = [[x(t)^2+y(t)^2 = 1, [diff(x(t), t) = -.5*(diff(x(t), t))-.5*(diff(y(t), t)), diff(y(t), t) = -.5*(diff(x(t), t))-.5*(diff(y(t), t))]]]); a := odeplot(res, [x(t), y(t)], 0 .. 100, axes = boxed); b := implicitplot(x^2+y^2 = 1, x = -1 .. 1, y = -1 .. 1); display(a, b)

Markiyan Hirnyk
Сообщения: 1366
Зарегистрирован: Вс дек 04, 2011 11:07 pm

Закон Гюйгенса

Сообщение Markiyan Hirnyk » Пн фев 24, 2014 6:10 pm

Алексей Борисович!
Спасбо за внимание к моему вопросу и Ваш труд. Код надо подправить : на рисунке угол отражения не равен углу падения.
Изображение

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

Re: Ответ

Сообщение алексей_алексей » Пн фев 24, 2014 6:14 pm

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

restart: with(plots):
 sys := {diff(x(t), t, t) = 0, diff(y(t), t, t) = 0, x(0) = 0., y(0) = -0., (D(x))(0) = .5, (D(y))(0) = -.5}; res := dsolve(sys, numeric, events = [[x(t)^2+y(t)^2 = 1, [diff(x(t), t) = -.5*(diff(x(t), t))-.5*(diff(y(t), t)), diff(y(t), t) = -.5*(diff(x(t), t))-.5*(diff(y(t), t))]]]); a := odeplot(res, [x(t), y(t)], -100 .. 100, axes = boxed); b := implicitplot(x^2+y^2 = 1, x = -1 .. 1, y = -1 .. 1); display(a, b)

Markiyan Hirnyk
Сообщения: 1366
Зарегистрирован: Вс дек 04, 2011 11:07 pm

Почти то

Сообщение Markiyan Hirnyk » Пн фев 24, 2014 6:18 pm

Качество рисунка оставляет желать лучшего:
Изображение
Алексей Борисович!
Пожалуйста, ответьте и на MaplePrimes.

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

Re: Почти то

Сообщение алексей_алексей » Пн фев 24, 2014 6:28 pm

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

restart;
with(plots):
sys := {diff(x(t), t, t) = 0, diff(y(t), t, t) = 0, x(0) = 0., y(0) = -0., (D(x))(0) = 1, (D(y))(0) = -1}; res := dsolve(sys, numeric, events = [[x(t)^2+y(t)^2 = 1, [diff(x(t), t) = -.5*(diff(x(t), t))-.5*(diff(y(t), t)), diff(y(t), t) = -.5*(diff(x(t), t))-.5*(diff(y(t), t))]]]);
a := odeplot(res, [x(t), y(t)], -150 .. 150, axes = boxed, numpoints = 5000); b := implicitplot(x^2+y^2 = 1, x = -1 .. 1, y = -1 .. 1);
display(a, b);


Markiyan Hirnyk
Сообщения: 1366
Зарегистрирован: Вс дек 04, 2011 11:07 pm

То же самое

Сообщение Markiyan Hirnyk » Пн фев 24, 2014 6:30 pm

Угол падения не равен углу отражения.

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

Сообщение Kitonum » Пн фев 24, 2014 8:11 pm

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

restart;
A:=plots[animate](plot,[[2/Pi*abs(arcsin(sin(Pi/2*(0.5+0.2*t)))), 2/Pi*abs(arcsin(sin(Pi/2*(0.5+sqrt(3)/5*t)))), t=0..a], thickness=5, color=red, numpoints=1000], a=0..20, frames=200):
B:=plottools[rectangle]([0,1], [1,0], thickness=2, color=green):
plots[display](A, B, scaling=constrained, axes=none);


Результат:
http://s019.radikal.ru/i640/1402/20/74830c66a048.gif

Markiyan Hirnyk
Сообщения: 1366
Зарегистрирован: Вс дек 04, 2011 11:07 pm

Re

Сообщение Markiyan Hirnyk » Пн фев 24, 2014 9:26 pm

Kitonum писал(а):Для упругого шара, если нет потери энергии, задача по-сути носит чисто геометрический характер. Вот решение, не использующее дифуры. Оно легко обобщается на любой прямоугольный стол:

restart;
A:=plots[animate](plot,[[2/Pi*abs(arcsin(sin(Pi/2*(0.5+0.2*t)))), 2/Pi*abs(arcsin(sin(Pi/2*(0.5+sqrt(3)/5*t)))), t=0..a], thickness=5, color=red, numpoints=1000], a=0..20, frames=200):
B:=plottools[rectangle]([0,1], [1,0], thickness=2, color=green):
plots[display](A, B, scaling=constrained, axes=none);


Результат:
http://s019.radikal.ru/i640/1402/20/74830c66a048.gif
Спасибо. Однако цель вопроса - освоить примененение events.

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

Re: Re

Сообщение Kitonum » Пн фев 24, 2014 10:49 pm

Markiyan Hirnyk писал(а):...Спасибо. Однако цель вопроса - освоить применение events.

Осваивайте, Вам никто не мешает. Однако знание явного уравнения траектории, я думаю, также полезно.

Markiyan Hirnyk
Сообщения: 1366
Зарегистрирован: Вс дек 04, 2011 11:07 pm

Сообщение Markiyan Hirnyk » Ср фев 26, 2014 8:43 am

алексей_алексей писал(а):

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

restart: with(plots):
 sys := {diff(x(t), t, t) = 0, diff(y(t), t, t) = 0, x(0) = 0., y(0) = -0., (D(x))(0) = .5, (D(y))(0) = -.5}; res := dsolve(sys, numeric, events = [[x(t)^2+y(t)^2 = 1, [diff(x(t), t) = -.5*(diff(x(t), t))-.5*(diff(y(t), t)), diff(y(t), t) = -.5*(diff(x(t), t))-.5*(diff(y(t), t))]]]); a := odeplot(res, [x(t), y(t)], -100 .. 100, axes = boxed); b := implicitplot(x^2+y^2 = 1, x = -1 .. 1, y = -1 .. 1); display(a, b)
Сделано на MaplеPrimes двумя способами: переходя к полярным координатам и в декартовых координатах. Вторая реализация сложнее, т. к. нелинейно изменяются обе производные и приходится вводить дополнительную переменную. Ее преимущество состоит в том, что таким образом можно задавать бильярды в произвльных выпуклых областях, ограниченных гладкой кривой и в произвольных выпуклых многоугольниках. Алексей Борисович, может Вы ради спортивного (и не только спортивного) интереса попробуете создать бильярд в выпуклом четырехугольнике, не являющемся прямоугольником.

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

не могу

Сообщение алексей_алексей » Ср фев 26, 2014 11:18 am

Вставить условием достижение границы для кривой и для выпуклого многоугольника, думаю, можно. В последнем случае это видится рутинным перечислением всех прямых. Не знаю, как задать изменение направления скорости относительно нормали к кривой в точке, поменяв угол alpha на “-”alpha пересечения с ней (или к касательной на Pi+2*(Pi/2-alpha)). Точно так же не знаю, как поступить в аналогичном случае с прямыми, задаваемыми сторонами выпуклого многоугольника.

Markiyan Hirnyk
Сообщения: 1366
Зарегистрирован: Вс дек 04, 2011 11:07 pm

Re: не могу

Сообщение Markiyan Hirnyk » Ср фев 26, 2014 11:22 am

алексей_алексей писал(а): Не знаю, как задать изменение направления скорости относительно нормали к кривой в точке, поменяв угол alpha на “-”alpha пересечения с ней (или к касательной на Pi+2*(Pi/2-alpha)). Точно так же не знаю, как поступить в аналогичном случае с прямыми, задаваемыми сторонами выпуклого многоугольника.
Это указано на MaplePrimes в комментарии Preben Alsholm. Алексей Борисович, Вы заглядывали туда, прежде чем скоропалительно отвечать?

Markiyan Hirnyk
Сообщения: 1366
Зарегистрирован: Вс дек 04, 2011 11:07 pm

Моя неудачная попытка

Сообщение Markiyan Hirnyk » Чт фев 27, 2014 8:05 pm

описать бильярд в трапеции О(0,0), А(0,1), В(1,2), С(1,0):

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

restart; sys := {diff(x(t), t, t) = 0, diff(y(t), t, t) = 0, x(0) = .5, y(0) = .5, (D(x))(0) = .2, (D(y))(0) = .1*sqrt(3)};
res := dsolve(sys, numeric, events = [[x(t) = 0, diff(x(t), t) = -(diff(x(t), t))], [x(t) = 1, diff(x(t), t) = -(diff(x(t), t))], [y(t) = 0, diff(y(t), t) = -(diff(y(t), t))], [y(t)-1-x(t) = 0, [temp = diff(x(t), t), diff(x(t), t) = -2*((diff(x(t), t))/sqrt(2)-(diff(y(t), t))/sqrt(2))+diff(x(t), t), diff(y(t), t) = -2*(temp/sqrt(2)-(diff(y(t), t))/sqrt(2))+diff(y(t), t)]]]);
plots:-odeplot(res, [x(t), y(t)], 0 .. 100, axes = boxed, frames = 100)

Использованы формулы с http://www.gamedev.ru/code/forum/?id=14 ... e=4#m50%20