Кривая
|
Источник света
|
Каустика
|
Окружность
|
На плоскости
|
Кардиоида
|
Окружность
|
Не на плоскости
|
Улитка Паскаля
|
Окружность
|
Бесконечность
|
Нефроида
|
Парабола
|
Лучи,
перпендикулярные директрисе
|
Кривая третьего
порядка Чирнгаузена
|
Кривая
Чирнгаузена
|
Фокус
|
Полукубическая
парабола
|
Циссоида Диоклеса
|
Фокус
|
Кардиоида
|
Кардиоида
|
Касп
|
Нефроида
|
Четырёхлистник
|
Центр
|
Астроида
|
Дельтоид
|
Бесконечность
|
Астроида
|
Логарифмическая
спираль
|
Центр
|
Логарифмическая
спираль
|
Циклоида
|
Лучи,
перпендикулярные линии через каспы
|
1/2 Циклоиды
|
y = Ex
|
Лучи,
перпендикулярные y-лучу
|
Цепная линия
|
1.4
Математическая модель каустики
Теория огибающих имеет важное приложение в
оптике. Пусть точечный источник света расположен в точке A плоскости.
Лучи света, исходящие из A, падают на кривую Z, которая является
зеркалом, и, отражаются от нее по правилу равенства углов падения и отражения.
Огибающую всех отраженных лучей называют катакаустикой или фокальной
линией отражения.
Теорема: Пусть
точка Aимеет координаты (). Тогда
уравнение прямой, на которой лежит луч, отраженный от точки (X,Y) зеркала
Zимеет вид:
(1.4)
Где
§ N=0 - уравнение в
декартовых координатах (x,y) касательной к кривой zв точке (X,Y);
§ T=0 - уравнение в
декартовых координатах (x,y) нормали к кривой zв точке (X,Y);
§ и - результат подстановки в выражения для N и T.
Если зеркало Z задано параметрическими
уравнениями:
(1.5)
то уравнение
(1.6)
определяет прямую, на которой лежит луч,
отраженный от зеркала в точке (X,Y).
Для получения огибающей этого семейства
прямых нужно продифференцировать это уравнение по . Два получившихся
уравнения крайне громоздки и исключить из них переменную в общем случае
затруднительно. Тем не менее, можно сравнительно легко найти параметрическое
представление дискриминантной кривой. Дело в том, что оба уравнения являются линейными
по переменным и , и получившуюся систему
можно разрешить - например, по формулам Крамера.
Пример. Найти катакаустику при отражении от параболы и при источнике света, расположенном в точке x=0,y=0.
Рисунок 1.12 - Отражение лучей от кривой
На рисунке 1.12 падающие лучи показаны желтым, отраженные -
красным.
Решение. Имеем и уравнение семейства отраженных прямых
(1.7)
Для нахождения уравнения огибающей данного
семейства, вычисляем дискриминант этого полинома, рассматриваемого относительно
переменной t:
(1.8)
Это уравнение определяет одну изолированную точкуx=0, y=1/4 (фокус параболы) и некоторую кривую, симметричную
относительно оси . Мы выделим только ее кусок, лежащий внутри параболы (рис.1.13):
Рисунок 1.13 - Каустика
Теперь проиллюстрируем другой способ представления каустики - в
параметрическом виде. Продифференцируем уравнение семейства по параметру t:
(1.9)
и разрешим получившуюся систему
относительно x и y:
(1.10)
Заметим, что если бы мы двигали источник света по оси по направлению к фокусу параболы, то отраженные от нее лучи
пересекались бы все дальше и дальше, пока не стали бы параллельными этой оси.
На этом факте основано применение параболы: отражающая поверхность прожектора
делается в виде параболоида вращения, в фокус которого помещается источник
света.
Можно и "обратить" процесс, заставив лучи, приходящие
"из бесконечности" собираться в фокусе - тогда мы получаем
параболическую антенну спутникового телевидения, или же боевое оружие
античности, первое появление которого принято связывать с именем Архимеда.
Рисунок 1.14 - Катакаустика
Пример. Катакаустика (рис.1.14) при отражении от полуокружности
при источнике света, находящемся "на бесконечности":
дает представление о картинке, наблюдаемой при отражении солнечных
лучей от поверхности чашки (рис.1.15):
Рисунок 1.15 - Каустика в чашке
Получим теперь формулы для каустики. Для окружности и источнике света, находящемся в бесконечности на оси , уравнение семейства отраженных лучей имеет вид:
(1.11)
Дифференцируем это равенство по t и
находим y, затем подстановкой в исходное устанавливаем x:
(1.12)
Из этих уравнений можно найти и неявное
представление каустики в виде:
(1.13)
Эта кривая известна как нефроида и
относится к типу эпициклоид.
Пример. Катакаустика при отражении от
окружности при источнике света, находящемся внутри нее (рис.1.16):
Рисунок 1.16 - Катакаустика (источник
света внутри).
Для окружности и при имеем уравнение каустики в виде:
(1.14)
При источнике на самой окружности - например, при , - получим уравнение каустики в виде:
(1.15)
оно определяет замкнутую кривую, которая
называется кардиоидой (рис.1.17).
Рисунок 1.17 - Кардиодида
1.5 Примеры
каустик
Радуга - это частный случай каустики, игр света. Каустика
представляет собой сложную, а иногдаи очень красивую картину, создаваемую
сходящимися световыми лучами в результате их (многократных) преломлений и
отражений на поверхностях раздела сред с различной оптической плотностью.
Простейшими случаями каустики являются: яркая точка света в фокусе собирающей
линзы; тонкий луч прожектора, в фокусе параболического зеркала которого имеется
точечный источник света; сложная световая фигура (кардиоида) внутри полого,
открытого сверху цилиндра (например, чашки) при отражении света от его
внутренней поверхности (рис.1.18); дрожащие ячейки света на дне неглубокого
водоёма (рис.1.19); узкие лучи, полученные в результате отражений от
поверхностей2-го и более высоких порядков; солнечные и лунные дорожки на
поверхности воды и, наконец, радуга на небе.
Рис.1.18. Каустика в виде кардиоиды
Рис.1.19. Каустика на дне освещённого бассейна
Значение каустик во вселенной
В настоящее время распределение вещества во Вселенной
чрезвычайно неоднородно (существуют планеты. Солнце, звезды, галактики,
скопления галактик и т.д.). Современная астрофизика считает, что на ранних
этапах развития Вселенной таких неоднородностей не было. Как же они
сформировались? Я.Б. Зельдович в 1970 г. предложил объяснение образования
скоплений пылевидной материи, математически эквивалентное анализу возникновения
особенностей каустик, начатому в 1963 г.Е.М. Лифшицем, Халатниковым и
Судаковым.
Рассмотрим бесстолкновительную среду, т.е. среду настолько
разреженную, что ее частицы проходят друг "сквозь" друга, не
сталкиваясь. Предположим, для простоты, что частицы не взаимодействуют и
движутся по инерции: через какое-то время t частица, находившаяся в
точке х, перейдет в точку х + vt.
Предположим, что в начальный момент скорость частицы,
находящейся в точке х, была v0 (х); векторное поле v0 называется начальным
полем скоростей среды. С течением времени частицы будут двигаться, и поле
скоростей будет меняться (хотя скорость каждой частицы и не меняется, в
следующий момент времени эта частица находится на новом месте).
На рисунке1.20показано начальное поле скоростей одномерной
среды v0 и получающиеся из него через время t = 1, 2 и 3 ноля v1, v2, v3. Мы наблюдаем, что,
начиная с определенного момента, более быстрые частицы начинают обгонять более
медленные; в результате чего, поле скоростей становится трехзначным: через одну
точку пространства три потока частиц проходят с разной скоростью.
Рисунок 1.20 - Эволюция поля скоростей
бесстолкновительной среды
Движение нашей среды можно описать как однопараметрическое семейство
отображений прямой на прямую. Собственно для каждого t определено отображение gt, переводящее начальное
положение частицы (х) в конечное: gt = х + v0 (х) t.
Отображение g0 является тождественным преобразованием,
оставляющим каждую точку на месте. Отображения, соответствующие моментам t, близким к 0, взаимно
однозначны и не имеют особенностей. После момента первого обгона отображение gt имеет две складки.
Рисунок 1.21 - Особенности плотности после обгона
Пусть в начальный момент плотность среды в точке х была ρ0 (х), Современем плотность изменится. Легко догадаться, что после
обгона график плотности будет иметь вид (на расстоянии 8 от точки складки
плотность оказывается порядка 1/√ε).
Исходя из этого, небольшие отличия начального ноля скоростей
от постоянного приводят через достаточно большое время к образованию скоплений
частиц (в местах бесконечно большой плотности).
Глава II.
Maple
Для построения каустик, используется система компьютерной
математикиMaple (ver.17) о которой и будет идти речь в этой главе.
2.1 Общие
сведения о СКМ Maple
Maple-это система компьютерной математики созданная
для символьных вычислений, работы с математическими выражениями в разных
формах.
Система компьютерной математики является мощным инструментом
для выполнения вычислений и работы с разными формами математических и
аналитических расчетов. Данная система может быть применена в образовательных и
научных целях.
2.2
Программирование в СКМ Maple
Оператор цикла
Как и в другихязыках программмирования, в операторе цикла
надо указать переменную цикла, ее начальное значение (from), шаг (by), и ее конечное значение
(to) и условие для ее
выполнения (while). Тело цикла, если даже тело состоит не из большого количества
операторов, располагается между do и od. Команда для закрытия скобки odможет быть убрана и
вместо нее могут стоять команды end или enddo. Не все команды оператора цикла могут быть
включены в цикл, это является не обязательным условием. Обычно за начальное
условие применяется единица и тогда это выглядит следующим образом from: =1 (это можно не
писать). Как и начальное условие, шаг цикла будет равен единице,a условие while, и, если дано конченое
число для заданной переменной данного цикла, нужно указывать только для того
чтобы досрочно выйти из цикла. Аналогично, если будет дано условие, то конечное
число можно не вносить. Пример:
>k: =0;
> for i from 1 by 2 to 12 while
(k<3) do k: =i od;
При выполнении данного цикла (в данном случае он выйдет из
цикла при значении k равным 4) на экран выводиться всего две записиk: =1 и k: =3, фиксирующие
присвоение заданной переменной k чисел 1 и 3.
До выполнения данного цикла нашей переменной следует
присвоить любое значение, если мы этого не сделаем система выдаст сообщение об
ошибке и не сможет выполнить проверку этого условия (Error, cannot determine if this expression is true or false). Наша переменная цикла
после выполнения присваивает и сохраняет последнее пройденное значение. У нас в
примере им будет число 5. Чтобы программа нам вывела на экран значение i не нужно использовать
обычные операторы print или while. Нужно написать имя нашей переменной и точку с запятой:
> i;
Если вместо точки с запятой поставить двоеточие, то на экран
значение переменной выводиться не будет. А если цикл двойной,
например,
> for i to 5 do for j to 6 do i*j od;
od;
То в независимости от того что мы поставим: точку с запятой
или же двоеточие, результат промежуточных вычислений мы не сможем увидеть, пока
не применим специальный оператор вывода, например,print.
Условный оператор
Простейший условный оператор включает три основные служебные
команды: if, then и end. Логическое выражение (условие) ставиться после оператора if. Всегда в конце
условного оператора следует поставить: end,fi или endif. Например: if k<5 then k5end. Если взять усложненную
форму оператора, то есть возникает возможность выбора из нескольких вариантов: if k <5 then k5else k0 end. При выполнении условия k <5, на экран будет
выводиться значение переменной k5, если оно не будет выводиться, то на экран выводиться
условие k0.
Также условный оператор может выполнять команды как оператор выбора:
>ifk=1 thenk1
>elif k=2 then k2
>elif k=4 then k4
>elif k=5 then k5 else k0 end if;
Опираясь на значение переменной kпроизводиться вывод на
экран одного и заданных выражений k1, k2, k3 или k5. Но если не удовлетворяется ни одно из условий, то оператор
выдаст значение k0. Есть еще условный оператор, который имеет следующую структуру
обращений: ‘if’ (K, M, L). Тут K - заданное условие, за Mвыступает выражения, возвращающее оператор. Ели
заданное условие является истинным, за Lвыступает выражение,
возвращающее условие, если данное условие ложно. Данный оператор удобно
использовать определенно в формулы, в качестве некоторой функции какой - либо
переменной K
с обоими ее значениями Mи L.
Процедуры
Для описания простейшей процедуры нужно описать ее имя,
перечислить аргументы и задать тело процедуры. В ниже приведенной процедуре FcF с аргументами x и y производиться вычисление
произведения xy:
> F: =proc (x,y) x*y end;
Для того чтобы обратиться к процедуре может выглядеть
следующим образом: F (5,9).
Функции
В системе компьютерной математики Mapleфункция является
упрощенным типом процедуры. Как и в процедурах нужно указывать название,
обозначать аргументы и тело заданной функции. Все аргументы нужно указывать в
скобках. Если аргумент единственный, то скобки можно не писать. Ниже приведен
пример, помогающий вычислить мощность объединения и пересечения множеств:
> A: ={a,b,c}:
> B: ={a,b,d}:
> f: = (x,y) - > (nops (x union
y),nops (x intersect y));
> f (A,B);: = (x, y) → (nops (x ∪ y), nops (x ∩
y))
4, 2.
2.3 Работа с
графикой
Для того чтобы построить график или нарисовать рисунок,
включая анимированных графиков и рисунков, в СКMMaple существует большое
количество разнообразных команд. Множество из них требует того, чтобы был
вызван пакет, включающий данную команду. Несложные графики можно построить
тривиальным оператором plot. Для оператора plotне нужно подключение
каких-либо команд. Для него не нужно дополнительных средств.
Приведем некоторые простые варианты для построения графиков
функций.
График элементарных функций; plot (cos). В данной команде
указание аргумента не обязательно. (рис.2.1).
Например, такие не очень привычные выражения, как plot (sin+cos). Такие математические
высказывания случат либо как демонстративный материал, либо как не очень
сложный справочный. Довольно нечасто такие виды обращений используют для
оператора plot.
График (рис.2.2) функции y = sin (x) с обозначением пределов x=0.2*Piи обрезанием системы
координат (0.0.8):
plot (sin (x),x=0.2*Pi,0.0.8).
Если не фиксировать ограничения по оси, то на графике будет
показан весь график (если он не определен на бесконечности).
Рисунок 2.1-График функций y = sin (x) Рисунок2.2 - График функций y = sin (x) с ограничениями
Глава III
3.1
Построение модели каустики для кривой второго порядка, заданной в
параметрическом виде
В данной работе рассматривается два случая:
) Когда источником света выступают точки находящиеся на
заданном расстоянии друг от друга на заданной прямой
) Когда источником света выступает пучок света.
Для каждого из этих случаев рассматриваются следующие модели,
состоящие из:
. Источника света (Пучок или параллельные лучи)
. Отражающей поверхности (кривая) заданной в
параметрическом виде.
. Лучей света, описываемых законами геометрической
оптики.
Модель источников света
Для начала рассмотрим случай, когда источником света
выступает прямая на которой расположены точки, на прямой параллельной Oy (x=const) и отдалённой от вершины
заданной кривой на расстоянии L. Источники света равномерно располагаются на заданной прямой
с шагом h.
Все источники на заданной прямой излучают свет в одинаковом направлении, и под
определенным углом к оси абсцисс.
Рисунок 3.1 - Построение отражения луча
В качестве источника света (рис.3.1) возьмем прямую (1), из
которой будут выходить лучи света, заданные через направляющий вектор (2). В
работе я буду рассматривать различные вариации направления вектора и начальных
данных прямой. Исследую и рассмотрю различные прямые и различные источники
света, а также различное количество лучей падающий на заданную кривую.
Геометрическая оптика
Геометрическая оптика - это раздел оптики, рассматривающий
законы излучения света в различных средах и отражения света от различного вида
поверхностей.
Закон отражения света
Луч света в определенной среде прямолинеен до того момента,
пока он не достигнет границы данной среды с другой средой. На границе этих сред
луч светаизменяет свое направление. Часть света (а в ряде случаев и весь свет)
возвращается в первую среду. Это явление называется отражением света.
Одновременно свет частично проходит во вторую среду, меняя при этом направление
своего распространения - преломляется.
Вычисление отражённого
луча
Если использовать вектора, то решение
находиться в более простом виде (рис.3.2):
Рисунок 3.2 - Геометрическое представление отражения
луча
Для того чтобы найти вектор отражения нам нужно из вектора
направления света вычесть удвоенное произведение вектора нормали на отношение
скалярного произведения вектора направления света и нормали к скалярному
произведению вектора нормали самого на себя:
(3.1)
Делит на скалярное произведение нормали самой на себя не
нужно, если данная нормаль уже нормирована.
Если прямая задана точкой (, ) и вектором направления (, ), то вектор нормали (ненормированный). В некоторых случаях нужно
знать знак нормали, чтобы знать, какая сторона прямой "внешняя".
Алгоритм построения для первого источника света (параллельные
лучи) (Рис.3.3):
1. Построение кривой. Кривая должна быть задана в параметрическом
виде.
. Задаем прямую на которой будут располагаться источники света.
. Задаем точку на данной прямой из которой будет исходить луч
света.
. Задаем направление вектора источника света.
. При помощи точки и направляющего вектора вычисляем уравнение
прямой.
. Строим отрезок соединяющий точку на прямой и точку пересечения.
. Находим нормаль.
. Строим направление отраженной от кривой луча.
. Строим прямую, определяющую направление отраженного луча.
Рисунок 3.3 - Параллельно направленные лучи
Рисунок 3.4 - Пучок света
Алгоритм построения для второго источника света
(Пучок) (Рис.3.4):
1. Построение кривой. Кривая должна быть задана в
параметрическом виде.
. Задаем точку источники света.
. Задаем прямую на которой будут располагаться координаты
направляющих векторов.
. Задаем направление вектора источника света.
. При помощи точки и направляющего вектора вычисляем
уравнение прямой.
. Строим отрезок соединяющий точку на прямой и точку
пересечения.
. Находим нормаль.
. Строим направление отраженной от кривой луча 9. Строим
прямую, определяющую направление отраженного луча
В данном примере будем рассматривать в качестве произвольной
кривой кривую второго порядка параболу.
Ее уравнение в параметрическом виде будет выглядеть следующим
образом:
, (3.2)
где t-параметр.
Уравнение прямой будет выведено при помощи начальной точки и
вектора направления:
, (3.3)
где k-параметр, - координата начальной точки на оси Ox, -координата начальной точки на оси Oy, -координата направления вектора на оси Ox,-координата направления вектора на оси Oy.
Для нахождения точек пересечение запишем наше уравнение в
векторном виде:
, (3.4)
Затем нужно найти точки пересечения, точки пересечение получим
приравняв по координатам части пар метрически заданных функций:
, (3.5)
При решении данной системы получим коэффициенты,,,. Далее нам нужно выбрать подходящие нам. Выбирать будем при
помощи вычисления расстояния от точки на прямой источников света до каждой из
полученных пары решений.
3.2
Построение модели каустики для плоскости, заданной в параметрическом виде
Построение каустик в пространстве производиться подобно
двумерному случаю.
Алгоритм
1. Строим геометрические фигуры, с определенными
ограничениями.
2. Задаем плоскость, которая будет выступать в роли
источника света.
. Задаем точку на плоскости, из которой будет исходить
луч света
. Задаем направление вектора при помощи второй точки.
. Через начальную точку, и направление вектора находим
точку пересечения луча с плоскостью.
. Строим отрезок соединяющий точку на прямой и точку
пересечения.
7. Находим нормаль.
8. Строим направление отраженного от плоскости луча.
. Строим прямую, определяющую направление отраженного
луча.
3.3 Написание
программных процедур СКМ Maple Двумерный случай
Для начала задаем параметрическое уравнение кривой второго
порядка:
r: = [t^2/ (2*3),t]:
Далее вводим координаты начальной точки:
x [0]: =10;
y [0]: =5;
Задаем направляющий вектор:
x_dir [0]: =-1;
y_dir [0]: =.3;
Выводим уравнение прямой, с учетом того вектор направления
может быть равен нулю:
ifx_dir [0] =0 then
a: = [x [0],y_dir [0] *k+y [0]]; else: =
[x_dir [0] *k+x [0],y_dir [0] *k+y [0]]::
Описываем уравнение наше кривой и ограничения на которых она
будет строиться:
plt [01]: =plot ([op
(r),t=-20.20],scaling=constrained):
Описываем точку, при помощи которой будем задавать
направление
луча:
plt [02]: =point ([x [0],y [0]]):
Описываем точку, которая будет задавать направления луча
света
выводим все значения, которые мы записали:
plt [03]: =point ([x [0] +x_dir [0],y [0]
+y_dir [0]]):
Описываем уравнение нашей прямой и ограничения на которых она
будет задаваться:
plt [04]: =plot ([op
(a),k=-20.20],scaling=constrained):
Выводим все значения, которые мы получили (Рис.3.4):
display ([plt [01],plt [02],plt [03],plt
[04]],scaling=constrained,axes=none);
Рисунок 3.4 - Построение кривой, начальной точки, и
вектора направления
Следующим шагом записываем систему уравнений, для нахождения
точки падения луча:
tmp1: =solve (eq,{t,k});
Находим точки пересечения прямой и кривой:
y [1]: =rhs (tmp1 [1] [2]):[1]: =subs
(t=rhs (tmp1 [1] [2]),r [1]):[2]: =rhs (tmp1 [2] [2]):[2]: =subs (t=rhs (tmp1
[2] [2]),r [1]):
Рассчитываем расстояние между корнями получившихся уравнений
и начальной точкой:
per1: =sqrt ( (x [1] - x [0]) ^2+ (y [1] -
y [0]) ^2):: =sqrt ( (x [2] - x [0]) ^2+ (y [2] - y [0]) ^2):
Выбираем подходящую нам точку при помощи задания условия:
if evalf (sqrt ( (x [1] - x [0]) ^2+ (y
[1] - y [0]) ^2)) < (sqrt ( (x [2] - x [0]) ^2+ (y [2] - y [0]) ^2))
then[11]: =x [1]x [11]: =x [2]if;evalf (sqrt ( (x [1] - x [0]) ^2+ (y [1] - y
[0]) ^2)) < (sqrt ( (x [2] - x [0]) ^2+ (y [2] - y [0]) ^2)) then[11]: =y
[1]y [11]: =y [2]if;
Для вычисления касательной находим производную:
dr: =diff (r,t);
Находим касательную:
kas: =evalf (subs (t=y [11],dr));
Задаем условие для вывода нормали и выводим нормаль:
norm1: = [kas [2],-kas [1]];norm1 [1] * (x
[0] - x [11]) +norm1 [2] * (y [0] - y [11]) <0norm1: = [-kas [2],kas
[1]];if;
Описываемточку определяющую вектор нормали:
plt [07]: =point ([x [11] +norm1 [1],y [11] +norm1 [2]]):
Выводим все значения, которые мы получили (Рис.3.5):
display ([plt [01],plt [02],plt [04],plt
[07]],scaling=constrained,axes=none);
Рисунок 3.5 - Построение кривой, начальной точки, точки
пересечения и направления нормали
По ранее выведенным формулам находим поэтапно вектор
отражения.
Первым этапом находим отношение скалярного произведения
вектора направления и вектору нормали к скалярному произведению вектора нормали
самого на себя:
n1: = (x_dir [0] *norm1 [1] +y_dir [0]
*norm1 [2]) / (norm1 [1] *norm1 [1] +norm1 [2] *norm1 [2]):
Затем находим удвоенное произведение вектора нормали на
получившуюся константу n1:
n2: = [norm1 [1] *2*n1,norm1 [2] *2*n1]:
Вычитаем из вектора направления вектор n2:
l: = [x_dir [0] - n2 [1],y_dir [0] - n2
[2]]:
Получившийся вектор l и будет вектором отражения. Запишем уравнение
отражённой прямой:
lk: =solve ( (x-x [11]) /l [1] = (y-y [11]) /l [2],y);
Описываем уравнение отраженной прямой:
plt [08]: =plot (lk,x=x [11] - 5. x [11],scaling=constrained):
Выводим все значения которые мы получили (Рис.3.6):
display ([plt [01],plt [02],plt [04],plt
[05],plt [08]],scaling=constrained,axes=none);
Рисунок 3.6 - Построение кривой, начальной точки, и
вектора направления
Далее пишем процедуру для некоторого количества лучей.
Для начала зададим входные параметры: кривую в
параметрическом виде и ограничения для нее:
r1: = [t^2/ (2*60),t]:
t1: =-10:
t2: =10:
1) Затем перейдем к первой процедуре.
Первая процедура будет заполнять заданный нами массив
определенными значениями.
В данной процедуре будут заполняться только две точки, первая
из которых является начальной, а вторая характеризует вектор направления:
prbRefLines: =proc (L,h,d1,d2,N)Line_,
r1;dh,dir, i;: = [d1,d2]:: = (h/N):i from 1 to N do_ [i]: = [[L,-h/2+ (i-1)
*dh],dir, [0,0], [0,0]];:
endproc:
2) Вторая процедура будет считать точку пересечения и точку
характеризующую вектор направления луча отражения:
prbRefArray: =proc (L,N)Line_,
r1:dir_,x,y,a_,tmp1,kas,per1,per2,norm,l,n1,n2,k, i,dr;
# вектор направленияi
from 1 to N do
# расчет пересечения:_:
= [k*Line_ [i] [2] [1] +Line_ [i] [1] [1],k*Line_ [i] [2] [2] +Line_ [i] [1]
[2]]:Line_ [i] [2] [1] =0 then_: =x=Line_ [i] [1] [1]:_: = [k*Line_ [i] [2] [1]
+Line_ [i] [1] [1],k*Line_ [i] [2] [2] +Line_ [i] [1] [2]]:if:({r1 [1] =a_
[1],r1 [2] =a_ [2] },{t,k}):: =solve ({r1 [1] =a_ [1],r1 [2] =a_ [2]
},{t,k}):[1]: =rhs (tmp1 [1] [2]):[1]: =subs (t=rhs (tmp1 [1] [2]),r1 [1]):[2]:
=rhs (tmp1 [2] [2]):[2]: =subs (t=rhs (tmp1 [2] [2]),r1 [1]):: =evalf (sqrt (
(x [1] - Line_ [i] [1] [1]) ^2+ (y [1] - Line_ [i] [1] [2]) ^2)):: =evalf (sqrt
( (x [2] - Line_ [i] [1] [1]) ^2+ (y [2] - Line_ [i] [1] [2]) ^2)):evalf (sqrt
( (x [1] - Line_ [i] [1] [1]) ^2+ (y [1] - Line_ [i] [1] [2]) ^2)) <evalf
(sqrt ( (x [2] - Line_ [i] [1] [1]) ^2+ (y [2] - Line_ [i] [1] [2]) ^2)) then_
[i] [3] [1]: =x [1]Line_ [i] [3] [1]: =x [2]if:evalf (sqrt ( (x [1] - Line_ [i]
[1] [1]) ^2+ (y [1] - Line_ [i] [1] [2]) ^2)) <evalf (sqrt ( (x [2] - Line_
[i] [1] [1]) ^2+ (y [2] - Line_ [i] [1] [2]) ^2)) then_ [i] [3] [2]: =y
[1]Line_ [i] [3] [2]: =y [2]:
# Нормаль в новой точке
dr: =diff (r1,t):
kas: =evalf (subs (t=Line_ [i] [3]
[2],dr)):: = [kas [2],-kas [1]]:evalf (norm [1] * (Line_ [i] [1] [1] - Line_
[i] [2] [1]) +norm [2] * (Line_ [i] [1] [2] - Line_ [i] [2] [2])) <0 then
norm: = [-kas [2],kas [1]]:if:
#Вектор направления отражения: =
(Line_ [i] [2] [1] *norm [1] +Line_ [i] [2] [2] *norm [2]) / (norm [1] *norm
[1] +norm [2] *norm [2]):: = [norm [1] *2*n1,norm [2] *2*n1]:: = [Line_ [i] [2]
[1] - n2 [1],Line_ [i] [2] [2] - n2 [2]]:
#точка отраженная
Line_ [i] [4] [1]: =l [1]:
Line_ [i] [4] [2]: =l [2]:do::
) Третья процедура будет демонстрировать лучи падения и лучи
отражения:
pprbLines: =proc
(L,h,N)r1,t1,t2:pltprb,pltLines,pltLines1;(r1);: =plot ([op (r1),t=t1.
t2],scaling=constrained,=black,=3):: =display (((_ [i] [1],_ [i] [3],=red),=1.
N),=false,=constrained
);: =display (((_ [i] [3],_ [i]
[4],=green),=1. N),=false,=constrained
);([pltLines,pltprb,pltLines1],scaling=constrained,axes=none);:
И самая главная процедура которая будет включать три
предыдущие и выводить наши отражения:
Mainproc: =proc (L,h,d1,d2,N)
global r1,t1,t2:(L,h,d1,d2,N):(L,N):(L,h,N):proc:(20,
20,-10,-.1,10);(Line_):
В результате получим следующую картинку (Рис.3.7).
Рисунок 3.7 - Каустика, при отражении
лучей от параболы (параллельно падающие)
Аналогично записывается написание процедур
для пучка света с поправками начальные условия. Процедуру будем задавать в
более удобной для нас форме. Процедура будет включать в себя три под процедуры,
у каждой из которых будет своя функция (рис.3.8).
Рисунок 3.8 - Каустика, при отражении
лучей от параболы (пучок)
Трехмерный случай
Отражение на плоскости будем рассматривать
на примере параболоида.
Задаем начальные условия для прямой:
x [0]: =0:
y [0]: =0:
z [0]: =0:
ax: =1:
ay: =-.3:
az: =10:
Задаем уравнение плоскости и уравнения прямой:
a: =x^2+y^2-z; b: =x-x [0] /ax=y-y [0]
/ay; c: =y-y [0] /ay= (z-z [0]) /az;
Рисуем прямую и параболоид:
pl_b: =spacecurve ([ax*t+x [0],ay*t+y
[0],az*t+z [0]],t=-10.10):_a: =implicitplot3d
(a,x=-10.10,y=-10.10,z=-10.10,numpoints=5000):
Находим точку пересечения прямой и плоскости:
sol: =evalf (solve ({a,b,c},{x,y,z}));_:
=evalf (subs (sol,x));_: =evalf (subs (sol,y));_: =evalf (subs (sol,z));_b:
=spacecurve ([ax*t+x [0],ay*t+y [0],az*t+z [0]],t=-1.1):_a: =implicitplot3d
(a,x=-10.10,y=-10.10,z=-10.10,numpoints=5000):_c: =pointplot3d ([x_,y_,z_],color=red);
Находим производные в точках:
x1: =subs ({x=x_,y=y_,z=z_},diff (a,x)):
y1: =subs ({x=x_,y=y_,z=z_},diff (a,y))::
=subs ({x=x_,y=y_,z=z_},diff (a,z)):
Вектор прямой задается следующим образом:
l: = [ax,ay,az];
Вектор нормали задается следующим образом:
n1: = [x1,y1,z1]; 4
Уравнение нормали и вывод:
pl_e: =spacecurve
([x1*t+x_,y1*t+y_,z1*t+z_],t=-8.4,color=green):
Аналогично двумерному случаю находим вектор отражения при
помощи вектора падения луча и вектора нормали:
k: =dotprod (l,n1) /dotprod (n1,n1);: =
[2*k*n1 [1],2*k*n1 [2],2*k*n1 [3]];
r: =l-p;
Вектор lбудет вектором отражения.
Задаемуравнениеотраженноголуча, и
обозначаем его:
b1: = (x-x_) /r [1] = (y-y_) /r [2];: =
(y-y_) /r [2] = (z-z_) /r [3];
pl_d: =spacecurve ([r [1]
*t+x_,r [2] *t+y_,r [3] *t+z_],t=-1.1,color=yellow):
Выводим луч падения (фиолетовый), нормаль (зеленая),
отраженный луч (желтый) (Рис.3.8).
Рисунок 3.8 - Отражение луча от
плоскости.
Заключение
В ходе данной работы был рассмотрен
теоретический материал по теме каустик и были изучены их свойства. Также в ходе
магистерской диссертации былаизучена среда СКМ Maple. Была создана процедура
для моделирования каустик для кривых второго порядка.
Таким образом, цели, поставленные в
квалификационной работе, выполнены.
Списокиспользованнойлитературы
1. Burridge,R. Asymptoticevaluationofintegrals related to
time-domain fields near caustics, SIAM J. Appl. Math., 55: 2, 1995.
2. Андреев,
А.Н. Каустики на плоскости/ А.Н. Андреев, А.А. Панов "Квант" - 2010.
- №3.
. Арнольд,
В.И. Особенности каустики волновых фронтов, В.И. Арнольд - М.: Фазис, 1996.
. Арнольд,
В.И. Теория катастроф / Арнольд В.И. М.: Наука, 1990.
. Баев,
А.В. Математическое моделирование волн в слоистых средах вблизи каустики/ А.В.
Баев // Матем. Моделирование-2013 - 25: 12-С.83-102
6. Гольдин,
С.В. Сейсмическое волновое поле в близи каустик: анализ во временной области/
С.В. Гольдин, А.А. Дучков // Физика Земли. - 2002. - №7. С.56-66.
7. Дьяконов,
В.П. Maple 9.5/10 в математике, физике и образовании/ В.П. Дьяконов - М.:
СОЛОН-Пресс, 2006. - 720 с.
. Кирсанов,
М.Н. Задачи по теоретической механике с решениями в Maple 11. М.: Физматлит,
2010, 264с.
. Кирсанов, М.Н. Maple и Maplet. Решения задач механики /
Кирсанов М.Н. - М.: Лань, 2012. - 512с.
. Лонгейр,
М. Крупномасштабная структура Вселенной/ Ред.М. Лонгейр, Я. Эйнасто. - М.: Мир,
1981.
. Постон,
Т. Теория катастроф и ее применения/ Постон Т., Стюарт И. - М.: Мир, 1980.
. Яновская,
Т.Б. Численный метод расчета поля поверхностной волны при наличии каустик /
Яновская Т.Б., Гейгер М.А. // Физика Земли. - 2007.
Приложения
Приложение А.
Программные процедуры моделирования каустик
#Отражение от произвольной кривой (параллельное
падение лучей)
>restart:(plots):(plottools):
>prbRefLines: =proc
(L,h,d1,d2,N)Line_;dh,dir, i;: = [d1,d2]:: = (h/N):i from 1 to N do_ [i]: =
[[L,-h/2+ (i-1) *dh],dir, [0,0], [0,0]];do:proc:
>prbRefLines (5,10,-1,-.1,50);(Line_):
>prbRefArray: =proc
(p,L,N)Line_:dir_,x,y,a_,r1,tmp1,kas,per1,per2,norm,l,n1,n2,k, i,dr;
#print (Line_ [1]);
# векторнаправленияi from 1 to N do
#print (Line_ [i] [2] [1]):_ [i] [2] [2]:
# расчетпересечения:_: = [k*Line_ [i] [2]
[1] +Line_ [i] [1] [1],k*Line_ [i] [2] [2] +Line_ [i] [1] [2]]:Line_ [i] [2]
[1] =0 then_: =x=Line_ [i] [1] [1]:_: = [k*Line_ [i] [2] [1] +Line_ [i] [1]
[1],k*Line_ [i] [2] [2] +Line_ [i] [1] [2]]:if:
#print (a_);: = [t^2/ (2*p),t]:({r1 [1]
=a_ [1],r1 [2] =a_ [2] },{t,k}):: =solve ({r1 [1] =a_ [1],r1 [2] =a_ [2] },{t,k});[1]:
=rhs (tmp1 [1] [2]):[1]: =subs (t=rhs (tmp1 [1] [2]),r1 [1]):[2]: =rhs (tmp1
[2] [2]):[2]: =subs (t=rhs (tmp1 [2] [2]),r1 [1]):: =sqrt ( (x [1] - Line_ [i]
[1] [1]) ^2+ (y [1] - Line_ [i] [1] [2]) ^2):: =sqrt ( (x [2] - Line_ [i] [1]
[1]) ^2+ (y [2] - Line_ [i] [1] [2]) ^2):
#print ([per1,per2]):evalf (sqrt ( (x [1]
- Line_ [i] [1] [1]) ^2+ (y [1] - Line_ [i] [1] [2]) ^2)) <evalf (sqrt ( (x
[2] - Line_ [i] [1] [1]) ^2+ (y [2] - Line_ [i] [1] [2]) ^2)) then_ [i] [3]
[1]: =x [1]Line_ [i] [3] [1]: =x [2]if:evalf (sqrt ( (x [1] - Line_ [i] [1]
[1]) ^2+ (y [1] - Line_ [i] [1] [2]) ^2)) <evalf (sqrt ( (x [2] - Line_ [i]
[1] [1]) ^2+ (y [2] - Line_ [i] [1] [2]) ^2)) then_ [i] [3] [2]: =y [1]Line_
[i] [3] [2]: =y [2]:
# Нормаль в новой точке
dr: =diff (r1,t):
kas: =evalf (subs (t=Line_ [i] [3]
[2],dr)):: = [kas [2],-kas [1]]:evalf (norm [1] * (Line_ [i] [1] [1] - Line_
[i] [2] [1]) +norm [2] * (Line_ [i] [1] [2] - Line_ [i] [2] [2])) <0 then
norm: = [-kas [2],kas [1]]:;
#norm1:
#Вектор направления отражения: = (Line_
[i] [2] [1] *norm [1] +Line_ [i] [2] [2] *norm [2]) / (norm [1] *norm [1] +norm
[2] *norm [2]):: = [norm [1] *2*n1,norm [2] *2*n1]:: = [Line_ [i] [2] [1] - n2
[1],Line_ [i] [2] [2] - n2 [2]]:
#точка отраженная
Line_ [i] [4] [1]: =l [1]:
Line_ [i] [4] [2]: =l [2]:do:proc:
>prbRefArray (60,10,50):
>pprbLines: =proc
(p,L,h,N)pltprb,pltLines,pltLines1,r1;
#: = [t^2/ (2*p),t]:(r1);: =plot ([op
(r1),t=-10.10],scaling=constrained,=black,=3):
#: =display (((_ [i] [1],_ [i]
[3],=red),=1. N),=false,=constrained
);: =display (((_ [i] [3],_ [i]
[4],=red),=1. N),=false,=constrained
);([pltLines,pltprb,pltLines1],scaling=constrained,axes=none);:
>pprbLines (60, 200,10,10);
#Отражение от произвольной кривой (пучок лучей)
> restart:(plots):(plottools):
> r1: = [t^2/ (200),t]:: =-20:: =20:
>prbRefLines: =proc (L,h,N)Line_,
r1;dh,dir, i,dL;: =L/2;: = (h/N):i from 1 to N do_ [i]: = [[L,0], [L-dL-L,-h/2+
(i-1) *dh], [0,0], [0,0]];do:proc:
>prbRefLines (5,10,5);
>prbRefArray: =proc (L,N)Line_, r1:dir_,x,y,a_,tmp1,kas,per1,per2,norm,l,n1,n2,k,
i,dr;
# векторнаправленияi from 1 to N do
# расчетпересечения:_: = [k*Line_ [i] [2]
[1] +Line_ [i] [1] [1],k*Line_ [i] [2] [2] +Line_ [i] [1] [2]]:Line_ [i] [2]
[1] =0 then_: =x=Line_ [i] [1] [1]:_: = [k*Line_ [i] [2] [1] +Line_ [i] [1]
[1],k*Line_ [i] [2] [2] +Line_ [i] [1] [2]]:if:({r1 [1] =a_ [1],r1 [2] =a_ [2]
},{t,k}):: =evalf (solve ({r1 [1] =a_ [1],r1 [2] =a_ [2] },{t,k})):_ [i] [3]
[1]: =subs (k=rhs (tmp1 [1]),a_ [1]):_ [i] [3] [2]: =subs (k=rhs (tmp1 [1]),a_
[2]):: =diff (r1,t):: =evalf (subs (t=Line_ [i] [3] [2],dr)):: = [kas [2],-kas
[1]]:evalf (norm [1] * (Line_ [i] [1] [1] - Line_ [i] [2] [1]) +norm [2] *
(Line_ [i] [1] [2] - Line_ [i] [2] [2])) <0 then norm: = [-kas [2],kas
[1]]:if:
#Векторнаправленияотражения: = (Line_ [i]
[2] [1] *norm [1] +Line_ [i] [2] [2] *norm [2]) / (norm [1] *norm [1] +norm [2]
*norm [2]):: = [norm [1] *2*n1,norm [2] *2*n1]:: = [Line_ [i] [2] [1] - n2
[1],Line_ [i] [2] [2] - n2 [2]]:
#точка отраженная
Line_ [i] [4] [1]: =l [1]:
Line_ [i] [4] [2]: =l [2]:do:proc:
>prbRefArray (10,5):
>pprbLines: =proc
(L,h,N)r1,t1,t2:pltprb,pltLines,pltLines1;: =plot ([op (r1),t=t1.
t2],scaling=constrained,=black,=3):: =display (((_ [i] [1],_ [i] [3],=red),=1.
N),=false,=constrained
);: =display (((_ [i] [3],
*Line_ [i] [4],=green),=1.
N),=false,=constrained
);([pltLines,pltprb,pltLines1],scaling=constrained,axes=none);proc:
>pprbLines (1, 20,1);
>pprbLines (200,.5,2):
>Mainproc: =proc
(L,h,N)r1,t1,t2:(L,h,N):(L,N):(L,h,N)::
>Mainproc (100, 20,10);
#Отражение от параметрически заданной плоскости
restart:(plots):(plottools):(linalg):[0]:
=0:[0]: =0:[0]: =0:: =1:: =-.3:: =10:: =x^2+y^2-z;: =x-x [0] /ax=y-y [0] /ay;:
=y-y [0] /ay= (z-z [0]) /az;(a,y);_b: =spacecurve ([ax*t+x [0],ay*t+y
[0],az*t+z [0]],t=-10.10,color=black):_a: =implicitplot3d
(a,x=-10.10,y=-10.10,z=-10.10,numpoints=5000):([pl_a,pl_b]);: =evalf (solve
({a,b,c},{x,y,z}));_: =evalf (subs (sol,x));_: =evalf (subs (sol,y));_: =evalf
(subs (sol,z));_b: =spacecurve ([ax*t+x [0],ay*t+y [0],az*t+z [0]],t=-1.1):_a:
=implicitplot3d (a,x=-10.10,y=-10.10,z=-10.10,numpoints=5000):_c: =pointplot3d
([x_,y_,z_],color=red);([pl_a,pl_b,pl_c]);: =subs ({x=x_,y=y_,z=z_},diff
(a,x)):: =subs ({x=x_,y=y_,z=z_},diff (a,y)):: =subs ({x=x_,y=y_,z=z_},diff (a,z))::
= [ax,ay,az];: = [x1,y1,z1];_e: =spacecurve
([x1*t+x_,y1*t+y_,z1*t+z_],t=-8.4,color=green):[1];: =dotprod (l,n1) /dotprod
(n1,n1);: = [2*k*n1 [1],2*k*n1 [2],2*k*n1 [3]];: =l-p;: = (x-x_) /r [1] =
(y-y_) /r [2];: = (y-y_) /r [2] = (z-z_) /r [3];_d: =spacecurve ([r [1] *t+x_,r
[2] *t+y_,r [3] *t+z_],t=-1.1,color=yellow):([pl_a,pl_b,pl_d,pl_e]);
#Частные случаи изображения каустики >restart:
with (plots):(plottools):
> r: = [u,u^2];[2];
>pl_a: =plot ([op
(r),u=-4.4],scaling=constrained);
> cx: =r [1] - (diff (r [2],u)) * (
(diff (r [1],u)) ^2+ (diff (r [2],u)) ^2) / ( (diff (r [1],u)) * (diff (r
[2],u$2)) - (diff (r [2],u)) * (diff (r [1],u$2)));: =r [2] + (diff (r [1],u))
* ( (diff (r [1],u)) ^2+ (diff (r [2],u)) ^2) / ( (diff (r [1],u)) * (diff (r
[2],u$2)) - (diff (r [2],u)) * (diff (r [1],u$2)));
> cx: =r [1] - (diff (r [2],u)) * (
(diff (r [1],u)) ^2+ (diff (r [2],u)) ^2) / ( (diff (r [1],u)) * (diff (r
[2],u$2)) - (diff (r [2],u)) * (diff (r [1],u$2)));: =r [2] + (diff (r [1],u))
* ( (diff (r [1],u)) ^2+ (diff (r [2],u)) ^2) / ( (diff (r [1],u)) * (diff (r
[2],u$2)) - (diff (r [2],u)) * (diff (r [1],u$2)));
> l: = [cx,cy];
>pl_b: =plot ([op
(l),u=-0.7.0.7],scaling=constrained,color=blue);
> display ([pl_a,pl_b]);
> restart:(plots):(plottools):
> a: =2;: =3;
> r: = [a*cos (u),b*sin (u)];[2];
>pl_a: =plot ([op
(r),u=1/2*Pi.3/2*Pi],scaling=constrained);([op
(r),u=1/2*Pi.3/2*Pi],scaling=constrained);
> cx: =r [1] - (diff (r [2],u)) * (
(diff (r [1],u)) ^2+ (diff (r [2],u)) ^2) / ( (diff (r [1],u)) * (diff (r
[2],u$2)) - (diff (r [2],u)) * (diff (r [1],u$2)));: =r [2] + (diff (r [1],u))
* ( (diff (r [1],u)) ^2+ (diff (r [2],u)) ^2) / ( (diff (r [1],u)) * (diff (r
[2],u$2)) - (diff (r [2],u)) * (diff (r [1],u$2)));
> l: = [cx,cy];
>pl_b: =plot ([op (l),u=-4.4],scaling=constrained,color=blue);
> display ([pl_a,pl_b]);
> restart:(plots):(plottools):
> a: =6;: =20;
> r: = [a*cosh (u),b*sinh (u)];[2];
>pl_a: =plot ([op
(r),u=-2*Pi.2*Pi],scaling=constrained);([op
(r),u=-2*Pi.2*Pi],scaling=constrained);
> cx: =r [1] - (diff (r [2],u)) * (
(diff (r [1],u)) ^2+ (diff (r [2],u)) ^2) / ( (diff (r [1],u)) * (diff (r
[2],u$2)) - (diff (r [2],u)) * (diff (r [1],u$2)));: =r [2] + (diff (r [1],u))
* ( (diff (r [1],u)) ^2+ (diff (r [2],u)) ^2) / ( (diff (r [1],u)) * (diff (r
[2],u$2)) - (diff (r [2],u)) * (diff (r [1],u$2)));
> l: = [cx,cy];
>pl_b: =plot ([op
(l),u=-2.2],scaling=constrained,color=blue);
> display ([pl_a,pl_b]);