МОДЕЛИРОВАНИЕ

Компьютерные модели школьных физических задач

Р.В. Бирих, Е.А. Еремин, В.И. Чернатынский,
г. Пермь

В этом году в серии “Жаркое лето” будет опубликовано наше пособие по курсу компьютерного моделирования. В качестве объектной области выбрана физика, предоставляющая множество проблем, которые могут быть эффективно решены методами математического моделирования. Предполагается не просто представить читателю набор отдельных задач, а изложить научный подход к моделированию физических процессов и продемонстрировать его приложение к школьной физике.

Публикация будет состоять из нескольких связанных частей: общая стратегия физического моделирования, особенности технологии математического моделирования на компьютере и, наконец, реализация этих принципов в школьных физических задачах. Будут рассмотрены вопросы постановки математической задачи на базе фундаментальных физических законов, выбора масштабов физических величин, анализа влияния погрешностей вычислений на достоверность физических результатов, особенности представления результатов на экране дисплея в графической и табличной формах и другие. Сформулированные общие положения будут проиллюстрированы на конкретных физических примерах. Далее описан один из них.

Постановка задачи. Задача о движении тела в центральном поле тяготения является хорошим примером, демонстрирующим возможности использования ПК для изучения поведения объекта, подчиняющегося некоторым общим физическим законам.

Процесс выведения спутника на орбиту обычно разбивается на два этапа. На первом этапе спутник поднимается над атмосферой практически вертикально на некоторую высоту. Затем обычно последняя ступень ракетоносителя придает спутнику необходимую горизонтальную скорость, и далее он двигается по инерции.

Рассмотрим инерционный полет небольшого тела (спутника) около притягивающего центра с большой массой (Земли). Будем интересоваться тем, какие траектории спутника возможны, какой должна быть его минимальная скорость вблизи поверхности Земли, чтобы он, двигаясь по круговой траектории, не упал на Землю (первая космическая скорость), какой должна быть минимальная начальная скорость спутника, чтобы получилась незамкнутая траектория и спутник ушел от Земли (вторая космическая скорость). В численном эксперименте можно также проверить законы Кеплера: 1) притягивающий центр расположен в одном из фокусов орбиты; 2) радиус-вектор от Солнца до планеты “заметает” равные площади за равные промежутки времени; 3) квадраты периодов обращения спутников вокруг притягивающего центра относятся как кубы больших полуосей эллипсов их орбит.

В основу модели мы положим закон всемирного тяготения. Будем считать, что на тело, движение которого рассматривается, действует только сила тяготения и уравнение Ньютона имеет вид:

.                         (1)

Здесь m и M — масса спутника и масса притягивающего центра, G — гравитационная постоянная, r — радиус-вектор, задающий положение спутника относительно притягивающего центра, a — ускорение спутника. Заметим, что хотя закон всемирного тяготения в форме (1) записан для материальных точек, он имеет такую же форму для сферически-симметричных тел.

Рис. 1. Система координат и возможная траектория спутника

Движение тела под влиянием центральной силы происходит в одной плоскости, положение которой определяется векторами r0 и v0, задающими начальное положение тела и его начальную скорость. Декартову систему координат с началом в центре тяготения и начало отсчета времени выберем так, чтобы движение происходило в плоскости Oxy и в начальный момент скорость тела была перпендикулярна оси x (рис. 1).

Тогда начальные условия можно записать в виде:

t = 0: x = x₀, y = 0, v_x = 0, v_y = v₀ (2)

Уравнение (1) вместе с условиями (2) полностью определяют траекторию спутника и все ее свойства.

Численная модель. Численный анализ задачи удобно проводить, используя в качестве единиц измерения характерные масштабы задачи. В качестве единицы длины удобно взять x0. Если разговор идет о спутнике Земли, то эта величина имеет порядок радиуса Земли R и равняется R + h, где h — высота спутника над поверхностью Земли. Всякое расстояние теперь будет задаваться числом, которое показывает, сколько раз в нем укладывается x0. Безразмерное х будет равняться х, измеренному в метрах, деленному на х0, также измеренному в метрах. Единицу времени удобно построить, используя гравитационную постоянную и характеристики притягивающего центра. Из уравнения (1) легко видеть, что множитель GM/r² имеет размерность ускорения (м/с²). Вместо расстояния r возьмем x₀ и сформируем выражение с размерностью времени (с): (GM/х₀³)^-1/2. Его и выберем в качестве единицы времени. В качестве единицы скорости тогда естественно взять х₀/(GM/х₀³)^-1/2, т.е. (GM/х₀)^1/2.

Измеренные в этих единицах проекции ускорения определяются следующими уравнениями (здесь и далее для безразмерных физических величин использованы те же обозначения, какие использовались для соответствующих размерных величин):

                            (3)

а начальные условия принимают вид:

t = 0: x = 1, y = 0, v_x = 0, v_y = V₀, (4)

где V₀ = v₀·(х₀GM)^1/2.

Все физические величины измеряются теперь в относительных единицах и будут одинаковыми для всех систем “спутник — притягивающий центр”. Уменьшилось также число параметров задачи. Единственный безразмерный параметр V₀, который остался в задаче, показывает, как соотносятся между собой кинетическая и потенциальная энергии спутника в начальный момент. Действительно, кинетическая энергия спутника в начальный момент равна  K = mv₀²/2, потенциальная энергия для ньютоновского закона тяготения равна П = GMm/x₀, и V₀² = 2К/П.

Для нахождения в различные моменты времени проекций скорости спутника и его координат на временной оси выберем дискретные точки t_n, отстоящие друг от друга на малые интервалы Dt. Тогда проекции скорости и в момент времени tn+1 будут приближенно (считаем, что ускорение на этом интервале времени не изменилось) представляться выражениями

(5) (6)

а координаты в этот момент будем вычислять, как при равномерном движении (опять считая, что интервал времени t мал, и скорость в течение него такая, как в конце интервала):

(7), (8)

В начальный момент времени проекции скорости и координаты спутника известны:

,

Система (5)–(8) позволяет шаг за шагом, при малом t, достаточно точно вычислить траекторию спутника и все ее характеристики.

Для классификации траекторий удобно вычислять их эксцентриситет е. В выбранной системе координат для вычисления эксцентриситета необходимо определить координату y = p точки на траектории, для которой . Тогда . Если е < 1, то траекторией является замкнутая кривая (окружность при е = 0 и эллипс при е 0 ), при e = 1 (в численном эксперименте это значение точно получено быть не может) траекторией является парабола, а при e > 1 траектория есть гипербола.

Площадь, “заметаемую” радиус-вектором за некоторое время , можно вычислить следующим образом:

(9)

Интересно сравнить эти площади, вычисленные на разных участках траектории. По второму закону Кеплера они должны быть равны.

По третьему закону Кеплера квадраты периодов обращения относятся как кубы больших полуосей орбит, т.е. для каждой траектории величина
К3 = Т ²/а³ должна оставаться постоянной. Для нахождения длины большой полуоси эллипса необходимо найти отрицательную координату x_мин при y = 0. Момент времени t_м, когда это произойдет первый раз, соответствует половине периода обращения спутника. Таким образом, величина большой полуоси равна: a = (1 + |x_мин|)/2, а период T = 2t_м.

Моделирование задачи в среде Delphi. Современная визуальная система программирования Delphi предоставляет удобное средство для создания среды моделирования любой физической задачи. Ее визуальная часть, будучи необычайно простой и естественной, позволяет быстро и без особых усилий создать на экране компьютера современный диалоговый интерфейс, который не уступает применяемому в профессиональном программном обеспечении. Лежащий в основе системы алгоритмический язык Паскаль является классическим учебным языком, так что запрограммировать на нем небольшую вычислительную задачу не представляет особой трудности. Наконец, встроенные в систему стандартные компоненты для поддержки таблиц и графики служат приятным дополнением к перечисленным выше достоинствам.

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

Перейдем к описанию среды моделирования нашей задачи. Как и любое Windows-приложение, она будет состоять из окон. Для наших целей естественным образом напрашивается создание четырех окон:

• главного с основными органами управления процессом численного моделирования,
• окна построения графиков
• и двух окон с таблицами наиболее важных физических характеристик: в первое помещаются интегральные характеристики процесса в целом, а во второе — промежуточные числовые результаты расчетов по каждому из графиков.

Возможный вид реализации виртуальной лабораторной среды приведен на рис. 2.

Рис. 2

Поскольку подавляющее большинство органов управления сосредоточено в главном окне, имеет смысл обратить основное внимание на его рассмотрение. Из приведенного рисунка устройство главного окна видно вполне отчетливо.

В левой части рассматриваемого окна располагается панель с 4 полями ввода, используя которые можно регулировать отображаемую в графическом окне область. Рядом со схематическим прямоугольником графического экрана предусмотрены поля для ввода минимального и максимального значения по каждой из осей X и Y. Отметим, что координаты, как и все остальные вводимые в главном окне параметры, указываются непосредственно в “физических” единицах; необходимый для рисования перевод в пиксели экрана берет на себя программа.

Следующая панель, занимающая основную часть главного окна в его центре, служит для задания параметров задачи — скорости запуска спутника V0 и шага по времени dt (t — в постановке задачи) — а также запуска и прерывания расчетов. Предусмотрено задание значения V₀ двумя способами. В первом случае строится график для одного значения параметра: оно задается в поле правее слов “Один график для” и программа запускается на счет кнопкой “Счет-1”. Во втором производится автоматическое вычисление сразу нескольких траекторий: вводятся начальное и конечное значения V₀, для которых необходимо произвести расчет, из списка выбирается количество требуемых траекторий (их может быть от 2 до 11) и нажимается кнопка “Счет-М”. Сочетание двух указанных режимов проведения расчетов позволяет производить численные эксперименты значительно удобнее, получая любые требуемые траектории. В обоих случаях можно в любой момент прервать расчет нажатием кнопки “Стоп”.

Примечание. В Delphi при вводе дробных чисел в компоненты полей ввода целую и дробную части следует разделять символом, который определяется национальными настройками ОС Windows, а именно: в русской версии необходимо использовать запятую, а в английской — точку. Для языка программирования подобное решение проблемы кажется несколько непоследовательным (согласно их синтаксису в вещественных числах всегда однозначно используется точка); тем не менее, если вспомнить электронные таблицы Excel, то там тоже имеет место аналогичный “двойственный” подход.

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

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

Окна с графиками функций и табл. 1 интегральных результатов для каждой траектории никаких органов управления не содержат. Зато для построения табл. 2 в правой части окна помещены два поля для ввода параметров построения данной таблицы. В первое вводится количество столбцов табл. 2, которое с точки зрения физического содержания задачи равняется числу точек траектории, для которых будет зафиксировано значение некоторой интересующей нас характеристики (в данной задаче это площадь сектора). Что касается второго поля, то в него задается значение времени, через которое эти значения будут заноситься в таблицу. Например, при указанных на рисунке значениях в таблицу для каждой траектории попадет по 10 точек для значений времени 0,1, 0,2, …, 1.

Конкретные детали реализации органов управления, которые определяются программами, обрабатывающими их нажатия, будут приведены в летнем выпуске.

Эксперимент на модели в среде Delphi. Проведем экспериментальное исследование движения спутника с помощью описанной выше модели. Для проведения численного эксперимента нужно прежде всего выбрать шаг по времени dt, с которым будем строить траекторию. По смыслу метода он должен быть достаточно малым, чтобы на этом интервале можно было пренебречь изменением ускорения при вычислении скорости и изменением скорости при вычислении новых координат. Шаг не должен быть и слишком маленьким, так как в этом случае для построения траектории потребуется большое время. Поскольку в качестве масштабов мы выбрали характерные параметры системы, то типичный временной интервал жизни системы равен 1. Малый интервал dt поэтому должен быть много меньше единицы; выберем его, например, равным 0,01. Второй параметр, который нужно задать, это начальная скорость. Опять, пользуясь тем, что расчет ведется в собственных масштабах задачи, для первого эксперимента естественно положить V₀ = 1. Теперь можно построить одну траекторию, нажав кнопку “Счет-1”. Посмотрим на полученную траекторию. Она похожа на окружность, но не замкнулась при обороте спутника вокруг притягивающего центра (внутренняя кривая на графиках рис. 3). Это, по-видимому, связано с недостаточной точностью счета. Давайте уменьшим величину dt в два раза: dt = 0,005 и изменим цвет кривой. Экран и таблицы чистить не будем, чтобы сравнить результаты, полученные для разных dt. Нажимаем кнопку “Счет-1”. Результат лучше, но требуется еще более мелкий шаг. Сделаем dt = 0,001.

Рис. 3. Состояние окон с результатами после эксперимента по выбору dt

Траектория получилась замкнутая. Если мы теперь посмотрим на интегральные характеристики построенных кривых (табл. 1 на рис. 3), то заметим, что эксцентриситет е с увеличением точности счета уменьшался и приблизился к нулю. Может быть не ноль, потому что еще грубо считаем? Уменьшим dt до 0,0001. Новый цвет кривой и кнопка “Счет-1”. (На рис. 3 приведено состояние окон после проведенных расчетов.) Вот теперь хорошо! Эксцентриситет равен нулю, большая полуось равна единице — значит, траектория окружность. Первая космическая скорость в принятых единицах равна 1, а период обращения —. В нижней таблице видно, как хорошо соблюдается для этой траектории второй закон Кеплера: площадь, “заметаемая” радиус-вектором за одно и то же время на разных участках траектории, одинакова.

Теперь можно посмотреть, какие еще бывают траектории спутника. Построим серию траекторий из 11 линий с диапазоном начальных скоростей от 1 до 2. Очистим графическое поле и таблицу, наберем необходимые параметры и кнопкой “Счет-М” запускаем программу. На рис. 4 приведены результаты этого эксперимента. Видно, что с увеличением начальной скорости эллипсы вытягиваются вдоль оси х. Притягивающий центр (начало координат) находится в правом фокусе эллипсов. Последняя замкнутая траектория получилась для V₀ = 1,4 (она так сильно вытянута, что не поместилась на экране). Для этой траектории эксцентриситет близок 1, а значит, соответствующая ей начальная скорость близка ко второй космической скорости. Анализ данных таблиц показывает, что второй и третий законы Кеплера хорошо выполняются.

Рис. 4. Состояние окон с результатами после расчета серии траекторий

Моделирование задачи в среде MS Excel. Рассмотрим теперь реализацию модели с помощью MS Excel. Общий алгоритм состоит из трех блоков. Первый — формирование таблицы значений координат спутника в последовательные моменты времени; второй — построение графика с помощью встроенного мастера диаграмм; третий — оформление интерфейса.

Первый блок содержит следующие шаги:

а) определение параметров таблицы,

б) определение структуры записи, т.е. данных, которые должны быть расположены в строке таблицы,

в) определение количества записей,

г) заполнение таблицы расчетными формулами и необходимыми значениями.

Шаги, образующие второй блок, однозначно определяются стандартным мастером диаграмм Excel.

Третий блок определяет удобное для проведения исследований расположение параметров, таблиц значений и графиков зависимостей.

Рассмотрим перечисленные шаги более подробно.

Параметры таблицы. В качестве параметров следует выбирать величины, задание которых определяет размеры таблицы и числовые значения данных, получаемых в ее ячейках. Предполагается, что значения параметров не изменяются в процессе вычислений; они размещаются в фиксированных ячейках, при обращении к которым используется абсолютная адресация. В сформулированной задаче значения координат спутника в любой момент времени однозначно определяются начальными условиями, поэтому их естественно определить в качестве параметров. Выбранный метод пошагового решения задачи содержит в качестве основной характеристики шаг по времени dt, который также удобно включить в список параметров. Таким образом, таблицу определяют четыре физических параметра — x₀, y₀, vx₀, vy₀ и один вычислительный — dt.

Структура записи. В ячейках строки таблицы (слева направо) следует расположить такие данные: время t, проекция скорости на ось x, проекция скорости на ось y, координата x, координата y, мгновенная секториальная скорость и две вспомогательные ячейки. Все данные (кроме секториальной скорости) вычисляются по формулам численного метода, которые представляют собой рекуррентные соотношения. С их помощью значение в данной ячейке определяется значениями величин в предыдущий момент времени, т.е. в предыдущей строке. Вспомогательные ячейки используются для выделения значений координат x и y, которые определяют такие параметры траектории, как эксцентриситет е и большая полуось а.

Количество записей. Размер таблицы, очевидно, определяется количеством шагов по времени, которое, в свою очередь, зависит от величины шага по времени dt. Поэтому мы должны сделать выбор до решения задачи. Практически шаг по времени определяется устойчивостью вычислительной процедуры и точностью решения, а максимальное число шагов задается в процессе решения каким-либо дополнительным условием (условием замыкания траектории, ограничением значений координат для незамкнутых траекторий и т.д.). Обойти эти трудности можно двумя путями. Первый состоит в изменении количества записей в зависимости от траектории, сохраняя постоянный шаг по времени. Другой способ заключается в создании большой таблицы для маленького шага по времени, используя для определения количества записей круговую траекторию в качестве эталона; обследовать траектории, соответствующие разным периодам обращения спутника, можно, изменяя шаг по времени.

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

Построение траекторий. Для построения траектории рекомендуется использовать диаграмму точечного типа, сглаженную, без координатных точек. Чтобы получить правильные пропорции кривой, следует определить оси системы координат одинакового размера, выбирая начало координат по горизонтали с учетом того, что максимальное значение этой координаты равно 1. Для увеличения полезной площади диаграммы легенду рекомендуется удалить.

Организация интерфейса. Для более удобной работы предлагается начальные данные с графиком траектории и таблицу результатов разместить на отдельном листе Книги MS Excel, а громоздкую таблицу с расчетами значений координат точек траектории — на следующем.

Следует обратить внимание на то, что математическая модель, создаваемая в MS Excel, оказывается менее гибкой и наглядной, чем модель, полученная с помощью среды программирования Delphi. В частности, в последней можно создать иллюзию движения спутника по траектории, представить для сравнения множество траекторий в одних осях, удобно провести вычисления характеристик траекторий для проверки законов Кеплера.

Эксперимент на модели в среде MS Excel. Расположим интерфейсную часть модели, т.е. параметры, диаграмму и характеристики траекторий, на первом листе книги, а таблицу вычисленных значений переменных на втором листе. Пусть параметры x₀, y₀, vx₀, vy₀ и dt в перечисленном порядке располагаются в ячейках B2, C2, D2, E2 и F2 Листа 1, соответственно. Тогда первые две строки таблицы значений на Листе 2, начиная со строки 4, заполняются формулами, приведенными в табл. 1.

Последующие записи получаются копированием строки 5, предварительно вычислив требуемое количество строк.

Размер таблицы определим из следующих соображений. При движении по окружности радиуса, равного 1, с характерной скоростью, равной 1, полный оборот совершается за промежуток времени, равный 6,28. Шаг по времени выберем равным 0,001, так что таблица будет содержать 6280 строк. Для заполнения такой большой таблицы рекомендуется использовать при копировании строк поле адреса строки ввода.

Для вычисления характеристик траектории и проверки законов Кеплера таблица значений дополнена тремя вспомогательными столбцами — F, G, H. В столбце F, начиная с первой строки таблицы, вычисляется комбинация проекций скорости и координат спутника для нахождения средней секториальной скорости на заданном промежутке времени. В столбце G с помощью функции ветвления ЕСЛИ и логического произведения И выделяются значения координаты y > 0 (столбец E4) в малой окрестности точки x = 0. Остальные ячейки столбца обнуляются. Таким образом, создается структура значений, позволяющая с помощью интервальной функции MAX найти приближенно координату точки, в которой траектория пересекает ось y, и вычислить эксцентриситет е. Аналогичная структура данных создается во вспомогательном столбце H для вычисления большой полуоси. При этом, поскольку выделяется малый диапазон отрицательных значений x, необходимо использовать функцию MIN. Соответствующие формулы приведены в табл. 2, являющейся фрагментом строки 4. Функции MIN и MAX используются в столбцах а и e таблицы результатов.

Вычисления указанных характеристик производятся в таблице результатов, которая также включается в интерфейс. Период обращения определяется как значение ячейки A6280, когда получена замкнутая траектория. Среднее значение секториальной скорости, вычисленное путем усреднения формулы (9) для 20 точек в окрестности значений горизонтальной координаты, соответствующих максимальному и минимальному удалению от силового центра, помещается в раздел К2 таблицы результатов.

Таблица 2

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

Рис. 5. Результаты эксперимента на модели MS Excel

Примеры вопросов для самостоятельного решения

1. Какие в проведенных исследованиях получаются размерные первая и вторая космические скорости для спутника Земли?

Ожидаемый ответ. Нужно вычислить единицу скорости для спутника у поверхности Земли: (GM/x₀)^1/2= (6,7 · 10-11м³/(кг ·   с²) ·  6,0 · 10²⁴ кг/(6,4 · 106 м))^1/2 = 7.9 км/с. Для первой космической скорости это число нужно умножить на 1, а для второй — на 1,4.

2. Какие траектории спутника получатся, если начальная безразмерная скорость его меньше единицы?

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

3. Как изменятся значение площади, заметаемой радиус-вектором, и средняя секториальная скорость, если интервал времени, за который эти величины считаются, увеличить вдвое? В 5 раз?

Ожидаемый ответ. Площадь увеличится в соответствующее число раз. Убедиться в этом можно, анализируя содержимое табл. 2 в разных экспериментах. Средняя секториальная скорость не изменится.

Примеры заданий для самостоятельного исследования

1. Построить 11 траекторий для значений начальной скорости в диапазоне от 0,8 до 1,8. Используя режим расчета одиночной траектории, нанести на рисунок красным цветом круговую траекторию, а розовым — самую вытянутую траекторию, целиком уместившуюся на экране.

2. Произведите расчет для значений начальной скорости 0,4. Получилась ли траектория замкнутой? Как добиться, чтобы она стала замкнутой (уменьшить шаг по времени)?

3. Постройте 6 траекторий для значений начальной скорости в диапазоне от 1,4 до 1,42. Определите, сколько из них будут замкнутыми. Найдите в таблице значение большой полуоси для самой вытянутой замкнутой траектории. Попробуйте подобрать масштабы по осям x и y, чтобы эта траектория поместилась в окне графика.

4. Подобрать такую начальную скорость спутника, чтобы период его обращения T был вдвое больше, чем при обращении по круговой орбите.

5. Исследовать, как влияет величина шага по времени на значение константы третьего закона Кеплера для разных траекторий.

Примечание. Количество траекторий при моделировании средствами MS Excel может быть уменьшено.