Решение физических задач в электронных таблицах Excel: учебное пособие [Роберт Валерьевич Майер] (pdf) читать онлайн

Майер Р. В. Решение физических задач в электронных таблицах Excel
Об издании 1, 2

Содержание

1

Вперёд

Майер Р. В. Решение физических задач в электронных таблицах Excel

МИНИСТЕРСТВО ОБРАЗОВАНИЯ И НАУКИ
РОССИЙСКОЙ ФЕДЕРАЦИИ
ФЕДЕРАЛЬНОЕ ГОСУДАРСТВЕННОЕ БЮДЖЕТНОЕ
ОБРАЗОВАТЕЛЬНОЕ УЧРЕЖДЕНИЕ
ВЫСШЕГО ПРОФЕССИОНАЛЬНОГО ОБРАЗОВАНИЯ

«ГЛАЗОВСКИЙ ГОСУДАРСТВЕННЫЙ ПЕДАГОГИЧЕСКИЙ ИНСТИТУТ
имени В. Г. Короленко»

Р. В. Майер
РЕШЕНИЕ ФИЗИЧЕСКИХ ЗАДАЧ
В ЭЛЕКТРОННЫХ ТАБЛИЦАХ EXCEL
Учебное пособие
Учебное электронное издание
на компакт-диске

Глазов
ГГПИ
2016

© Майер Р. В., 2016
© ФГБОУ ВПО «Глазовский государственный педагогический
институт имени В. Г. Короленко», 2016

ISBN 978-5-93008-211-1

2

Майер Р. В. Решение физических задач в электронных таблицах Excel
Назад

Содержание

Вперёд

1 – дополнительный титульный экран – сведения об издании

УДК 37.02:004.02
ББК 32.973.26
M14
Рекомендовано УМО по математике педвузов и университетов
Волго-Вятского региона в качестве учебного пособия
для студентов и преподавателей высших учебных заведений

Автор: Майер Роберт Валерьевич
Рецензенты:
В.А. Саранин, доктор физико-математических наук, профессор кафедры физики
и дидактики физики Глазовского государственного педагогического института;
Ю.А. Сауров, доктор педагогических наук, профессор кафедры физики и методики обучения физики Вятского государственного гуманитарного университета, членкорреспондент Российской Академии Образования
М14 Майер Р. В.
Решение физических задач в электронных таблицах Excel: учебное
пособие [Электронное учебное издание на компакт-диске]. – Глазов: Глазов.
гос. пед. ин-т, 2016. – 14,0 Мб.
Электронное учебное пособие посвящено использованию электронных таблиц
Excel для решения физических задач на занятиях по физике и компьютерному моделированию. В нем представлены физические задачи по кинематике, динамике
системы частиц, механическим и электрическим колебаниям, молекулярной физике, теплопроводности, электродинамике, волновому движению, оптике, динамическому хаосу, рассмотрены методы их решения в таблицах Excel. Пособие содержит
большое количество примеров, которые позволяют учащимся и студентам самостоятельно научиться создавать математические и компьютерные модели, выполнять вычислительные эксперименты.
Электронное учебное пособие предназначено преподавателям и студентам,
изучающим математику, физику и информатику, а также всем тем, кто интересуется
использованием электронных таблиц Excel для решения различных задач.
Системные требования: РС не ниже класса Pentium I; 128 Mb RAM; свободное
место на HDD 16 Mb; Windows 95/98/2000/XP/7/8; Adobe Acrobat Reader; дисковод
CD-ROM 2-х и выше; мышь.
© Майер Р. В., 2016
© ФГБОУ ВПО «Глазовский государственный педагогический институт
имени В. Г. Короленко», 2016

3

Майер Р. В. Решение физических задач в электронных таблицах Excel
Назад

Содержание

Вперёд

2 – дополнительный титульный экран – производственно-технические сведения

УЧЕБНОЕ ИЗДАНИЕ
Электронное учебное пособие Р. В. Майера «Решение физических задач в
электронных таблицах Excel» поставляется на одном CD-ROM и может быть использовано в локальном и сетевом режимах. В случае, когда система установлена на
одном из серверов вычислительной сети, к ней обеспечивается одновременный доступ нескольких пользователей.

Технический редактор, корректор В. В. Баженова
Оригинал-макет: И. С. Леус

Подписано к использованию 30.05.2016. Объём издания 24,3 Мб.
Тираж 12 экз. Заказ № 1754 – 2016.
ФГБОУ ВПО «Глазовский государственный педагогический институт
имени В. Г. Короленко»
427621, Россия, Удмуртская Республика, г. Глазов, ул. Первомайская, д. 25
Тел./факс: 8 (34141) 5-60-09, e-mail: izdat@mail.ru

4

Майер Р. В. Решение физических задач в электронных таблицах Excel

Назад

Содержание

1. Кинематика точки:
одномерное движение

Введение
Развитие цифровой педагогики предполагает разработку простых
компьютерных моделей физических, технических, биологических и социальных систем, которые были бы понятны студентам и школьникам.
Метод компьютерного моделирования [4–20] целесообразно применять
тогда, когда другие способы решения рассматриваемой задачи менее
результативны, либо слишком трудоемки и требуют сложных математических выкладок. Существующая методика компьютерного моделирования позволяет решить большое количество различных физических задач. Важно при этом использовать доступные программные средства,
освоение которых не вызывает сложностей. Школьный и вузовский курсы информатики предусматривают изучение табличного процессора MS
Excel. Он является мощным программным средством, которое объединяет в себе электронные таблицы, средства визуального программирования и графический модуль, позволяющий построить различные диаграммы, графики и поверхности. Поэтому при решении физических задач имеет смысл использовать именно этот программный продукт.
Хотя пакет MS Excel имеет меньше возможностей по сравнению со
специализированными пакетами (MathCad, MatLab, Math и т. д.), он позволяет реализовать простейшие алгоритмы численного решения диффуравнений и создать компьютерные модели. Известные книги и учебные пособия [2; 3; 21] не дают полного представления о возможностях
использования электронных таблиц при изучении физики. Поэтому проблема использования табличного процессора Excel для решения
физических задач остается актуальной. Так как с помощью пакета Excel можно промоделировать большое количество разнообразных физических систем, требующих численного решения дифференциальных
уравнений и нахождения определенных интегралов, то его возможно и
целесообразно использовать при изучении физики и основ компьютерного моделирования.
Методы решения задач в электронных таблицах
Для решения задач в табличном процессоре Excel будем использовать два способа. Первый способ предполагает табулирование функ-

5

Майер Р. В. Решение физических задач в электронных таблицах Excel

ций и построение графиков без использования макросов. Например, для
расчета движения камня, брошенного под углом к горизонту, в столбец A
вводят значения времени 0, 0,1, 0,2, 0,3 и т. д., в столбцы B и C вводят
формулы для расчета координат материальной точки x = x0 + υ0 cos α ⋅ τ ,

y = y0 + υ0 sin α ⋅ τ − aτ 2 / 2 , закодированные соответствующим образом, и
копируют их в ячейки, расположенные ниже. Используется относительная адресация, формулы вводятся в виде “=6+8*SIN(0,87)*A3–
–9,8*A3*A3/2”. При их копировании изменяются адреса ячеек. Параметры системы и начальные условия записываются в специальные ячейки
электронной таблицы, при этом в использующих их формулах применяется абсолютная адресация: $E$2, $F$3 и т. д. На основе получающейся таблицы стандартными средствами Excel строят графики исследуемых функций.

Рис. 1. Создание макроса для моделирования вынужденных колебаний
Для решения более сложных задач следует применять второй способ, предусматривающий написание макроса – небольшой программы,
создаваемой пользователем в специальном блоке. Макрос пишется на
языке Visual Basic и может содержать цикл по времени, координате, использоваться для нахождения интеграла, производной, многократном
пересчете элементов массива и т. д. В некоторых задачах с помощью
программ VBA (Visual Basic Application) рассчитываются координаты

6

Майер Р. В. Решение физических задач в электронных таблицах Excel

движущихся частиц, решается соответствующее уравнение в частных
производных, строится сечение Пуанкаре, проводится серия испытаний
и т. д. Для написания макроса достаточно выбрать: Вид → Панели инструментов → Элементы управления → Кнопка [2; 3]. Необходимо кнопку
Command Button1 перенести на таблицу и дважды кликнуть по ней. В
появившееся окно следует ввести или скопировать текст программы,
которая будет исполняться после запуска. Макрос считывает данные из
указанных в нем ячеек электронной таблицы и, произведя расчеты, создает таблицу результатов вычислений. На основе этой таблицы строится график изучаемой зависимости. Макрос может быть составлен так,
что при повторном нажатии на кнопку программа увеличивает время на
заданное время ∆τ , повторяет расчеты и строит новый график. В некоторых задачах физические величины измеряются в условных единицах.
На рис. 1 рассмотрен пример макроса, моделирующего вынужденные колебания физического маятника. После запуска программы получаются таблицы результатов решения соответствующего дифференциального уравнения второго порядка: Iϕ&& + rϕ& + mgl sin ϕ = M max sin ω1τ ,

где τ – время (рис. 2). В столбцы A, B, C, D выводятся время, угол отклонения, угловая скорость маятника и внешний вращающий момент,
вызывающий его колебания. На основе полученной таблицы строят
графики ϕ (τ ) , ω (τ ) и M (τ ) . При этом используют мастер диаграмм, вы-

бирают опцию “точечная диаграмма со значениями, соединенными
сглаживающими линиями”.

Рис. 2. Результаты моделирования вынужденных колебаний

7

Майер Р. В. Решение физических задач в электронных таблицах Excel

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

8

Майер Р. В. Решение физических задач в электронных таблицах Excel

Введение

2. Кинематика точки:
двумерное движение

Содержание

1. КИНЕМАТИКА ТОЧКИ: ОДНОМЕРНОЕ ДВИЖЕНИЕ
1.1. Точка движется вдоль оси Ox , ее скорость изменяется по закону: υ (t ) = 0,1 + 0,3τ + 0,2τ 2 − 0,1τ 3 (м/c). Постройте графики зависимости

x(τ ) и υ (τ ) . Найдите координату и скорость точки в момент τ ' = 3,2 c, а
также среднюю скорость за первую и третью секунды движения.
Запишем формулы: ∆x = υ (τ )∆τ , x t +1 = x t + υ t ∆τ . Координата x
находится как интеграл от скорости по времени, для ее вычисления
применяется метод прямоугольников. Шаг по времени ∆τ = 0,1 помещают в ячейку E1. В ячейку A1 таблицы вводят время τ = 0, в ячейку
A2 – ‘=A1+$E$1’, затем эту формулу копируют вниз. Высота таблицы
50 строк. В ячейку B1 вводят формулу для расчета скорости точки
‘=0,1+0,3*A1+0,2*A1*A1–0,1*A1*A1*A1’, а затем ее копируют вниз. В
ячейку С1 записывают формулу ‘=B1*$E$1’, а в ячейку C2 помещают
формулу ‘=C1+B2*$E$1’, которую также копируют вниз. В ячейках
столбца C вычисляется приращение координаты за время ∆τ и складывается с координатой точки в предыдущий момент времени. На
основе таблицы строят графики x(τ ) и υ (τ ) .
1.2. Поезд выехал из А в В со скоростью 70 км/ч. Через 3 часа после начала движения из B вышел второй поезд, движущийся со скоростью 80 км/ч по направлению к A. Расстояние от A до B 1560 км. Где и
когда поезда встретятся? Постройте графики движения x1 (τ ) и x2 (τ ) .
1.3. Две частицы движутся вдоль оси Ox в соответствии с уравнениями: x1 = −5 + 13τ (м) и x2 = 235 + 4τ − 2τ 2 (м). Постройте графики зависимостей координат и скоростей частиц от времени, определите время и место их встречи.
1.4. Точка в течение 2 с двигалась равномерно со скоростью 6 м/с,
2

затем 3,2 с двигалась с ускорением 1,2 м/с . Следующие 4,3 с она двигалась равномерно, а затем 3 – равнозамедленно с ускорением –1,6
2

м/с . Постройте графики зависимостей координаты и скорости точки от
времени, определите пройденный путь, вычислите среднюю скорость.
1.5. Точка 5 с движется по закону υ = 1,2 ⋅ τ (м/c), затем 10 с равномерно со скоростью 6 м/c, после чего скорость уменьшается по закону
9

Майер Р. В. Решение физических задач в электронных таблицах Excel

υ (τ ) = 6 − 0,7(τ − 15) (м/c). Постройте графики зависимости координаты и
скорости точки от времени. Определите скорость и координату точки в
моменты 9,3 с и 17,6 с. Вычислите среднюю скорость за первые 14 с.
1.6. Точка движется вдоль оси Ox c начальной скоростью 3 м/c.
Ускорение (в м/с 2 ) задано формулой:

если 0 ≤ τ < 5,
0,2τ ,


a x (τ ) =  1 + 0,3(τ − 5), если 5 ≤ τ < 8,
3,4 − 0,4(τ − 8), если 8 ≤ τ < 12;

Постройте графики зависимостей υ x (τ ) и x(τ ) , если начальная координата x0 = 1 м. Определите координату и скорость точки в момент

τ ' = 9,4 с, а также среднюю скорость за вторую секунду движения.
1.7. В момент τ = 0 камень падает без начальной скорости с высоты 23 м. В момент τ = 0,3 c второй камень бросают вверх со скоростью
8 м/c. Постройте графики зависимостей координат и скоростей тел от
времени. Определите момент времени, когда камни окажутся на одной
высоте h . Найдите высоту h .
1.8. Дробинка падает в вязком масле так, что ее скорость изменяется по закону υ (τ ) = 14(1 − e − 0,07⋅τ ) (м/с). Постройте графики зависимостей координаты, скорости и ускорения от времени. Определите момент
времени, когда скорость дробинки будет равна 8 м/с, а также среднюю
скорость за первую и третью секунды движения.
1.9. Точка колеблется по закону x(τ ) = 15 sin(4τ + 0,1) (см). Постройте графики зависимостей координаты, скорости и ускорения точки с течением времени, а также фазовую кривую в осях x и

υ x . Определите

координату, скорость и ускорение точки в момент τ ' = 3,7 c.
1.10. Точка совершает вынужденные колебания, описываемые
уравнением x(τ ) = 13(1 − e −0,05⋅τ ) sin(8τ + 0,8) (см). Постройте графики
зависимостей координаты, скорости и ускорения точки с течением времени, а также фазовую кривую в осях x и υ x . Определите координату,
скорость и ускорение точки в момент τ ' = 4,2 c.
1.11. Тело совершает затухающие колебания, описывающиеся
уравнением: x(t ) = Ae − β ⋅τ sin(ω ⋅ τ + ϕ 0 ) . Постройте графики зависимостей координаты x и проекции скорости υ x от времени, а также фазовую кривую при различных значениях амплитуды A , частоты ω , началь10

Майер Р. В. Решение физических задач в электронных таблицах Excel

ной фазы

ϕ 0 и коэффициента затухания β . Например, A = 10 см, ω = 5

рад/с, ϕ 0 = 1,3 рад, β = 0,1 с-1.
1.12. В результате сложения двух колебаний в одном направлении
точка движется по закону x(τ ) = 12 sin(3,3τ + 0,6) + 8 sin(3,6τ − 0,2) (см).
Постройте графики зависимости координаты, скорости и ускорения точки от времени. Определите значения x , υ x и ax в момент τ ' = 3,4 с.
1.13. Точка одновременно участвует в двух колебаниях, происходящих в одном направлении. Их уравнения: x1 (τ ) = 12 sin(ω1 ⋅ τ − 0,3) и

x2 (τ ) = 9 sin(ω2 ⋅ τ + 1,7) . Получите графики результирующих колебаний в
следующих случаях: 1) ω1 = 4,3 рад/c и ω2 = 4,4 рад/c; 2) ω1 = 4,3 рад/c и

ω2 = 4,8 рад/c; 3) ω1 = 4,3 рад/c и ω2 = 5,3 рад/c. Для каждого случая определите период биений.
1.14. Точка колеблется по закону x(τ ) = 12 sin(5τ − 1,3) (см). Определите первые три момента времени, при которых x = −3,2 см. Найдите
скорость и ускорение точки в эти моменты времени.
1.15. Точка колеблется по закону x(τ ) = 23e − 0,05τ sin(2,4τ − 0,3) см.
Определите первые три момента времени, при которых x = 8 см. Найдите скорость и ускорение точки в эти моменты времени.
1.16. Точка колеблется по закону x(τ ) = 15e − 0,04τ sin(3,2τ + 0,4) (см).
Определите первые три момента времени, при которых υ x = −16 см/с.
Найдите координату и ускорение точки в эти моменты времени.

11

Майер Р. В. Решение физических задач в электронных таблицах Excel
1. Кинематика точки:
одномерное движение

3. Динамика:
одномерное движение

Содержание

2. КИНЕМАТИКА ТОЧКИ: ДВУМЕРНОЕ ДВИЖЕНИЕ
2.1. Камень брошен со скоростью υ 0 = 17 м/c под углом α к горизонту. Рассчитайте его координаты x и y , проекции скорости υ x и υ y ,
модуль скорости υ , угол γ между вектором скорости и горизонталью,
тангенциальное aτ и нормальное an ускорения в момент τ '= 1,3 с.
Запишем формулы, позволяющие определить искомые величины:

x = υ0 cos α ⋅τ ,

y = υ0 sin α ⋅τ − gτ 2 / 2 ,

υ x = υ0 cos α , υ y = υ0 sin α − gτ , υ = υ x2 + υ y2 ,
γ = arctg (υ y / υ x ) ,

an = g cos γ ,

aτ = g sin γ .

Для того чтобы можно было изменять исходные данные, задаваемые значения начальной скорости υ0 и угла α следует ввести в
отдельные ячейки, например M2 и M3. Чтобы перевести градусы в
радианы, в ячейку N3 вводят формулу: “=M3*3,1415/180”. Для решения
задачи в Excel создается таблица; содержимое ее ячеек представлено
в табл. 2.1. Получающиеся результаты приведены на рис. 2.1.
Таблица 2.1
Ячейки

Величина

Содержимое

B1, B2, B3, …

τ

0; 0,2; 0,4; 0,6; 0,8; 1; 1,2 ...

C3

=$M$2*COS($N$3)*B3

D3

x
y

=$M$2*SIN($N$3)*B3-9,8*B3*B3/2

E3

υx

=$M$2*COS($N$3)

F3

υy

=$M$2*SIN($N$3)-9,8*B3

H3

υ
γ

=ATAN(F3/E3)*180/3,1415

I3

an

=9,8*COS(ATAN(F3/E3))

J3

aτ

=9,8*SIN(ATAN(F3/E3))

G3

=КОРЕНЬ(E3*E3+F3*F3)

12

Майер Р. В. Решение физических задач в электронных таблицах Excel

Рис. 2.1. Расчет движения камня, брошенного под углом к горизонту
2.2. Камень брошен со скоростью

υ0 = 13 м/c под углом α = 54 0 к

горизонту. Рассчитайте его координаты x и y , проекции скорости υ x и

υ y , модуль скорости υ , угол γ между вектором скорости и горизонталью, тангенциальное aτ и нормальное an ускорения в моменты времени

τ i = t ⋅ ∆τ , ∆τ = 0,2 с, t = 1, 2, 3,...

2.3. Камень брошен горизонтально со скоростью

υ0 = 14 м/c с

башни высотой 23 м. Рассчитайте его координаты x и y , проекции скорости υ x и υ y , модуль скорости υ , угол γ между вектором скорости и
горизонталью, тангенциальное aτ и нормальное an ускорения в моменты времени τ i = t ⋅ ∆τ , ∆τ = 0,2 с, t = 1, 2, 3,...

2.4. Стержень ОА длиной 14 см колеблется относительно точки О
по закону ϕ (τ ) = 1,5 sin(2 ⋅τ ) рад. Постройте графики зависимости модуля
линейной скорости точки A, а также нормального, тангенциального и
полного ускорений. Определите значения этих величин в момент, когда:
1) маятник проходит положение равновесия; 2) находится в крайнем положении; 3) отклонен на угол 0,7 рад. Нарисуйте рисунок, на котором
покажите эти вектора в заданные моменты времени.
2.5. Колесо радиусом 17 см вращается с ускорением ε = 0,3 + 0,1τ
(рад/c2), начальная скорость равна 0. Постройте графики зависимости

13

Майер Р. В. Решение физических задач в электронных таблицах Excel

угловой скорости ω , угла поворота ϕ , тангенциального aτ и нормального an ускорений точки А на ободе колеса от времени. Определите мо-

r

мент времени, когда угол между вектором полного ускорения a и радиу0

сом колеса равен 38 .
2.6. Точка одновременно участвует в двух гармонических колебаниях, происходящих во взаимно перпендикулярных направлениях. Амплитуды колебаний Ax , Ay , их частоты ω x , ω y и сдвиг фаз ∆ϕ заданы.
Необходимо определить траекторию движения точки и ее координаты x ,
y в произвольный момент τ ' .
Использующаяся программа 2.1 содержит цикл по времени; в нем
вычисляются координаты x(τ ) = Ax sin(ω x ⋅τ ) , y (τ ) = Ay sin(ω y ⋅ τ + ϕ 0 ) , а
результаты записываются в ячейки таблицы Excel. Значения амплитуд, частот колебаний и сдвига фаз считываются из ячеек H1–H5.
Получающаяся таблица и графики представлены на рис. 2.2.
Программа 2.1
Private Sub CommandButton1_Click()
Ax = Cells(1, 8): Ay = Cells(2, 8)
wx = Cells(3, 8): wy = Cells(4, 8)
fi = Cells(5, 8): dt = 0.2
While t < 60
t = t + dt: i = i + 1
x = Ax * Cos(wx * t)
y = Ay * Cos(wy * t + fi)
Cells(i, 1) = x: Cells(i, 2) = y
Wend
End Sub

Программа позволяет промоделировать движение точки при
различных амплитудах, частотах и сдвиге фаз. В случае когда частоты кратны целым числам (то есть ω x / ω y = n / m , где n, m ∈ Z ),
точка движется по замкнутой кривой, называемой фигурой Лиссажу.

14

Майер Р. В. Решение физических задач в электронных таблицах Excel

Рис. 2.2. Движения точки, колеблющейся во взаимно
перпендикулярных направлениях при различных отношениях ω x / ω y
2.7. Точка одновременно участвует в двух гармонических колебаниях, происходящих во взаимно перпендикулярных направлениях. Уравнения
колебаний
имеют
вид
x(τ ) = 8 sin(3,2 ⋅ τ − 0,3)
(см)
и

y (τ ) = 6 sin(2,8 ⋅ τ + 1,1) (см). Найдите модуль скорости точки и угол, который составляет вектор скорости с осью Ox в момент времени 2,8 с.
2.8. Точка движется по закону: x(τ ) = 5sh(0,2τ ) (м), y (τ ) = 3τ (м).
Постройте траекторию движения, определите скорость точки в момент
τ ' = 12,3 с.
Из определения гиперболического синуса получаем:

x(τ ) = 5

exp(0,2τ ) − exp(−0,2τ )
,
2

y (τ ) = 3τ .

Скорость точки вычисляется по формулам:

υ x = ( x t +1 − x t ) / ∆τ , υ y = ( y t +1 − y t ) / ∆τ , υ = υ x2 + υ y2 .

15

Майер Р. В. Решение физических задач в электронных таблицах Excel

Используется программа 2.2. После ее запуска получается таблица из
четырех столбцов: t , x , y , υ , на основе которой строят траекторию точки и определяют ее скорость υ в момент τ ' = 12,3 с.
Программа 2.2
Private Sub CommandButton1_Click()
a = 0.2: b = 0.15: dt = 0.1
While t < 20
t = t + dt: i = i + 1
x = 5 * (Exp(a * t) - Exp(-a * t)) / 2: y = 3 * t
Cells(i, 1) = t: Cells(i, 2) = x: Cells(i, 3) = y
vx = (x - x1) / dt: vy = (y - y1) / dt
v = Sqr(vx * vx + vy * vy)
Cells(i, 4) = v: x1 = x: y1 = y:
Wend
End Sub

2.9. Точка движется по закону: x(τ ) = 4ch(0,3τ ) (м), y (τ ) = 2 τ (м).
Постройте траекторию движения, определите скорость точки в момент
τ ' = 13,7 с.
Задача решается аналогично предыдущей. Формула для гиперболического косинуса имеет вид: ch(aτ ) = (exp(aτ ) + exp(− aτ )) / 2 .

16

Майер Р. В. Решение физических задач в электронных таблицах Excel
2. Кинематика точки:
двумерное движение

Содержание

4. Механические колебания

3. ДИНАМИКА: ОДНОМЕРНОЕ ДВИЖЕНИЕ
3.1. На точку массы 0,5 кг действует сила (в Н):

 10τ , если 0 ≤ τ < 6,

Fx (τ ) =  10, если 6 ≤ τ < 13,
− 40, если 13 ≤ τ < 20.

Исследуйте движение точки вдоль оси Ox , проанализируйте получившиеся графики зависимостей x(τ ) , υ x (τ ) , a x (τ ) .
Применяются следующие формулы:

a tx+1 = F t +1 / m , υ t +1 = υ t + a t +1∆τ , xt +1 = xt + υ t +1∆τ .
Используется программа 3.1, получающиеся результаты представлены на рис. 3.1. Решите задачу при другой зависимости Fx (τ ) .

Рис. 3.1. Движение точки под действием изменяющейся силы
Программа 3.1
Private Sub CommandButton1_Click()
dt = 0.001: m = 0.5
While t < 20
k = k + 1: t = t + dt: F = -40
If t < 13 Then F = 10: If t < 6 Then F = 10 * t
a = F / m: v = v + a * dt: x = x + v * dt

17

Майер Р. В. Решение физических задач в электронных таблицах Excel

If k Mod 100 = 0 Then
Cells(k / 100, 1) = t: Cells(k / 100, 2) = a
Cells(k / 100, 3) = v: Cells(k / 100, 4) = x / 10
End If
Wend
End Sub

3.2. Тело массой 10 кг начинает скользить вдоль оси Ox под действием постоянной силы 5 Н. Модуль силы сопротивления пропорционален скорости: F = 0,1 + 0,2υ (Н). Постройте график зависимости координаты, скорости и ускорения тела от времени.
3.3. Изучите движение материальной точки, падающей по вертикали в вязкой среде под действием силы тяжести при начальных условиях
x0 = 0 , υ0 x > 0 и υ0 x < 0 . Проанализируйте графики x(τ ) , υ x (τ ) , a x (τ ) .
Докажите, что время подъема камня, брошенного в вязкой среде вертикально вверх, меньше времени спуска.

r

r

r

Ось Ox направим вниз. Запишем формулы: F − rυ = ma ,

a tx+1 = ( Fxt − rυ xt ) / m , υ xt +1 = υ xt + a tx+1∆τ , xt +1 = x t + υ xt +1∆τ .
Используется программа 3.2, получающиеся результаты представлены на рис. 3.2. Можно повторить вычисления при меньшем значении
коэффициента сопротивления среды r . Как при этом изменится предельная скорость, с которой движется тело при τ → ∞ ?
Программа 3.2
Private Sub CommandButton1_Click()
dt = 0.001: m = 1: r = 0.5: F = 10: v = -30
While t < 10
k = k + 1: t = t + dt
a = (F - r * v) / m: v = v + a * dt: x = x + v * dt
If k Mod 100 = 0 Then
Cells(k / 100, 1) = t: Cells(k / 100, 2) = a
Cells(k / 100, 3) = v: Cells(k / 100, 4) = x / 5
End If
Wend
End Sub

18

Майер Р. В. Решение физических задач в электронных таблицах Excel

Рис. 3.2. Расчет падения камня в вязкой среде
3.4. Рассчитайте движение ракеты, удаляющейся по прямой от
Земли, если масса горючего уменьшается со скоростью b =
= ∆m / ∆τ = const . Масса корпуса равна m0 , масса заправленной ракеты

m > m0 . Скорость вытекания газов относительно ракеты равна υ отн . Радиус R Земли известен. Постройте графики зависимостей координаты
x(τ ) , скорости υ x (τ ) , массы m(τ ) ракеты от времени.
Программа 3.3 содержит цикл по времени, в котором вычисляются масса ракеты m , реактивная сила F р = υ отн ∆m / ∆τ , сила притяжения Земли F , ускорение a x , скорость υ x и координата x в моменты τ i = ∆τ ⋅ t : mt +1 = mt − b∆τ ,

F t +1 = GMmt +1 /( x t +1 ) 2 ,

Fрt +1 = υотн (mt − mt +1 ) / ∆τ ,
a tx+1 = ( Fрt +1 − F t +1 ) / m t +1 ,

υ xt +1 = υ xt + a xt +1∆τ ,

xt +1 = xt + υ xt +1∆τ .

19

Майер Р. В. Решение физических задач в электронных таблицах Excel

Рис. 3.3. Результат моделирования движения ракеты
Программа 3.3
Private Sub CommandButton1_Click()
dt = 0.01: x = 50: m1 = 10: b = 1: vot = 4: m0 = 1
While t < 20
k = k + 1: t = t + dt
m1 = m1 - b * dt: m2 = m1 - b * dt:
If m2 < m0 Then b = 0
Fr = vot * (m1 - m2) / dt: F = 500 * m1 / x / x
a = (Fr - F) / m1: v = v + a * dt: x = x + v * dt
If k Mod 10 = 0 Then
Cells(k / 10, 1) = t: Cells(k / 10, 2) = m1 * 10
Cells(k / 10, 3) = v * 10: Cells(k / 10, 4) = x
End If
Wend
End Sub

Получающиеся графики приведены на рис. 3.3. Видно, что по мере вытекания газов масса ракеты с топливом уменьшается до значения m0 , скорость сначала возрастает, а после остановки двигателя
в момент τ ' ≈ 9 с начинает убывать из-за притяжения Земли. Ракета
продолжает удаляться.
3.5. Тележка с телом движется с некоторой скоростью. В тело попадает пуля (рис. 3.4) и застревает в нем (или проходит сквозь него). Как
изменяются скорости тележки и пули в процессе взаимодействия? Постройте графики υ1 (τ ) и υ 2 (τ ) .

20

Майер Р. В. Решение физических задач в электронных таблицах Excel

Рис. 3.4. Движение тела до и после попадания пули
Предположим, что сила взаимодействия пули и тела пропорциональна скорости их относительного движения. Запишем формулы,
позволяющие рассчитать ускорение пули и тела, их скорости и координаты во время взаимодействия:
t +1
t +1
t +1
Fтр
= r | υ 2t − υ1t | , a1t x+1 = Fтр
/ m1 , a2t +x1 = − Fтр
/ m2 ,

υixt +1 = υixt + aixt +1∆τ ,

xixt +1 = xixt + υixt +1∆τ ,

i = 1,2.

Рис. 3.5. Графики зависимостей скоростей тела и пули от времени
Когда пуля находится вне тела, она с ним не взаимодействует.
На рис. 3.5 приведены графики зависимостей скоростей пули и тела
от времени в следующих случаях: 1) пуля и тело двигались навстречу
друг другу, пуля прошла навылет (рис. 3.5.1); 2) пуля догнала тело и
застряла в нем (рис. 3.5.2).
3.6. На шарик массой m1 кладут шарик массой m2 и сбрасывают их
с высоты h на горизонтальную поверхность. Рассчитайте высоту подскока второго шарика h2 после того, как первый шарик отскочит от поверхности и соударится с ним. Постройте график зависимости высоты
h2 от его массы m2 , считая все соударения абсолютно упругими.

21

Майер Р. В. Решение физических задач в электронных таблицах Excel

Рис. 3.6. Падение двух шариков на горизонтальную поверхность
До взаимодействия с горизонтальной поверхностью шарики
имеют скорость υ1 = υ 2 = 2 gh = υ (рис. 3.6.1). После удара о пол первый шарик начинает двигаться вверх с той же по модулю скоростью
υ . Скорости тел после центрального абсолютно упругого удара:

r
r
r 2m2υ 2 + (m1 − m2 )υ1
u1 =
,
m1 + m2

r
r
2m1υ1 + (m2 − m1 )υ 2
r
u2 =
.
m1 + m2

Проекции скоростей:

υ1 y = υ , υ 2 y = −υ , u1 y =

(−3m2 + m1 )υ
(3m1 − m2 )υ
, u2 y =
.
m1 + m2
m1 + m2

Высота подскока h2 = u 22 / 2 g . Используется программа 3.4, результаты ее работы представлены на рис. 3.6.2. В нашем случае m1 = 0,3 кг.
Видно, что, когда m2 / m1 = 3 , первый шарик после соударения со вторым останавливается, передавая ему всю свою энергию.
Программа 3.4
Private Sub CommandButton1_Click()
h = 1: m1 = 0.03: v = Sqr(2 * 9.8 * h)
While m2 < 0.05
i = i + 1: m2 = 0.001 * i: u2 = (3 * m1 - m2) * v / (m1 + m2)
u1 = (-3 * m2 + m1) * v / (m1 + m2): h2 = u2 * u2 / 2 / 9.8
Cells(i, 1) = m2: Cells(i, 2) = u2
Cells(i, 3) = h2: Cells(i, 4) = u1

22

Майер Р. В. Решение физических задач в электронных таблицах Excel

Wend
End Sub

3.7. Две материальные точки совершают одномерное движение в
промежутке между двумя параллельными стенками, находящимися на
расстоянии L = 120 см (рис. 3.7). Между частицами действуют силы отталкивания по закону F = k / r 2 (Н), где k = 1000 Н/см 2 . При ударе о
стенку частица изменяет свою скорость по направлению. Необходимо
рассчитать координаты и скорость частиц в заданный момент времени
τ ' = 2,7 с, построить графики движения x1 (τ ) и x2 (τ ) .

Рис. 3.7. Движение отталкивающихся частиц между стенками
Рассмотрим движение двух частиц m1 и m2 в одномерном сосуде.
Запишем основные уравнения: F = 1000 / r 2 , r =| x1 − x2 | , a1 = − F / m1 ,
a2 = F / m2 . При абсолютно упругом ударе о стенку проекция скорости
частицы на ось Ox изменяет свой знак. Используется программа 3.5,
получающиеся графики x1 (τ ) и x2 (τ ) представлены на рис. 3.8. Координаты и скорости частиц в момент τ ' = 2,7 с определяются из таблицы. Изучите движение системы, когда m1 значительно больше m2 .
Программа 3.5
Private Sub CommandButton1_Click()
dt = 0.005: m1 = 1.2: m2 = 1: x1 = 0: x2 = 100: v1 = 10: v2 = -5
For i = 1 To 50000
t = t + dt: r = Abs(x1 - x2): F = 1000 / r / r: F1 = -F: F2 = F
a1 = F1 / m1: a2 = F2 / m2
v1 = v1 + a1 * dt: v2 = v2 + a2 * dt
x1 = x1 + v1 * dt: x2 = x2 + v2 * dt
If x1 < 0 Then v1 = -v1: If x2 > 120 Then v2 = -v2
If i Mod 50 = 0 Then Cells(i / 50, 1) = t / 50
If i Mod 50 = 0 Then Cells(i / 50, 2) = x1
If i Mod 50 = 0 Then Cells(i / 50, 3) = x2
Next
End Sub

23

Майер Р. В. Решение физических задач в электронных таблицах Excel

Рис. 3.8. Графики x1 (τ ) и x2 (τ ) движения частиц между стенками
3.8. Три материальные точки совершают одномерное движение в
промежутке между двумя стенками, расположенными на расстоянии

L = 120 см (рис. 3.9). Частицы отталкиваются по закону F = k / r 2 (Н), где
k = 1000 Н/см 2 . Рассчитайте координаты и скорости частиц в момент
времени τ ' = 3,9 с, постройте графики движения x1 (τ ) , x2 (τ ) и x3 (τ ) .

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

Рис. 3.10. Движение трех частиц в одномерном сосуде

24

Майер Р. В. Решение физических задач в электронных таблицах Excel

Программа 3.6
Private Sub CommandButton1_Click()
dt = 0.005: m1 = 0.6: m2 = 0.5: m3 = 0.4
x1 = 10: x2 = 50: x3 = 80: v1 = 5: v3 = -5
For i = 1 To 50000
t = t + dt: r12 = Abs(x1 - x2): r23 = Abs(x2 - x3)
F12 = 500 / r12 / r12: F23 = 500 / r23 / r23
a1 = -F12 / m1: a2 = (F12 - F23) / m2: a3 = F23 / m3
v1 = v1 + a1 * dt: v2 = v2 + a2 * dt: v3 = v3 + a3 * dt
x1 = x1 + v1 * dt: x2 = x2 + v2 * dt: x3 = x3 + v3 * dt
If x1 < 0 Then v1 = -v1
If x3 > 120 Then v3 = -v3
If i Mod 50 = 0 Then Cells(i / 50, 1) = t / 50
If i Mod 50 = 0 Then Cells(i / 50, 2) = x1
If i Mod 50 = 0 Then Cells(i / 50, 3) = x2
If i Mod 50 = 0 Then Cells(i / 50, 4) = x3
Next
End Sub

3.9. Между двумя вертикальными стенками движутся две частицы
массами m1 и m2 , скользящие по горизонтальному стержню и связанные
пружиной жесткостью k (рис. 3.11). Начальные координаты и скорости
частиц заданы. При ударе о стенку частица изменяет свою скорость по
направлению. Постройте графики зависимостей x1 (τ ) и x2 (τ ) .

Рис. 3.11. Движение связанных пружиной частиц между стенками
На первую частицу со стороны второй действует сила упруго-

r

r

сти F1x = Fx = k (( x2 − x1 ) − b) . Так как F1 = − F2 , то F2 x = − F1x и ускорения частиц a1x = Fx / m1 , a2 = − Fx / m2 . Используется программа 3.7,
результат моделирования на рис. 3.12. Видно, что система частиц
перемещается как единое целое, отскакивая от одной стенки к другой. При этом частицы совершают свободные колебания относительно общего центра масс. Повторите вычисления при других m1 ,

m2 , k , начальных координатах и скоростях.

25

Майер Р. В. Решение физических задач в электронных таблицах Excel

Программа 3.7
Private Sub CommandButton1_Click()
dt = 0.001: m1 = 1.8: m2 = 1.5: k = 3: b = 30
x1 = 20: x2 = 60: v1 = 4: v2 = -2
For i = 1 To 400000
t = t + dt: F = k * ((x2 - x1) - b): a1 = F / m1: a2 = -F / m2
v1 = v1 + a1 * dt: v2 = v2 + a2 * dt
x1 = x1 + v1 * dt: x2 = x2 + v2 * dt
If x1 < 0 Then v1 = -v1
If x2 > 120 Then v2 = -v2
If i Mod 200 = 0 Then
Cells(i / 200, 1) = t / 200: Cells(i / 200, 2) = x1
Cells(i / 200, 3) = x2
End If
Next
End Sub

Рис. 3.12. Получающиеся графики x1 (τ ) и x2 (τ )
3.10. Три материальные точки совершают одномерное движение в
промежутке между двумя закрепленными частицами, отстоящими друг
от друга на расстояние L . Между всеми пятью частицами действуют
силы отталкивания по закону F = k / r 2 . Рассчитайте координаты и скорости движущихся частиц в момент времени τ ' , постройте графики их
движения x1 (τ ) , x2 (τ ) и x3 (τ ) .

26

Майер Р. В. Решение физических задач в электронных таблицах Excel
3. Динамика:
одномерное движение

Содержание

5. Динамика: двумерное
движение в силовом поле

4. МЕХАНИЧЕСКИЕ КОЛЕБАНИЯ
4.1. Тело совершает затухающие колебания, описывающиеся
уравнением: x(τ ) = A exp(− βτ ) sin(ωτ + ϕ 0 ) . Постройте графики зависимости координаты x = x (τ ) и проекции скорости υ x = υ x (τ ) от времени, а
также фазовую кривую при различных значениях амплитуды A , частоты
ω , начальной фазы ϕ 0 и коэффициента затухания γ .
Решим эту задачу с помощью электронных таблиц, не используя
модуль создания макросов. Значения амплитуды, частоты, начальной
фазы и коэффициента затухания запишем в ячейки H1–H4. При обращении к этим ячейкам следует использовать абсолютную адресацию.
Таблица 4.1
Ячейки
A1, A2, A3, …
B2
C2

Величина

τ
x(τ )
υ x (τ )

Содержимое
0; 0,1; 0,2; 0,3; 0,4; 0,5; 0,6 ...
=$H$1*EXP(-$H$4*A2)*SIN($H$2*A2+$H$3)
=(B3-B2)/0,1

Рис. 4.1. Результат моделирования затухающих колебаний

27

Майер Р. В. Решение физических задач в электронных таблицах Excel

4.2. Пружинный маятник состоит из тела массы m , подвешенного
на пружине с жесткостью k в среде с вязкостью r . Систему выводят из
равновесия и сообщают начальную скорость. Промоделируйте затухающие колебания при различных параметрах системы и начальных
условиях x0 и υ 0 . Постройте фазовую кривую в осях x и υ x .
По второму закону Ньютона: max = − kx − rυ x . Получаем:

m&x& + r&x& + kx = 0 , a tx+1 = (−kx t − rυ xt ) / m ,

υ xt +1 = υ xt + a tx+1∆τ , xt +1 = x t + υ xt +1∆τ .
В столбец A вводят время τ с шагом 0,01 с. Параметры m , k и
r записывают в ячейки F2, G2, H2 соответственно. В ячейку B3 следует ввести “=(-$G$2*D2-$H$2*C2)/$F$2”, в ячейку C3 вводят
“=C2+B3*0,01”, а в ячейку D3 – ”=D2+C2*0,01”. После этого формулы копируют вниз, создавая таблицу. Результаты моделирования представлены на рис. 4.2. Изменяя параметры колебательной системы,
можно убедиться в том, что: 1) при увеличении массы период колебаний растет; 2) при увеличении жесткости период уменьшается;
3) при увеличении коэффициента сопротивления затухание происходит быстрее. Используя значения x(τ ) и υ x (τ ) , постройте фазовую

кривую в осях x и υ .

Рис. 4.2. Затухающие колебания пружинного маятника

28

Майер Р. В. Решение физических задач в электронных таблицах Excel

4.3. Напишите программу, моделирующую вынужденные колебания пружинного маятника. Промоделируйте колебания при различных
параметрах системы m , k , r , частоте ω вынуждающей силы и начальных условиях x0 и υ 0 . Постройте фазовую кривую в осях x и υ x .
По второму закону Ньютона: ma x = − kx − rυ x + Fm sin(ωτ ) ,

F t = Fm sin(ω ⋅τ ) ,

a tx+1 = ( F t − kxt − rυ xt ) / m ,

υ xt +1 = υ xt + a tx+1∆τ , x t +1 = x t + υ xt +1∆τ , τ = t ⋅ ∆τ .
Используется программа 4.1, результаты на рис. 4.3.
Программа 4.1
Private Sub CommandButton1_Click()
m = 0.2: k = 10: r = 0.2
w = 5: dt = 0.005
While t < 20
t = t + dt: i = i + 1: F = 2 * Sin(w * t)
a = (F - k * x - r * v) / m: v = v + a * dt
x = x + v * dt: Cells(i, 1) = t: Cells(i, 2) = a
Cells(i, 3) = v: Cells(i, 4) = x
Wend
End Sub

Рис. 4.3. Вынужденные колебания при различных частотах ω
4.4. Используя решение предыдущей задачи, проведите серию вычислительных экспериментов при различных значениях частоты ω вынуждающей силы и постройте резонансную кривую A = A(ω ) .
Механическим резонансом называется резкое увеличение амплитуды вынужденных колебаний A , когда частота вынуждающей
силы ω совпадает с собственной частотой колебательной системы

ω0 = k / m . Необходимо рассчитать частоту собственных колебаний

ω0 и с помощью программы 4.1 (по таблице или графику) определить
установившуюся амплитуду колебаний A , когда частота вынуждаю-

29

Майер Р. В. Решение физических задач в электронных таблицах Excel

щей силы ω изменяется вблизи

ω0 . На основе полученных результа-

тов постройте график зависимости амплитуды колебаний от частоты вынуждающей силы A = A(ω ) .
4.5. Напишите программу, моделирующую затухающие колебания
пружинного маятника. Изучите затухающие колебания при слабом и
сильном затухании. Постройте фазовую кривую в осях x и υ x .
Для решения этой задачи необходимо в программе 4.1 приравнять вынуждающую силу F (τ ) к 0 и задать начальные значения координаты x0 и скорости υ x 0 .
4.6. К шкиву массой m радиуса R , установленному на горизонтальной оси, с помощью невесомого стержня длиной l прикреплен груз
массой m1 (рис. 4.4.1). На шкив намотана нить, к ее концу привязан груз

m2 . При τ = 0 система покоилась. Рассчитайте угловые ускорение ε ,
скорость ω и перемещение ϕ шкива в последующие моменты времени,

r

r

если в подшипнике действует тормозящий момент M тр = − rω .

Рис. 4.4. К задаче о вращении шкива с грузом
Из законов динамики: Iε = TR − m1 gl sin ϕ − rω , T = m2 g − m2 a , где

I = mR 2 / 2 + m1l 2 – момент инерции шкива с грузом, a = Rε – ускорение тела m2 . Отсюда: Iε = m2 gR − m2 R 2ε − m1 gl sin ϕ − rω . Угловое ускорение шкива равно:

ε=

m2 gR − m1 gl sin ϕ − rω
.
I + m2 R 2

Для моделирования рассматриваемого движения необходимо создать
цикл по времени τ , в котором рассчитывались бы угловые ускорение

30

Майер Р. В. Решение физических задач в электронных таблицах Excel

ε , скорость ω и перемещение ϕ (программа 4.2). Из графиков
(рис. 4.4.2 и 4.4.3) видно, что при заданных параметрах системы ускорение ε и скорость ω колеблются относительно некоторых значений, угловое перемещение ϕ шкива возрастает. Если масса тела m2
мала, то шкив не совершит ни одного оборота, в системе будут происходить затухающие колебания. Изучите эту ситуацию самостоятельно.
Программа 4.2
Private Sub CommandButton1_Click()
m1 = 0.05: m2 = 0.2: m = 0.1: R = 0.03: l = 0.1: dt = 0.02
Iner = m * R * R / 2 + m1 * l * l
While t < 50
t = t + dt: i = i + 1
eps = (m2 * 9.8 * R - m1 * 9.8 * l * Sin(fi) - 0.02 * w)
/ (Iner + m2 * R * R)
w = w + eps * dt: fi = fi + w * dt
Cells(i, 1) = t: Cells(i, 2) = eps
Cells(i, 3) = w: Cells(i, 4) = fi
Wend
End Sub

4.7. На физический маятник действует внешний момент сил, изменяющийся по гармоническому закону M (τ ) = M m sin(ω1τ ) . Напишите
дифференциальное уравнение движения маятника, промоделируйте его
колебания и постройте фазовую кривую.
Вынужденные колебания маятника описываются уравнением:
Iϕ&& + rϕ& + mgl sin ϕ = M sin(ω1τ ) .
Используется программа 4.3. При малых амплитудах sin ϕ ≈ ϕ , поэтому колебания почти гармонические. По мере увеличения амплитуды
они становятся негармоническими. Можно так подобрать амплитуду
внешнего воздействия M m и ее частоту ω1 , что маятник будет переворачиваться через точку подвеса. Из рис. 4.5 видно, что маятник
совершил один полный оборот и колеблется относительно положения равновесия ϕ p = 2π ≈ 6,28 рад. Получающиеся фазовые кривые
представлены на рис. 4.6. При определенных M m и ω1 в спектре колебаний появляются новые гармоники, при этом фазовая кривая содержит несколько циклов (рис. 4.6.2).

31

Майер Р. В. Решение физических задач в электронных таблицах Excel

Рис. 4.5. Вынужденные колебания физического маятника
Программа 4.3
Private Sub CommandButton1_Click()
Iner = 0.1: Mm = 0.6: w1 = 1.9: r = 0.05
m = 0.04: l = 0.1: g = 9.8: dt = 0.005
While t < 50
t = t + dt: i = i + 1
eps = (Mm * Sin(w1 * t) - r * w - m * g * l * Sin(fi)) / Iner
w = w + eps * dt: fi = fi + w * dt
If i Mod 10 = 0 Then Cells(i / 10, 1) = t
If i Mod 10 = 0 Then Cells(i / 10, 2) = fi
If i Mod 10 = 0 Then Cells(i / 10, 3) = w
Wend
End Sub

Рис. 4.6. Фазовые кривые для вынужденных колебаний
С помощью компьютерной модели можно установить, что при
τ → ∞ фазовая кривая стремится к некоторому множеству точек
(замкнутой кривой), которое называется аттрактором. Форма и раз-

32

Майер Р. В. Решение физических задач в электронных таблицах Excel

меры аттрактора определяются свойствами колебательной системы и не зависят от начальных условий. Самостоятельно убедитесь в
этом, повторив моделирование при других начальных условиях.
4.8. Два пружинных маятника с жесткостями пружин k и массами
m1 и m2 соединены упругой связью с жесткостью q (рис. 4.7). Один из
маятников сместили из положения равновесия и отпустили. Промоделируйте движение системы, определите координаты и скорости маятников
в момент времени τ ' = 1,2 с.

Рис. 4.7. Система из двух маятников, связанных пружиной
Если маятники находятся в положении равновесия ( x1 = x2 = 0) ,
то связывающая их пружина не деформирована и не создает силу упругости Fупр = q ( x2 − x1 ) = 0 . Из второго закона Ньютона:

m1a1 = Fупр − k1 x1 − rυ1 , m2 a2 = − Fупр − k2 x2 − rυ 2 , Fупр = q ( x2 − x1 ) .
Программа должна содержать цикл по времени τ , в котором вычисляются ускорения, скорости и координаты маятников по формулам:
t +1
Fупр
= − q( x1t − x2t ) ,
t +1
a1t x+1 = ( Fупр
− rυ1t − kx1t ) / m1 , υ1t x+1 = υ1t x + a1t x+1∆τ , x1t +1 = x1t + υ1t +1∆τ ,

t +1
a2t +x1 = (− Fупр
− rυ 2t x − kx2t ) / m2 , υ 2t +x1 = υ 2t x + a2t +x1∆τ , x2t +1 = x2t + υ 2t +1∆τ .

Используется программа 4.4; на основе результатов вычислений
строят графики x1 (τ ) , x2 (τ ) (рис. 4.8).
Программа 4.4
Private Sub CommandButton1_Click()
m1 = 1: m2 = 1.1: r = 0.2: k = 100: q=10: dt = 0.001: x1 = 0.2
While i < 15000
i = i + 1: F = -q * (x1 - x2)
a1 = (F - r * v1 - k * x1) / m1
a2 = (-F - r * v2 - k * x2) / m2
v1 = v1 + a1 * dt: v2 = v2 + a2 * dt
x1 = x1 + v1 * dt: x2 = x2 + v2 * dt
If i Mod 10 = 0 Then Cells(i / 10, 1) = i * dt

33

Майер Р. В. Решение физических задач в электронных таблицах Excel

If i Mod 10 = 0 Then Cells(i / 10, 2) = x1
If i Mod 10 = 0 Then Cells(i / 10, 3) = x2
Wend
End Sub

Рис. 4.8. Графики x1 (τ ) и x2 (τ ) колебаний системы из двух маятников
4.9. Автоколебательная система состоит из груза массой m , подвешенного на пружине жесткостью k , и клапана, регулирующего поступление энергии от источника. При прохождении грузом положения равновесия ( | x |< 0,5 ) в направлении оси Ox на него действует постоянная сила F в направлении движения. Необходимо рассчитать состояние системы в произвольный момент времени τ ' , построить график автоколебаний и фазовую кривую.
Математическая модель системы:
a x = dυ x / dt = ( Fx − kx − rυ x ) / m , υ x = dx / dτ ,

5, если x < 0,5 и υ x ≥ 0,
F ( x ,υ x ) = 
 0, в противном случае.
Используемая программа 4.5 состоит из цикла по времени τ , в котором рассчитываются значения координаты x , проекции скорости υ x
и ускорения a x в следующий момент времени τ + ∆τ . Результаты
приведены на рис. 4.9. Видно, что фазовая кривая стремится к неко-

34

Майер Р. В. Решение физических задач в электронных таблицах Excel

торому множеству точек фазовой плоскости, называемому аттрактором. Его форма не зависит от начальных условий x0 и υ 0 .
Программа 4.5
Private Sub CommandButton1_Click()
m = 1.1: k = 1: r = 0.05: dt = 0.02
While t < 50
t = t + dt: i = i + 1
If (Abs(x) < 0.5) And (v >= 0) Then F = 5 Else F = 0
a = (F - r * v - k * x) / m
v = v + a * dt: x = x + v * dt : Cells(i, 1) = t
Cells(i, 2) = x: Cells(i, 3) = v: Cells(i, 4) = a
Wend
End Sub

Рис. 4.9. Результаты моделирования автоколебательной системы
4.10. Промоделируйте автоколебательную систему, описывающуюся уравнением Ван-дер-Поля: m&x& + kx = b(1 − x 2 ) x& .
Автоколебательная система состоит из колебательной системы, источника энергии, клапана и цепи обратной связи. Колебательная система воздействует посредством положительной обратной связи на ключевой элемент, регулируя поступление энергии от
источника. Амплитуда колебаний растет до тех пор, пока энергия
потерь за период не будет равна энергии, поступающей от источника за то же время. Получаем уравнения:

a tx+1 = (b(1 − ( x tx ) 2 )υ xt − kxtx ) / m , υ xt +1 = υ xt + a tx+1∆τ , xt +1 = x t + υ xt +1∆τ .
Используется программа 4.6, результаты – на рис. 4.10.
Программа 4.6
Private Sub CommandButton1_Click()
dt = 0.005: m = 0.8: k = 0.5: b = 0.2: x = 0.01
While t < 150
t = t + dt: i = i + 1

35

Майер Р. В. Решение физических задач в электронных таблицах Excel

a = (b * (1 - x * x) * v - k * x) / m
v = v + a * dt: x = x + v * dt
If i Mod 20 = 0 Then
Cells(i / 20, 1) = t: Cells(i / 20, 2) = a
Cells(i / 20, 3) = v: Cells(i / 20, 4) = x: End If
Wend
End Sub

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

Рис. 4.10. Моделирование генератора Ван-дер-Поля

36

Майер Р. В. Решение физических задач в электронных таблицах Excel
4. Механические колебания

Содержание

6. Задача двух тел

5. ДИНАМИКА: ДВУМЕРНОЕ ДВИЖЕНИЕ В СИЛОВОМ ПОЛЕ
5.1. Камень массой m брошен с некоторой начальной скоростью
υ 0 под углом α к горизонту. Рассчитайте траекторию его движения, координаты и скорость в заданный момент времени τ , если со стороны

r

r

воздуха на него действует сила сопротивления FC = − rυ , где r – коэффициент сопротивления воздуха.

Рис. 5.1. Двумерное движение камня в вязкой среде
При отсутствии силы сопротивления камень движется по параболе, а при ее наличии – по более сложной кривой. На рис. 5.1 показаны

r r
r
действующие силы. Получаем: mg + FC = ma ,

a tx+1 = − rυ xt / m , a ty+1 = (−mg − rυ ty ) / m , υ xt +1 = υ xt + a xt +1∆τ ,

υ ty+1 = υ ty + a ty+1∆τ , xt +1 = x t + υ xt +1∆τ , y t +1 = y t + υ ty+1∆τ .
Программа 5.1 содержит цикл по времени t , в котором рассчитываются ускорение, проекции скорости и координаты камня в следующий
момент времени t + 1, а результаты вычислений выводятся в таблицу. Проекции скоростей в начальный момент времени и коэффициент
сопротивления считываются программой из ячеек I1, I2, I3. На основе
результатов вычислений строится траектория y = y (x ) (рис. 5.2).
При r = 0 камень движется по параболе.
Программа 5.1
Private Sub CommandButton1_Click()
dt = 0.005: m = 1: r = Cells(3, 9)
vx = Cells(1, 9): vy = Cells(2, 9)

37

Майер Р. В. Решение физических задач в электронных таблицах Excel

While y >= 0
t = t + dt: ay = (-10 * m - r * vy) / m
ax = -r * vx / m: vx = vx + ax * dt
vy = vy + ay * dt: x = x + vx * dt
y = y + vy * dt: i = i + 1
If i Mod 5 = 0 Then
Cells(i / 5, 1) = t: Cells(i / 5, 2) = ax: Cells(i / 5, 3) = vx
Cells(i / 5, 4) = ay: Cells(i / 5, 5) = vy: Cells(i / 5, 6) = x
Cells(i / 5, 7) = y: End If
Wend
End Sub

Рис. 5.2. Движение камня с учетом сопротивления воздуха
5.2. Пользуясь таблицей, полученной при решении задачи 5.1, определите и сравните время подъема τ п , время спуска τ с и общее время

τ движения камня в случае, когда: 1) коэффициент сопротивления r = 0 ;

2) r = 0,5 − 1,5 .

Чтобы найти время подъема следует учесть, что в наивысшей
точке траектории проекция скорости на ось y меняет свой знак на
противоположный. При вычислении общего времени движения τ дв расчеты производятся до тех пор, пока y не станет меньше нуля. При
наличии сопротивления воздуха время подъема τ п меньше времени
спуска τ с = τ дв − τ п . Если r = 0 , то τ с = τ п .
5.3. Камень массой m брошен с некоторой начальной скоростью
под углом α к горизонту. Рассчитайте траекторию его движения, коор-

38

Майер Р. В. Решение физических задач в электронных таблицах Excel

динаты и скорость в последовательные моменты времени, если вдоль

r

r

оси Ox дует ветер со скоростью υ В . Учтите, что FC = − rυотн , где
r
r r
υотн = υ − υ В – скорость камня относительно воздуха.

5.4. На диске, вращающемся с угловой скоростью ω , установлен
маятник Фуко, состоящий из небольшого шарика A на нити. Маятник колеблется в вертикальной плоскости. Рассчитайте траекторию движения
конца маятника A в системе отсчета, связанной с диском.

Рис. 5.3. К задаче о расчете траектории движения маятника Фуко
Если математический маятник совершает малые колебания, то
возвращающая сила пропорциональна его смещению из центра диска
O. Поэтому его можно заменить материальной точкой, способной
двигаться в горизонтальной плоскости, на которую действует ква-

r

r

зиупругая сила F упр = − kr . Относительно неинерциальной системы
отсчета, связанной с вращающимся диском, на материальную точку

r
r
A также действуют центробежная сила Fц = ω 2 r и сила Кориолиса
r
r
r
Fк = 2m[υ от , ω ] (рис. 5.3.1). Будем считать, что воздух вращается
вместе с диском с той же угловой скоростью ω , поэтому действующая на колеблющийся маятник сила сопротивления равна

r
r
Fc = − rcυ от . Из основного закона динамики относительное ускорение
r t +1 r t +1 r t +1 r t +1
r t +1
в момент t + 1 равно: a от
= ( F упр
+ Fц + Fк + Fc ) / m . Относиr t +1 r t
r t +1
тельная скорость и радиус-вектор точки A: υ от
= υ от + a от
∆τ ,
r t +1
r t +1 r t
rот
= rот + υ от
∆τ . Получаем уравнения:
r = ( x 2 + y 2 ) 0,5 , cos α = x / r , sin α = y / r , υ от = (υ x2 + υ y2 ) 0,5 ,

39

Майер Р. В. Решение физических задач в электронных таблицах Excel

cos β = υ x / υ от , sin β = −υ y / υ от , Fк = 2mυ от ω , Fкx = − Fк sin β ,
Fкy = − Fк cos β , F упр, x = −kx , F упр, y = −ky , Fц , x = ω 2 x , Fц , y = ω 2 y ,
t +1
t +1
t +1
t +1
Fc, x = −rcυ x , Fc, y = − rcυ y , a tx+1 = ( Fупр
, x + Fц , x + Fк , x + Fc , x ) / m ,
t +1
t +1
t +1
t +1
t +1
t
t +1
a ty+1 = ( Fупр
, y + Fц , y + Fк , y + Fc , y ) / m , υ x = υ x + a x ∆τ ,

υ ty+1 = υ ty + a ty+1∆τ , xt +1 = xt + υ xt +1∆τ , y t +1 = y t + υ ty+1∆τ .
Используется программа 5.2, результаты – на рис. 5.3.2 и 5.4 при различных rc и ω .
Программа 5.2
Private Sub CommandButton1_Click()
dt = 0.0005: m = 0.2: k = 10: w = 2.5: rr = 0.08: x = 0.2
While t < 10
i = i + 1: t = i * dt
r = Sqr(x * x + y * y): cosa = x / r: sina = y / r
votx = vx + w * r * sina: voty = vy - w * r * cosa
vot = Sqr(votx * votx + voty * voty)
cosb = votx / vot: sinb = -voty / vot:
Fk = 2 * m * vot * w: Fkx = -Fk * sinb: Fky = -Fk * cosb
Fux = -k * x: Fuy = -k * y: Fcx = w * w * x: Fcy = w * w * y
ax = (Fux + Fcx + Fkx - rr * vx) / m
ay = (Fuy + Fcy + Fky - rr * vy) / m
vx = vx + ax * dt: vy = vy + ay * dt
x = x + vx * dt: y = y + vy * dt
If i Mod 20 = 0 Then Cells(i / 20, 1) = x: Cells(i / 20, 2) = y
End If
Wend
End Sub

Рис. 5.4. Получающиеся траектории движения маятника Фуко

40

Майер Р. В. Решение физических задач в электронных таблицах Excel

Рис. 5.5. К задаче о движении точки в поле центральной силы
5.5. Планета (или комета) движется в гравитационном поле Солнца. Их массы m и M . Рассчитайте траекторию движения планеты, ее
координаты и проекции скорости в последовательные моменты времени. Решите задачу при различных массах m и M , начальных координатах x0 , y0 и скоростях υ x 0 , υ y 0 .
Из законов механики следует (рис. 5.5):

F = GmM / r 2 ,

r
r
r
r
F = F ⋅ [−( x / r ) ⋅ i − ( y / r ) ⋅ j ] = ma ,

Fx = − F cos α = − F ⋅ x / r ,

Fy = − F sin α = − F ⋅ y / r ,

a tx+1 = Fxt / m , υ xt +1 = υ xt + a xt +1∆τ , xt +1 = xt + υ xt +1∆τ ,
a ty+1 = Fyt / m , υ ty+1 = υ ty + a ty+1∆τ , y t +1 = y t + υ ty+1∆τ .
Программа для расчета движения планеты содержит цикл, в котором пересчитываются проекции действующей силы, ускорения,
скорости, координаты в последовательные моменты времени
τ i = i ⋅ ∆τ и рисуется траектория (рис. 5.6).
Программа 5.3
Private Sub CommandButton1_Click()
dt = 0.0005: m = 1: x = -20: y = 0: vy = 8: r1 = Abs(x)
For i = 1 To 25000
t = i * dt: r = Sqr(x * x + y * y): cosa = x / r: sina = y / r
F = -2000 * m / r / r: ax = F * cosa / m: ay = F * sina / m
vx = vx + ax * dt: vy = vy + ay * dt
x = x + vx * dt: y = y + vy * dt
If i Mod 50 = 0 Then
Cells(i / 50, 1) = t: Cells(i / 50, 2) = x
Cells(i / 50, 3) = y
End If
Next
End Sub

41

Майер Р. В. Решение физических задач в электронных таблицах Excel

Рис. 5.6. Траектория движения планеты (кометы) вокруг звезды
5.6. В гравитационном поле Земли движется тело. Его начальная
r
скорость υ0 перпендикулярна направлению на центр Земли. Убедитесь,
что при малых скоростях υ0 тело движется по замкнутой эллиптической
орбите, а при увеличении υ0 эллипс становится более вытянутым, превращается в параболу (критический случай), а затем в гиперболу.
5.7. Промоделируйте движение точки в поле центральной силы
притяжения, не подчиняющейся закону обратных квадратов. Рассмотрите случаи: 1) F = k / r ; 2) F = k / r ; 3) F = k / r 3 / 2 .
По теореме Бертрана, частица движется по замкнутой траек-

r

r

тории в двух случаях: 1) в поле квазиупругой силы F = −kr ; 2) в поле
силы притяжения, которая обратно пропорциональна квадрату расстояния r до центра O. Для решения задачи нужно в программу 5.3
внести незначительные изменения. На рис. 5.7 представлены результаты расчетов движения частицы в центрально-симметричных полях, не подчиняющихся закону обратных квадратов. Траекторией является незамкнутая кривая (розетка).

42

Майер Р. В. Решение физических задач в электронных таблицах Excel

Рис. 5.7. Движение точки в полях: 1) F = k / r ; 2) F = k / r 3 / 2
5.8. Измените программу так, чтобы она моделировала движение
частицы в центральном поле силы Fr (r ) = k1 / r 4 − k 2 / r 2 .
Подобное уравнение описывает силу взаимодействия между
атомами двухатомной молекулы: при больших r преобладают силы
притяжения Fr < 0 , а при малых – силы отталкивания Fr > 0 . Промоделируйте несколько ситуаций при различных начальных условиях и
коэффициентах k1 и k2 , получающиеся траектории движения
(рис. 5.8) нарисуйте в тетради.

Рис. 5.8. Движение в центральном поле F (r ) = 5 ⋅105 / r 4 − 5000 / r 2
5.9. По направлению к массивному положительно заряженному ядру движутся альфа-частицы (опыт Резерфорда). Рассчитайте траекторию движения одной альфа-частицы, учитывая, что между ней и ядром
действуют силы отталкивания F = k / r 2 .

43

Майер Р. В. Решение физических задач в электронных таблицах Excel

Рис. 5.9. Движение частицы в центральном поле сил отталкивания
Используется программа 5.4. Вначале задается прицельный параметр y0 = 30. Действуют силы отталкивания, поэтому в программе
сила F положительна. Траекториями частиц являются гиперболы
(рис. 5.9).
Программа 5.4
Private Sub CommandButton1_Click()
dt = 0.001: m = 1: x = -100: y = 25: vx = 15
For i = 1 To 15000
t = i * dt: r = Sqr(x * x + y * y): cosa = x / r: sina = y / r:
F = 2000 / r / r: ax = F * cosa / m: ay = F * sina / m
vx = vx + ax * dt: vy = vy + ay * dt:
x = x + vx * dt: y = y + vy * dt:
If i Mod 100 = 0 Then
Cells(i / 100, 1) = t: Cells(i / 100, 2) = x
Cells(i / 100, 3) = y: ‘Cells(i / 100, 4) = vy / vx
End If
Next
End Sub

5.10. В опыте Резерфорда к неподвижному ядру атома золота подлетает альфа-частица и, отталкиваясь от него, отклоняется на некоторый

44

Майер Р. В. Решение физических задач в электронных таблицах Excel

угол от первоначального направления движения (рис. 5.9). Проведите
серию вычислительных экспериментов при различных прицельных параметрах y0 и по полученным данным постройте график зависимости
угла отклонения от прицельного параметра α ( y 0 ) .
Программу 5.4 следует дополнить так, чтобы она вычисляла угол
отклонения направления движения частицы после взаимодействия.
Если активизировать оператор Cells(i/100,4)=vy/vx, то после
запуска программы в столбце D появится тангенс угла α между векr
тором скорости υ и осью Ox (рис. 5.9). Чтобы получить значение угла в градусах, в ячейку E1 записывают “=ATAN(D1)*180/3,1415”, а затем копируют эту формулу в остальные ячейки столбца E. Искомое
значение угла α будет находиться в нижних строчках таблицы. Повторите вычисления при y0 = 0, 1, 2, 3, 5, 10, 15, 20, 25 и т. д., каждый
раз определяя угол отклонения частицы α ; постройте график α ( y 0 )
стандартными средствами Excel. Получится кривая, которая с ростом y0 монотонно убывает от 180 до 0.
5.11. Предложите метод расчета секториальной скорости планеты,
движущейся по эллиптической орбите. Убедитесь в том, что секториальная скорость планеты остается постоянной (второй закон Кеплера).
Секториальной скоростью ω s называется отношение площади

r

сектора ∆S , заметаемого радиус-вектором планеты r , к соответствующему промежутку времени ∆τ . Пусть за время ∆τ планета
перемещается из A( x , y ) в B( x1 , y1 ) (рис. 5.10.1). Радиус-вектор пла-

r

неты r заметает площадь треугольника OAB. Получаем:

| BC |2 =| AB |2 − | AC |2 , | AB |2 = ( x − x1 ) 2 + ( y − y1 ) 2 ,

| AC |= r − r1 , | BC |= ( x − x1 ) 2 + ( y − y1 ) 2 − (r − r1 ) 2 ,
∆S = r | BC | / 2 , секториальная скорость: ωs = ∆S / ∆τ .
На рис. 5.10.2 представлены графики зависимости расстояния r
от планеты до Солнца, модуля линейной скорости υ и секториальной
скорости ω s от времени. Видно, что секториальная скорость не изменяется, это и подтверждает второй закон Кеплера.

45

Майер Р. В. Решение физических задач в электронных таблицах Excel

Рис. 5.10. К подтверждению второго закона Кеплера
Программа 5.5
Private Sub CommandButton1_Click()
dt = 0.0005: m = 1: x = -20: y = 0: vy = 6: r1 = Abs(x)
For i = 1 To 25000
t = i * dt: r = Sqr(x * x + y * y): cosa = x / r: sina = y / r
F = -2000 * m / r / r: ax = F * cosa / m: ay = F * sina / m
vx = vx + ax * dt: vy = vy + ay * dt
x = x + vx * dt: y = y + vy * dt
If i Mod 50 = 0 Then
Cells(i / 50, 1) = t: Cells(i / 50, 2) = x
Cells(i / 50, 3) = y: ‘Cells(i / 50, 4) = dS / 50 / dt
‘LL = (x - x1) * (x - x1) + (y - y1) * (y - y1)
‘BC = Sqr(LL - (r - r1) * (r - r1))
‘dS = r * BC / 2: x1 = x: y1 = y: r1 = Sqr(x1 * x1 + y1 * y1)
End If
Next
End Sub

5.12. Подтвердите третий закон Кеплера: для любой планеты Солнечной системы отношение квадрата периода обращения к кубу большой полуоси ее орбиты остается постоянным: T 2 / a 3 = const . Проведите
серию вычислительных экспериментов, изменяя начальную скорость
планеты υ0 y и расстояние x0 до Солнца и определяя период вращения

T и большую полуось a орбиты.
Если начальные условия задать так: x0 = –20, y0 = 0, υ x = 0,

υ y = 6 (все величины в условных единицах), то планета начнет свое
движение из точки А( x0 ,0) , лежащей на оси Ox левее нуля, со скоростью, параллельной оси Oy (рис. 5.11). Через половину периода

t1 = T / 2 она оказывается в точке B( x1 ,0). Большая ось AB эллипса
имеет длину x1 − x0 и совпадает с осью Ox . Поэтому большая полуось
орбиты a = ( x1 − x0 ) / 2 . Значение x1 и соответствующий ему момент

46

Майер Р. В. Решение физических задач в электронных таблицах Excel

времени t1 могут быть найдены из таблицы, получающейся в результате работы программы 5.3 или 5.5. Период обращения планеты
T = 2t1 . Выполнив 5–10 численных экспериментов при различных начальных координатах x0 и скоростях υ y , заполните таблицу 5.1. Убедитесь, что коэффициент k = T 2 / a 3 остается постоянным.

Рис. 5.11. Расчет движения планеты вокруг звезды
Таблица 5.1

5.13. Несколько раз промоделируйте движение планеты при одинаковых x0 и различных υ 0 = υ 0 y , в каждом случае определяя большую
полуось a = ( x1 − x0 ) / 2 орбиты, минимальное расстояние q до Солнца
(оно равно | x0 | или x1 ) и эксцентриситет e = 1− q / a . Результаты оформите в виде таблицы, содержащей υ 0 y , x0 , x1 , a , q , e , проанализируйте
ее и сделайте вывод.

47

Майер Р. В. Решение физических задач в электронных таблицах Excel
5. Динамика: двумерное
движение в силовом поле

Содержание

7. Динамика твердого тела

6. ЗАДАЧА ДВУХ ТЕЛ
6.1. Две частицы массами m1 и m2 движутся навстречу друг к другу
со скоростями υ1 и υ2 . Промоделируйте абсолютно упругое взаимодействие этих частиц, если между ними действуют: а) силы отталкивания; б)
силы притяжения. В обоих случаях F = k / r 2 .

Рис. 6.1. Нецентральное взаимодействие двух частиц
Расстояние между частицами r = (( x1 − x2 ) 2 + ( y1 − y2 ) 2 ) 0,5 , где

( x1, y1 ) и ( x2 , y2 ) – их координаты (рис. 6.1). Проекции силы отталкивания: F1x = − F cos α , F1 y = − F sin α , F2 x = − F1x , F2 y = − F1 y где F = k / r 2 ,
cos α = ( x2 − x1 ) / r , sin α = ( y2 − y1 ) / r . Проекции ускорений, скоростей и
координаты частицы равны:

aixt +1 = Fixt +1 / mi , υixt +1 = υ ixt + aixt +1∆τ , xit +1 = xit + υixt +1∆τ ,
aiyt +1 = Fiyt +1 / mi , υiyt +1 = υ iyt + aiyt +1∆τ , yit +1 = yit + υ iyt +1∆τ , i = 1,2.
Для решения этой задачи двух тел используется программа 6.1. Результаты моделирования в случае притяжения и отталкивания частиц с массами m1 / m2 = 1 / 3 представлены на рис. 6.2. Видно, что при
нецентральном упругом взаимодействии двух частиц, между которыми действуют силы отталкивания, частицы движутся по гиперболам. Чтобы получить правильные траектории, необходимо использовать одинаковый масштаб по осям Ox и Oy . Эта компьютерная модель также позволяет изучить центральный удар.

48

Майер Р. В. Решение физических задач в электронных таблицах Excel

Программа 6.1
Private Sub CommandButton1_Click()
dt = 0.002: m1 = 1: m2 = 3: x1 = -100: x2 = 100
y1 = 18: v1x = 6: v2x = -4
For i = 1 To 20000
t = t + dt
r = Sqr((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2))
sinal = (y2 - y1) / r: cosal = (x2 - x1) / r
F = -1000 / r / r: F1x = -F * cosal: F2x = F * cosal
F1y = -F * sinal: F2y = F * sinal
a1x = F1x / m1: a1y = F1y / m1
a2x = F2x / m2: a2y = F2y / m2
v1x = v1x + a1x * dt: v1y = v1y + a1y * dt
v2x = v2x + a2x * dt: v2y = v2y + a2y * dt
x1 = x1 + v1x * dt: y1 = y1 + v1y * dt
x2 = x2 + v2x * dt: y2 = y2 + v2y * dt
If i Mod 50 = 0 Then Cells(i / 50, 1) = t / 50
If i Mod 50 = 0 Then Cells(i / 50, 2) = x1
If i Mod 50 = 0 Then Cells(i / 50, 3) = y1
If i Mod 50 = 0 Then Cells(i / 50, 4) = x2
If i Mod 50 = 0 Then Cells(i / 50, 5) = y2
Next
End Sub

Рис. 6.2. Взаимодействие двух частиц c m2 / m1 = 3 :
1 – притяжение; 2, 3 – отталкивание
6.2. К неподвижной частице массой m1 подлетает частица массой

r
m2 , движущаяся со скоростью υ 2 . Частицы взаимодействуют с силами
отталкивания Fr = k / r 2 , удар абсолютно упругий и нецентральный. Необходимо рассчитать движение частиц, определить их углы разлета по
отношению к первоначальному направлению движения при различных
прицельных параметрах.
Используемая программа 6.2 содержит цикл по времени, в котором, исходя из координат и скоростей частиц в момент t , вычисляются их новые значения в момент t + 1 и записываются в таблицу.

49

Майер Р. В. Решение физических задач в электронных таблицах Excel

Программа 6.2
Private Sub CommandButton1_Click()
dt = 0.0005: m1 = 1: m2 = 3: x1 = -100: x2 = 0
y1 = 12: v1x = 6: v2x = 0
For i = 1 To 70000
r = Sqr((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2))
sinal = (y2 - y1) / r: cosal = (x2 - x1) / r
F = 200 / r / r: F1x = -F * cosal: F2x = F * cosal
F1y = -F * sinal: F2y = F * sinal: a1x = F1x / m1
a1y = F1y / m1: a2x = F2x / m2: a2y = F2y / m2
v1x = v1x + a1x * dt: v1y = v1y + a1y * dt
v2x = v2x + a2x * dt: v2y = v2y + a2y * dt
x1 = x1 + v1x * dt: y1 = y1 + v1y * dt
x2 = x2 + v2x * dt: y2 = y2 + v2y * dt
tga1 = v1y / v1x: tga2 = v2y / v2x: t = t + dt
If i Mod 200 = 0 Then
Cells(i / 200, 1) = t / 200: Cells(i / 200, 2) = x1
Cells(i / 200, 3) = y1: Cells(i / 200, 4) = x2
Cells(i / 200, 5) = y2: Cells(i / 200, 6) = tga1
Cells(i / 200, 7) = tga2: End If: Next
End Sub

Рис. 6.3. Движение частиц при абсолютно упругом ударе
При нецентральном упругом взаимодействии частицы движутся
по гиперболам (рис. 6.3.1). Можно провести серию вычислительных
экспериментов, изменяя прицельный параметр ρ = y1 (0) и определяя
r

r

углы α1 и α 2 между осью Ox и скоростями υ1 и υ 2 частиц после взаимодействия. Для этого в столбцы F и G вводятся значения
tgα1 = υ1 y / υ1x и tgα 2 = υ 2 y / υ 2 x , исходя из которых определяются значения α1 и α 2 в градусах. Получающиеся графики α1 ( ρ ) и α 2 ( ρ ) представлены на рис. 6.3.2.

50

Майер Р. В. Решение физических задач в электронных таблицах Excel

6.3. Используя решение задачи 6.2, убедитесь в том, что частицы с
равными массами после абсолютно упругого удара разлетаются под углом 90 0 .
6.4. Две звезды массами m1 и m2 притягиваются друг к другу силами гравитационного притяжения Fr = − k / r 2 и вращаются вокруг общего
центра масс. Рассчитайте траектории их движения (рис. 6.4).

Рис. 6.4. Вращение двух звезд вокруг общего центра масс
6.5. Частица m1 движется в поле притяжения частицы m2 , причем

m2 >> m1 . Рассчитайте траектории движения частиц, если действуют
силы притяжения Fr = −k / r a , где a > 0 и a ≠ 2 .
Получающиеся результаты представлены на рис. 6.5. Видно,
что более легкая частица движется по незамкнутой траектории,
называемой розеткой.

Рис. 6.5. Движение частицы в центрально-симметричном поле Fr (r )
более массивной частицы: 1) F = k / r ; 2) F = k / r ; 3) F = k / r 3 / 2

51

Майер Р. В. Решение физических задач в электронных таблицах Excel
6. Задача двух тел

Содержание

8. Аналитическая механика

7. ДИНАМИКА ТВЕРДОГО ТЕЛА
7.1. Имеется неоднородный шар радиусом R = 1 м, плотность которого задана уравнением:

ρ (r ) = 103 /(r + 0,2) кг/м 3 , где r – расстояние

до центра O в метрах. Необходимо найти его массу.
Разобьем шар на N = 200 элементарных масс ∆mi , имеющих
форму сферических слоев радиусом ri и толщиной ∆r . Объем каждого
слоя ∆Vi = 4π ⋅ ri2 ∆r , а его масса ∆mi = ρ ( ri ) ∆Vi = 10 ∆Vi /( ri + 0,2) .
Общая масса шара
3

103
m ≈ ∑ ∆mi =∑ ρ (ri )∆Vi = ∑
4π ⋅ ri2 ∆r .
i =1
i =1
i =1 ri + 0,2
N

N

N

Точное значение массы шара равно интегралу:
R

4π ⋅103 2
m = ∫ ρ (r )dV = ∫
r dr .
r
+
0
,
2
V
0
Для решения этой задачи численным методом в среде Excel создается таблица из 200 строк. В столбец A следует ввести значения

ri = 0,005i , где i = 0, 1, 2, ... В столбце B вычисляется элементарный
объем ∆Vi сферического слоя толщиной ∆r = 0,05 м, для этого в ячейку B2 записывают формулу ‘=4*3,1415*A2*A2*0,05’. В ячейку C2 вводят
‘=1000/(A2+0,2)’, в ячейку D2 записывают ‘=B2*C2’, в ячейку E3 вводят
‘=D3+E2'. После этого содержимое ячеек B2, C2, D2 и E3 копируют
вниз, заполняя все 200 строк таблицы. В последнем столбце суммируются все элементарные массы ∆mi . Искомое значение массы шара
находится в самой нижней ячейке столбца E и равно 46966 кг.
7.2. Имеется неоднородный шар радиусом R = 0,4 м, плотность которого равна

ρ (r ) = 2 ⋅103 /(r 2 + 0,15) кг/м 3 , где r – расстояние до цен-

тра O в метрах. Вычислите массу шара.
7.3. Определите момент инерции неоднородного диска
относительно оси Oz , проходящей через центр масс перпендикулярно
диску. Толщина диска 5 см, радиус 0,5 м, плотность равна

ρ (r ) = 4 ⋅103 /( r + 0,12) кг/м 3 , где r – расстояние до оси Oz в метрах.

52

Майер Р. В. Решение физических задач в электронных таблицах Excel

7.4. Параллельно оси Ox лежит неоднородный стержень диаметром 1,5 см, концы которого имеют координаты 0 и l = 0,8 м. Плотность
стержня равна ρ ( x ) = 6,5 ⋅ 103 (1 + x / 2) кг/м , где x измеряется в метрах и
3

отсчитывается от левого конца. Определите массу стержня.
7.5. Найдите положение центра масс стержня, рассмотренного в
предыдущей задаче.
7.6. Определите момент инерции стержня относительно оси, проходящей через его левый конец перпендикулярно стержню. Параметры
стержня взять из задачи 7.4.
7.7. Имеется прямоугольная пластина размером 1,2 х 0,7 м 2 , толщина которой зависит от координат так: h( x, y ) = 10 −3 (2 + 0,8 x + 0,6 y ) м.
3

Определите массу пластины, если ее плотность 7800 кг/м .
7.8. Определите положение центра масс пластины, рассмотренной
в предыдущей задаче.
7.9. Невесомый жесткий стержень, соединяющий две материальные точки m1 и m2 , устанавливают на горизонтальной поверхности так,
чтобы он образовывал некоторый угол α 0 с вертикалью, и отпускают.
Необходимо рассчитать координаты его концов при падении в последовательные моменты времени.

Рис. 7.1. К задаче о падающем стержне
Рассмотрим систему из двух частиц m1 и m2 , связанных невесомым упругим стержнем с жесткостью k (рис. 7.1). Длина стержня в
53

Майер Р. В. Решение физических задач в электронных таблицах Excel

недеформированном состоянии равна l0 , при его сжатии на частицы
A и B действуют силы упругости F1 и F2 . Будем считать, что нижний
конец стержня A скользит по горизонтальной поверхности, не отрываясь от нее ( y1 = 0 ). При этом на него действует сила вязкого трения FТР , пропорциональная скорости и направленная в сторону, противоположную движению. Проекции сил, действующих на материальные точки m1 и m2 , вычисляются из формул:

x2 − x1
− rυ1x ,
l
y −y
= F sin α − m2 g = F 2 1 − m2 g .
l

F1 = F2 = F = k (l0 − l ) , R1x = − F cos α − rυ1x = − F
R2 x = F cos α =

x2 − x1
, R2 y
l

При этом y1 = 0 . Ускорения и скорости частиц вычисляются так:

a1t x+1 = R1tx+1 / m1 , a2t +x1 = R2t +x1 / m2 , a2t +y1 = R2t +y1 / m2 ,

υ1t x+1 = υ1t x + a1t x+1∆τ , υ 2t +x1 = υ 2t x + a2t +x1∆τ , υ 2t +y1 = υ 2t y + a2t +y1∆τ ,
x1t +1 = x1t + υ1t x+1∆τ , x2t +1 = x2t + υ 2t +x1∆τ , y2t +1 = y2t + υ 2t +y1∆τ .

Рис. 7.2. Результаты моделирования падения стержня

54

Майер Р. В. Решение физических задач в электронных таблицах Excel

В программе 7.1 в последовательные моменты времени пересчитываются силы, действующие на материальные точки, их ускорения,
скорости и координаты. Результаты моделирования падения стержня представлены на рис. 7.2. При необходимости можно рассчитать
траекторию движения центра масс стержня, зависимость его угла
наклона от времени и т. д.
Программа 7.1
Private Sub CommandButton1_Click()
m1 = 1: m2 = 2: k = 100: R = 0.02: l0 = 300: dt = 0.001
x1 = 0: y1 = 0: x2 = 1: y2 = 300
For i = 1 To 100000
l = Sqr((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2))
F = k * (l0 - l): Fx1 = -F * (x2 - x1) / l - R * vx1
Fx2 = F * (x2 - x1) / l: Fy2 = F * (y2 - y1) / l - 9.8 * m2
vx1 = vx1 + Fx1 / m1 * dt: vy1 = vy1 + Fy1 / m1 * dt
vx2 = vx2 + Fx2 / m2 * dt: vy2 = vy2 + Fy2 / m2 * dt
x1 = x1 + vx1 * dt: x2 = x2 + vx2 * dt: y2 = y2 + vy2 * dt
t = t + dt
If y2 < 1 Then
vx2 = 0: vy2 = 0: End If
If i Mod 500 = 0 Then
Cells(i / 500, 1) = t
Cells(i / 500, 2) = x1: Cells(i / 500, 3) = y1
Cells(i / 500, 4) = x2: Cells(i / 500, 5) = y2
Cells(i / 500, 6) = ((y2 - y1) / (x2 - x1)): End If
Next
End Sub

7.10. Жесткий невесомый стержень с двумя материальными точками на концах прислоняют к стене так, чтобы он составлял с горизонталью угол α 0 , и отпускают. Рассчитайте падение стержня, если известно,
что его концы скользят по поверхностям стены и пола, не отрываясь от
них, а действующая на них сила трения пропорциональна скорости.
Рассмотрим систему из двух частиц массами m1 и m2 , соединенных невесомым упругим стержнем жесткостью k и длиной l0 в недеформированном состоянии. Эти точки как бы скользят по вертикальной и горизонтальной направляющим, при этом на них действуют силы вязкого трения, прямо пропорциональные скорости и направленные в сторону противоположную движению (рис. 7.3). Найдем
проекции равнодействующих всех сил, приложенных к точкам A и B:
F1 = F2 = F = k (l0 − l ) , x1 = 0 , y 2 = 0 ,

R1 y = F sin α − m1 g − rυ1 y = F

55

y1 − y 2
− m1 g − rυ1 y ,
l

Майер Р. В. Решение физических задач в электронных таблицах Excel

R2 x = F cos α − rυ 2 x = F

x2 − x1
− rυ 2 x .
l

υ1t y+1 = υ1t y + R1t y+1∆τ / m1 , υ 2t +x1 = υ 2t x + R2t +x1∆τ / m2 ,
y1t +1 = y1t + υ1t +y 1∆τ ,

x2t +1 = x2t + υ 2t +x1∆τ .

Рис. 7.3. Падение стержня, приставленного к стене
Задача решается так же, как предыдущая. Используется программа,
содержащая цикл по времени, в котором вычисляются силы, действующие на материальные точки, определяются их ускорения, скорости и координаты, а также координаты центра масс.
7.11. На горизонтальной поверхности покоится кольцо массой mк , к
внутренней стороне которого прикреплен груз mг . Расстояние l =| OM |
от оси кольца до его центра масс известно. Кольцо смещают из положения равновесия и отпускают. Получите графики колебаний при разных
амплитудах.
Рассмотрим качение кольца по горизонтальной поверхности
(рис. 7.4.1). Запишем формулы для расчета расстояния lC от центра
кольца O до центра масс C, момента M силы тяжести, момента
инерции I относительно мгновенной оси вращения A и углового ускорения ε z : lC = mг l /(mг + mк ) , d 2 =| AM |2 = R 2 + l 2 − 2 Rl cos ϕ ,

M z = −(mг + mк ) glC sin ϕ z ,

I = 2mк R 2 + mг d 2 ,

ε z = M z / I = −(mг + mк ) glC sin ϕ z /(2mк R 2 + mг d 2 ) .
56

Майер Р. В. Решение физических задач в электронных таблицах Excel

r

Здесь учитывается, что плечо силы тяжести (mг + mк ) g равно
lC sin ϕ , а момент инерции кольца по теореме Штейнера относительно оси вращения A равен 2mк R 2 .

Рис. 7.4. Нелинейные колебания кольца
со смещенным центром тяжести
Зная момент силы и момент инерции, можно определить угловое ускорение тела в последовательные моменты времени, вычислить угловую скорость и угол поворота, а также горизонтальную координату центра кольца. При этом используется программа 7.2, результаты представлены на рис. 7.5. Чтобы промоделировать качение кольца со смещенным центром масс, следует задать начальную скорость.

Рис. 7.5. Графики колебаний кольца при разных амплитудах

57

Майер Р. В. Решение физических задач в электронных таблицах Excel

Программа 7.2
Private Sub CommandButton1_Click()
dt = 0.005: g = 9.8: R = 60: l = 55
mgr = 10: x = -10: fi = 0.25: w = 0
For i = 1 To 25000
lc = mgr * l / (mgr + mt): M = (mgr + mt) * g * lc * Sin(fi)
dd = R * R + l * l - 2 * R * l * Cos(fi):
Iner = mgr * dd + 2 * mt * R * R: w = w - M / Iner * dt
fi = fi + w * dt: x = x - R * w * dt: t = t + dt
If i Mod 100 = 0 Then
Cells(i / 100, 1) = t: Cells(i / 100, 2) = fi: End If
Next
End Sub

7.12. Кольцо со смещенным центром масс колеблется на горизонтальной плоскости. Используя решение предыдущей задачи, постройте
график зависимости периода колебаний от амплитуды.
Компьютерная модель (программа 7.2) позволяет провести серию вычислительных экспериментов, определяя период колебаний
кольца T при различных амплитудах A (начальных угловых координатах

ϕ z 0 ). Для нахождения T необходимо по таблице найти время, в

течение которого система из крайнего положения переходит в положение равновесия, и умножить его на 4. Получающийся график зависимости периода колебаний от амплитуды представлен на рис. 7.4.2.
7.13. Шар радиусом 1 м и плотностью ρ = 800 кг/м 3 плавает вблизи
поверхности жидкости плотностью ρ 0 = 1000 кг/м 3 . Определите глубину
погружения центра шара.

Рис. 7.6. Определение глубины погружения шара в жидкость

58

Майер Р. В. Решение физических задач в электронных таблицах Excel

Шар погрузится на такую глубину h , при которой сила тяжести
FT = mg = ρVg уравновешивается силой Архимеда: ρ 0 gVпчт = ρ Vg ,
где V = 4π ⋅ R 3 / 3 – полный объем шара, Vпчт – объем его части, погруженной в жидкость. Чтобы рассчитать Vпчт , разрежем шар горизонтальными плоскостями на диски радиусом ri = ( R 2 − yi2 )0,5 , толщиной ∆y и объемом ∆Vi = Si ∆y = π ⋅ ri2 ∆y (рис. 7.6.1). Суммарный объем:
N

Vпчт = ∑ π ( R 2 − yi2 )∆y .
i =1

Используемая программа 7.3 рассчитывает силу тяжести и силу
Архимеда, действующую на шар, при различных значениях глубины
погружения h ; результат выводится в виде таблицы. На ее основе
строят график FA = FA (h) и находят точку пересечения с прямой

F = mg (рис. 7.6.2). Точное значение h можно определить, уменьшив
шаг ∆y и найдя в таблице строку с примерно равными FA и mg .
Программа 7.3
Private Sub CommandButton1_Click()
rho = 800: rho0 = 1000: R = 1: dy = 0.01: y = -R
While y < R
y = y + dy: rr = R * R - y * y
V = V + 3.14 * rr * dy
FA = rho0 * 9.8 * V: i = i + 1
VV = 4 * 3.14 * R * R * R / 3
Cells(i, 1) = y: Cells(i, 2) = FA
Cells(i, 3) = rho * VV * 9.8
Wend
End Sub

7.14. Конус с углом при вершине 60 0 , радиусом основания 10 см и
плотностью ρ плавает на поверхности жидкости плотностью ρ 0 вершиной вниз. Определите глубину погружения вершины конуса.

59

Майер Р. В. Решение физических задач в электронных таблицах Excel
7. Динамика твердого тела

Содержание

9. Динамический хаос

8. АНАЛИТИЧЕСКАЯ МЕХАНИКА
8.1. Система состоит из N = 10 материальных точек, которые имеют известные массы mi и координаты xi , yi , zi . Необходимо рассчитать
положение центра масс.
Координаты центра масс системы вычисляются так:

r
1 N r
1 N
1 N
1 N
rc = ∑ mi ri , xc = ∑ mi xi , yc = ∑ mi yi , zc = ∑ mi zi .
m i =1
m i =1
m i =1
m i =1
Используется программа 8.1. Она считывает исходные данные mi , xi ,

yi , zi из таблицы, вычисляет координаты центра масс и результаты
выводит в ячейки C13, D13 и E13 (рис. 8.1).

Рис. 8.1. Расчет центра масс системы из 10 частиц
Программа 8.1
Private Sub CommandButton1_Click()
Dim m(10) As Single: Dim x(10) As Single
Dim y(10) As Single: Dim z(10) As Single
For i = 1 To 10
m(i) = Cells(i + 1, 2): x(i) = Cells(i + 1, 3)
y(i) = Cells(i + 1, 4): z(i) = Cells(i + 1, 5): Next
For i = 1 To 10
mm = mm + m(i): Next
For i = 1 To 10
xc = xc + x(i) / mm: yc = yc + y(i) / mm: zc = zc + z(i) / mm
Next
Cells(13, 3) = xc: Cells(13, 4) = yc: Cells(13, 5) = zc
End Sub

60

Майер Р. В. Решение физических задач в электронных таблицах Excel

8.2. Кривошип ОА радиусом R = 12 см равномерно вращается со
скоростью ω = 1,5 рад/с вокруг точки O (рис. 8.2). К точке А прикреплен
шатун АВ длиной L = 0,16 см, другой конец которого соединен с ползуном, колеблющимся вдоль оси Ox . Изучите зависимость координаты и
скорости ползуна от времени и найдите их значения в момент τ ' = 2,7 с.

Рис. 8.2. Движение шатуна и ползуна

Рис. 8.3. Результаты расчетов движения шатуна
Для решения задачи используются формулы:

ϕ = ω ⋅τ , x = R cos ϕ + L2 − R 2 sin 2 ϕ , υ x = ( x2 − x1 ) / ∆τ .
Скорость ползуна может быть найдена путем взятия производной
υ x = x't (τ ) численным методом. Необходимо создать таблицу, в которой вычислялись бы значения ϕ , x и υ в моменты τ = i∆τ , где

∆τ = 0,02 с. На ее основе строятся графики x(τ ) , υ x (τ ) (рис. 8.3).
8.3. В эпициклическом механизме (рис. 8.4.1) центр неподвижной
шестерни радиуса r1 = 12 см соединен кривошипом с шестернейсателлитом радиусом r2 = 5 см, находящейся во внешнем зацеплении,

61

Майер Р. В. Решение физических задач в электронных таблицах Excel

на которой установлен диск радиусом L =| O ' A |= 3 см. Рассчитайте траекторию движения точки A на краю диска. Найдите ее координаты и скорость в момент τ ' = 1,6 с.
Для расчета координат точки A используются формулы:

ϕ1 = ωτ , ϕ 2 = ω (r1 + r2 )τ / r2 , x' = (r1 + r2 ) cos ϕ1 ,
y ' = (r1 + r2 ) sin ϕ1 , x A = x'+ L cos ϕ 2 , y A = y '+ L sin ϕ 2 .
Здесь x ' и y ' – координаты точки O ' . Для определения скорости точки A необходимо вычислить производные от x A и y A по времени:
2
2
υ Ax = ( x tA+1 − x tA ) / ∆τ , υ Ay = ( y tA+1 − y tA ) / ∆τ , υ A = υ Ax
+ υ Ay
.

Используется программа 8.2. На рис 8.5.1, 8.5.2 и 8.5.3 представлены
результаты расчетов траектории движения точки A при следующих
параметрах: 1) r1 = 12 см, r2 = 5 см, L = 3 см; 2) r1 = 17 см, r2 = 3 см,

L = 4 см; 3) r1 = 6 см, r2 = 8 см, L = 4 см.

Рис. 8.4. К решению задач 8.3 и 8.4
Программа 8.2
Private Sub CommandButton1_Click()
dt = 0.1: r1 = 12: r2 = 5: L = 3: w = 0.4
While t < 30
t = t + dt: i = i + 1
f1 = w * t: f2 = w * (r1 + r2) * t / r2
x1 = (r1 + r2) * Cos(f1): y1 = (r1 + r2) * Sin(f1)
x2 = x1 + L * Cos(f2): y2 = y1 + L * Sin(f2)
Cells(i, 1) = x1: Cells(i, 2) = y1
Cells(i, 3) = x2: Cells(i, 4) = y2
Wend
End Sub

62

Майер Р. В. Решение физических задач в электронных таблицах Excel

Рис. 8.5. Траектории движения точки A при различных r1 , r2 , L
8.4. Решите предыдущую задачу для случая, когда неподвижная
шестерня и шестерня-сателлит находятся во внутреннем зацеплении
(рис. 8.4.2).
8.5. Камень массой m = 1 кг брошен вертикально вверх со скоростью υ = 10 м/c так, что в момент τ1 = 0 его координата y1 = 0 , а в момент τ 2 = 1 c она равна y2 = 5 м. Численным методом найдите действие
S для действительного пути. Варьируя коэффициенты в уравнении
движения y = C1τ 2 / 2 + C2τ + C3 , определите действие S ' при движении
по окольным путям и убедитесь в том, что оно всегда больше.
Камень движется с ускорением g = 10 м/с 2 , направленным вниз,
ось Oy направим вверх. В начальный момент времени τ1 = 0 камень
имеет направленную вверх скорость 10 м/с и координату y (0) = 0 .
Дважды интегрируя уравнение &y& = − g = C1 , получаем:

y& = C1τ + C2 ,

y = C1τ 2 / 2 + C2τ + C3 . Постоянные интегрирования определяются из
условий: y& (0) = 10 м/с, y (0) = 0 , y (1) = 5 м. Отсюда следует, что

C1 = −10 м/с 2 , C2 = 10 м/с, C3 = 0 , то есть точка движется по закону
y& = −10τ + 10 , y = −10τ 2 / 2 + 10τ .
Будем изменять функцию y = y (τ ) так, чтобы в моменты времени τ 1 = 0 и τ 2 = 1 c система находилась в положениях y = 0 и y = 5 м
соответственно. То есть должны выполняться условия

C1τ 12 / 2 + C2τ 1 = 0 ,

C1τ 22 / 2 + C2τ 2 = 5 .

Так как τ1 = 0 , то первое уравнение приводит к тождеству. Подставляя во второе уравнение τ 2 = 1 с, получаем уравнение, связывающее

63

Майер Р. В. Решение физических задач в электронных таблицах Excel

коэффициенты C1 и C 2 : C 2 = 5 − C1 / 2 . Если в цикле с некоторым шагом изменять C 2 , то можно, вычисляя C1 , построить графики

y ′ = y ′(τ ) , соответствующие окольным путям движения системы при
фиксированных y (0) = 0 и y (1) = 5 м (рис. 8.6.1).
Программа 8.3
Private Sub CommandButton1_Click()
dt = 0.00001: m = 1: g = 10: v = 10
y = 0: S = 0: C1 = Cells(2, 1)
C2 = 5 - C1 / 2: Cells(2, 2) = C2
While t < 1
t = t + dt
y = C1 * t * t / 2 + C2 * t
v = C1 * t + C2
S = S + (m * v * v / 2 - m * g * y) * dt
Wend
Cells(2, 5) = S
End Sub

Рис. 8.6. Действительный и окольные пути. График S = S (C1 ) .
Программа 8.3 позволяет, исходя из значения C1 , определить C2
и вычислить действие S по формуле:

64

Майер Р. В. Решение физических задач в электронных таблицах Excel
τ2

S=

∫
τ

1

∑ (my&
N

L( y, y& , t )dt ≈

2

)

/ 2 − mgy ∆τ ,

i =1

где L( y , y& , t ) = my& 2 / 2 − mgy – функция Лагранжа системы. График зависимости S = S (C1 ) действия от коэффициента C1 показан на
рис. 8.6.2. Видно, что минимального значения –16,667 Дж ⋅ с действие
достигает при C1 = −10 м/с 2 и C2 = 10 м/с, что соответствует движению по действительному пути. При движении системы по окольным путям значение S больше. Малые вариации действительного
пути приводят к изменениям действия на величину более высокого
порядка малости, как того требует принцип Гамильтона.
8.6. Известны главные центральные моменты инерции тела

I1 = 5 кг ⋅ м 2 , I 2 = 7 кг ⋅ м 2 , I 3 = 3 кг ⋅ м 2 . Рассчитайте момент инерции
относительно оси OA, если заданы координаты точки A ( x A , y A , z A ). Оси
координат повернули на угол ϕ вокруг оси Oz . Запишите тензор инерции и снова рассчитайте момент инерции тела относительно оси OA.

Рис. 8.7. К задаче на расчет тензора инерции
Пусть оси координат x , y и z совпадают с главными осями
инерции тела (рис. 8.7). Осевые моменты инерции I x = I1 , I y = I 2 ,

I z = I 3 , а центробежные – I xy = I xz = I yz = 0 . Направляющие косинусы
| OA |= ( x A2 + y A2 + z 2A ) 0,5 , cos α = x A / | OA | , cos β = y A / | OA | ,
cos γ = z A / | OA | . Зная главные центральные моменты инерции тела,

оси OA:

можно рассчитать момент инерции тела относительно оси OA:
I OA = I1 cos 2 α + I 2 cos 2 β + I 3 cos 2 γ .

65

Майер Р. В. Решение физических задач в электронных таблицах Excel

При x A = 1,0 м, y A = 0,8 м, z A = 0,6 м получающееся значение момента
инерции 5,28 кг ⋅ м 2 .
Допустим, система координат повернулась вокруг оси Oz на
угол ϕ = 400 = 0,698 рад; при этом взаимное расположение тела, его
эллипсоида инерции и оси OA остались неизменными. В новой системе
координат отрезок OA повернулся на − 0,698 радиан. Матрица поворота вокруг оси Oz имеет вид:

 cos ϕ

M z (ϕ ) =  − sin ϕ
 0


sin ϕ
cos ϕ
0

0

0 .
1 

Рис. 8.8. Расчет тензора инерции до и после поворота осей координат
Новые значения координат точки A вычисляются по формулам:
x' A = x A cos ϕ + y A sin ϕ , y ' A = − x A sin ϕ + y A cos ϕ , z ' A = z A .
Направляющие косинусы оси OA в новой системе координат:
cos α ' = x' A / | OA | ,
cos β ' = y ' A / | OA | ,

cos γ ' = z ' A / | OA | , | OA |= ( x'2A + y '2A + z '2A ) 0,5 .
Элементы тензора инерции в новой системе координат Ox' y ' z ' :

I x ' x ' = I1 cos 2 ϕ + I 2 cos 2 (π / 2 − ϕ ) + I 3 cos 2 (π / 2) = I1 cos 2 ϕ + I 2 sin 2 ϕ ,

66

Майер Р. В. Решение физических задач в электронных таблицах Excel

I y ' y ' = I1 cos 2 (π / 2 − ϕ ) + I 2 cos 2 ϕ + I 3 cos 2 (π / 2) = I1 sin 2 ϕ + I 2 cos 2 ϕ ,

I z ' z ' = I 3 , I x ' y ' = I y ' x ' = ∫ x' y ' dm = ( I 2 − I1 ) sin ϕ cos ϕ , I x ' z ' = I y ' z ' = 0 .
Исходя из тензора инерции в новой системе координат Ox' y ' z ' , можно
рассчитать момент инерции тела относительно оси OA по формуле:

I OA = I x ' x ' cos 2 α '+ I y ' y ' cos 2 β '+ I z ' z ' cos 2 γ '−

− 2 cos α ' cos β '⋅I x ' y ' − 2 cos α ' cos γ '⋅I x ' z ' − 2 cos β ' cos γ '⋅I y ' z ' .
Используется документ Excel, в котором перечисленные выше величины вычисляются по формулам в ячейках таблицы. В результате
получается такое же значение момента инерции 5,28 кг ⋅ м 2 , как и в
предыдущем случае (рис. 8.8).
8.7. По внутренней поверхности вращающейся трубы радиуса R
катится цилиндр радиуса r , ось которого соединена со стержнем OA
(рис. 8.9). Скорость вращения трубы изменяется по заданному закону
ω (τ ) . Между поверхностью трубы и цилиндром, а также в подшипниках
O и A действует сила вязкого трения. Рассчитайте движение системы.
Система имеет две степени свободы. Введем две обобщенные
координаты q1 и q2 , измеряемые в радианах. Учтем, что: 1) сила вязкого трения между поверхностями трубы и цилиндра пропорциональна скорости движения точек одной поверхности относительно другой; 2) тормозящий момент, действующий в подшипниках, пропорционален угловой скорости вращения.

Рис. 8.9. Движение цилиндра внутри вращающейся трубы

67

Майер Р. В. Решение физических задач в электронных таблицах Excel

Рис. 8.10. Результаты расчетов движения цилиндра
Запишем систему уравнений Лагранжа:

d ∂L ∂L
−
= Qkтр ,
dt ∂q& k ∂qk
L = T −U =

k = 1, 2 ,

m 2
I
q&1 ( R − r ) 2 + q& 22 + mg ( R − r ) cos q1 ,
2
2

mq&&1 ( R − r ) 2 + mg ( R − r ) sin q1 = FТР R − rq&1 , Iq&&2 = FТР r − rq& 2 .
Здесь FТР = rS (ω ⋅ R − q&1R − q& 2 r ) – сила вязкого трения между внутренней поверхностью трубы и цилиндром. Пусть угловая скорость вращения трубы сначала равна 0, при τ > 1 с увеличивается от 0 до 20
рад/с, затем остается постоянной, после чего при τ = 4 с резко обращается в 0. Используется программа 8.4. На рис 8.10.1 представлен
график ω1 (t ) = q&1 (t ) ; на рис. 8.10.2 – ω2 (t ) = q&2 (t ) для случая, когда начальная координата q1 = 1,2 рад, q2 = 0. Из графиков следует, что когда труба покоится (τ < 1 c и τ > 4 с), цилиндр совершает затухающие
колебания. В момент τ = 3,6 c цилиндр вращается практически с постоянной угловой скоростью ω2 .

68

Майер Р. В. Решение физических задач в электронных таблицах Excel

Программа 8.4
Private Sub CommandButton1_Click()
dt = 0.0001: X0 = 320: Y0 = 240: R = 180: rc = 20: w = 0
m = 0.8: rs = 3: rr = 1: Iner = 0.5: q1 = 1.2: i = 1
While t < 5
t = t + dt: i = i + 1
If t > 1 Then w = 10 * (t - 1): If w > 20 Then w = 20
If t > 4 Then
w = 0: Ftr = rs * (w * R - w1 * R - w2 * rc)
End If
e1 = (Ftr * R - rr * w1 –
m * 980 * (R - rc) * Sin(q1))/ m / Sqr(R - rc)
e2 = (Ftr * rc - rr * w2) / Iner
w1 = w1 + e1 * dt: w2 = w2 + e2 * dt
q1 = q1 + w1 * dt: q2 = q2 + w2 * dt: fi = fi + w * dt
If i Mod 100 = 1 Then Cells(i / 100, 1) = t
If i Mod 100 = 1 Then Cells(i / 100, 2) = w * 20
If i Mod 100 = 1 Then Cells(i / 100, 3) = w1
If i Mod 100 = 1 Then Cells(i / 100, 4) = w2
If i Mod 100 = 1 Then Cells(i / 100, 5) = q1
Wend
End Sub

69

Майер Р. В. Решение физических задач в электронных таблицах Excel
8. Аналитическая механика

Содержание

10. Молекулярная физика

9. ДИНАМИЧЕСКИЙ ХАОС
9.1. Биллиард Синая состоит из горизонтально расположенного
ящика с изогнутыми стенками, в котором установлены один или несколько вертикальных стержней. По дну ящика движется шарик, который упруго соударяясь со стенками и стержнями и отскакивает от них. Изучите
движение шарика при различных начальных координатах и скоростях,
убедитесь в том, что его траектория очень чувствительна к начальным
условиям.
Рассмотрим движение шарика в биллиарде прямоугольной формы с наклонными стенками и конусообразным выступом, вершина A
которого имеет координаты x1 и y1 . При приближении шарика к точке A на расстояние l , меньшее r , на него действует постоянная по

r
F , направленная от A. Ее проекции равны:
Fx = F cos α = F ( x − x1 ) / l , Fy = F sin α = F ( y − y1 ) / l . При столкновении с

модулю

сила

наклонной стенкой на шарик действует постоянная сила, направленная перпендикулярно стенке. С помощью программы 9.1 можно рассчитать траекторию движения шарика и убедиться в том, что небольшие вариации начальных условий приводят к тому, что траектория движения шарика сильно изменяется (рис. 9.1). Чтобы промоделировать ситуацию, когда на дне биллиарда имеется конусообразная впадина, необходимо изменить знаки в выражениях для Fx и Fy .

Рис. 9.1. Хаотическое движение шарика в биллиарде Синая.

70

Майер Р. В. Решение физических задач в электронных таблицах Excel

Программа 9.1
Private Sub CommandButton1_Click()
m = 0.001: x1 = 130: y1 = 150: x = 10: y = 100
vx = 1: vy = 10: dt = 0.01
While t < 80
i = i + 1: t = t + dt: Fx = 0: Fy = 0
S = Sqr((x1 - x) * (x1 - x) + (y1 - y) * (y1 - y))
If S < 40 Then F = 0.04
If x < 80 Then Fx = 0.1: If x > 220 Then Fx = -0.1
If y < 80 Then Fy = 0.1: If y > 220 Then Fy = -0.1
Fx = Fx + F * (x – x1) / S: Fy = Fy + F * (y – y1) / S
ax = Fx / m: ay = Fy / m: vx = vx + ax * dt: vy = vy + ay * dt
x = x + vx * dt: y = y + vy * dt
If i Mod 10 = 0 Then
Cells(i / 10, 1) = x: Cells(i / 10, 2) = y
End If
Wend
End Sub

9.2. Размножение насекомых на изолированном острове описывается квадратичным отображением x t +1 = ax t (1 − x t ) . Изучите поведение
этой функции в зависимости от параметра a , изменяющегося в интервале от 0 до 4, при t > 200. Постройте график x = x (a ) и найдите точки
бифуркации, в которых происходит расщепление пути эволюции этой
системы. Определите a , при которых наступает динамический хаос.
С одной стороны, численность насекомых x t +1 в ( t + 1)-м поколении пропорциональна числу насекомых xt в предыдущем t -м поколении. С другой стороны, с ростом x t насекомые начинают испытывать нехватку энергетического ресурса (пищи) и их коэффициент
размножения уменьшается. Когда численность насекомых xt достигает своего предела A , коэффициент размножения уменьшается до
0. Получаем уравнение: x t +1 = rx t ( A − x t ) .

Рис. 9.2. Квадратичное отображение x

71

t +1

= axt (1 − xt )

Майер Р. В. Решение физических задач в электронных таблицах Excel

Используется программа 9.2, которая вычисляет по 100–1000
значений xi при t > 200 и различных значениях параметра a из интервала [0; 4] и выводит результаты на экран в графическом виде. Результат ее работы представлен на рис. 9.2. Видно, что пока a < 3,
каждому значению параметра a соответствует одно значение x .
Дальнейшее увеличение a приводит к каскаду бифуркаций удвоения
периода: сначала одному значению a отвечают два значения x , затем 4, затем 8 и т. д. Последовательность бифуркаций, начинаясь
при a = 3, переходит в хаос при a > 3,6, которому соответствует уже
не непрерывная линия, а фрактальное множество точек. Полосы хаоса перемежаются с "окнами периодичности".
Программа 9.2
Private Sub CommandButton1_Click()
x = 0.001: a = 0.01: j = 1: b = 200
While a < 3.9
a = a + 0.01: x = 0.001
For i = 0 To 250
x = x * a * (1 - x)
If i > b Then j = j + 1: Cells(j, 1) = a: Cells(j, 2) = x
Next i
Wend
End Sub

9.3. Маятник Дафинга представляет собой шарик, находящийся
внутри потенциальной ямы с двумя углублениями так, что его потенциальная энергия U ( x) = k ( x 4 / 4 − x 2 / 2) . Изучите хаотическиеколебания
шарика в случае, когда на него действует периодически изменяющаяся
сила Fx (τ ) = Fm cos(ω ⋅ τ ) .
На шарик внутри потенциальной ямы действует возвращающая
сила: F ' x = −∂U / ∂x = − k ( x 3 − x ) . Вынужденные колебания маятника
Дафинга описываются уравнением: m&x& + rx& + k ( x − x) = Fm cos(ω ⋅ τ ) .
3

Для его решения используется программа 9.3.
Программа 9.3
Private Sub CommandButton1_Click()
F = 1: k = 4: r = 0.5: m = 1: w = 2.3: dt = 0.001
fi = 2.64: faza = 1
While t < 80
j = j + 1: t = t + dt: z = Cos(w * t + faza)
a = (F * Cos(w * t) - k * (x * x * x - x) - r * v) / m

72

Майер Р. В. Решение физических задач в электронных таблицах Excel

v = v + a * dt: x = x + v * dt
If j Mod 20 = 0 Then
Cells(j / 20, 1) = t: Cells(j / 20, 2) = x
Cells(j / 20, 3) = v
End If
Wend
End Sub

Рис. 9.3. Маятник Дафинга: графики x(t ) , υ x (t ) и фазовая кривая
Получающиеся графики x(τ ) и

υ x (τ ) и фазовая кривая представ-

лены на рис. 9.3. Видно, что при заданных параметрах системы, амплитуде и частоте вынуждающей силы происходят хаотические колебания. Шарик, совершив несколько колебаний в одном углублении,
перескакивает в другое. Его движение неустойчивое, сильно зависит
от начальных условий; фазовая кривая имеет многочисленные самопересечения.
9.4. Для маятника Дафинга, который совершает вынужденные колебания, постройте сечение Пуанкаре при различных фазах ϕ вынуждающей силы Fx (τ ) = Fm cos(ω ⋅ τ + ϕ 0 ) .
Необходимо решить уравнение m&x& + rx& + k ( x − x) = Fm cos(ω ⋅ t )
3

тем же методом, как и в задаче 9.3. Состояние системы можно охарактеризовать координатой x , проекцией импульса шарика p x и проекцией силы Fx (τ ) . Сечение Пуанкаре, соответствующее фазе

ϕ вы-

нуждающей силы, строится так: в момент, когда cos(ω ⋅ τ + ϕ ) обращается в 0, на фазовой плоскости ставится точка с координатами x
и p x . Используется программа 9.4, результаты моделирования для
различных ϕ приведены на рис. 9.4. Видно, что сечение Пуанкаре имеет фрактальную структуру, что характерно для хаотических коле73

Майер Р. В. Решение физических задач в электронных таблицах Excel

баний. В программе время τ изменяется дискретно с шагом ∆τ , поэтому условие cos(ω ⋅ τ + ϕ ) = 0 практически никогда не выполняется.
Вместо него в программе проверяется другое условие, когда одновременно cos(ω ⋅ (τ − ∆τ ) + ϕ ) < 0 и cos(ω ⋅ τ + ϕ ) > 0 .

Рис. 9.4. Сечение Пуанкаре для маятника Дафинга
при фазах ϕ = 0, 1, 2, 4 рад
Программа 9.4
Private Sub CommandButton1_Click()
F = 1: k = 4: r = 0.5: m = 1: w = 2.3: dt = 0.001
fi = 2.64: faza = Cells(1, 5)
While t < 5000
t = t + dt: z = Cos(w * t + faza)
a = (F * Cos(w * t) - k * (x * x * x - x) - r * v) / m
v = v + a * dt : x = x + v * dt
If (z > 0) And (zz < 0) Then
j = j + 1: Cells(j, 1) = x: Cells(j, 2) = v
End If
zz = z
Wend
End Sub

9.5. Маятник Дафинга совершает свободные колебания. Промоделируйте перемешиваемость фазового объема в случаях, когда коэффициент затухания: 1) равен 0; 2) отличен от нуля.
Представим себе совокупность из 1600 одинаковых маятников
Дафинга, совершающих свободные колебания, которые отличаются
только начальными условиями x0 , p0 . Движение каждого маятника
описывается уравнением: m&x& + rυ + k ( x − x) = 0 . Пусть в момент
3

τ = 0 фазовые точки, характеризующие начальное состояние маятников, находятся внутри прямоугольника, ограниченного интервалами [ x, x + ∆x ] и [ p, p + ∆p ] . Программа 9.5 рассчитывает состояние
каждого маятника в заданный момент времени τ ' , который заранее
вводят в ячейку E1, и записывает результаты вычислений в столбцы

74

Майер Р. В. Решение физических задач в электронных таблицах Excel

А и B. Получается таблица из 1600 строк, на ее основе строится фазовый портрет данного ансамбля маятников в момент τ ' .
Программа 9.5
Private Sub CommandButton1_Click()
m = 1: k = 1: r = 0.02: dt = 0.002
For i = 1 To 40: For j = 1 To 40
x = 0.04 * i: v = 0.04 * j
t = 0: s = s + 1
Vremya = Cells(1, 5)
While t < Vremya
t = t + dt
a = (-k * (x * x * x - x) - r * v) / m
v = v + a * dt: x = x + v * dt
Wend: Cells(s, 1) = x: Cells(s, 2) = v
Next: Next
End Sub

На рис. 9.5 представлены результаты вычислений для моментов
τ ' = 10, 20, 60. Видно, что при τ ' → ∞ происходит расползание фазового
объема, его перемешивание в фазовом пространстве, что свидетельствует о хаотичности колебаний. Если колебания затухают ( r > 0), то фазовый объем через некоторое время уменьшается до 0. В случае незатухающих колебаний ( r = 0) фазовый объем ведет себя как несжимаемая жидкость (то есть сохраняет свою величину), как того требует теорема Лиувилля.

Рис. 9.5. Перемешивание фазового объема для маятника Дафинга
в моменты τ ' = 10, 20 и 60
9.6. Промоделируйте странный аттрактор Лоренца, уравнения которого имеют вид: x& = a ( − x + y ) , y& = bx − y − xz , z& = −cz + xy .
Эта система дифференциальных уравнений описывает движение конвекционных потоков в атмосфере. Конечно-разностные уравнения выглядят так:

υ xt +1 = a (− x t + y t ) , υ ty+1 = bx t − y t − x t z t , υ zt +1 = −cz t + x t y t ,

75

Майер Р. В. Решение физических задач в электронных таблицах Excel

x t +1 = x t + υ xt +1∆τ , y t +1 = y t + υ ty+1∆τ , z t +1 = z t + υ zt +1∆τ .
Используется программа 9.6. Получающиеся фазовая кривая и график
хаотических колебаний представлены на рис. 9.6. Система эволюционирует так, что фазовая кривая описывает закрученную восьмерку,
захватывающую две большие области фазового пространства, образованного осями x , y и z . Аттрактор Лоренца называется странным потому, что его фазовая кривая имеет фрактальную структуру.
Программа 9.6
Private Sub CommandButton1_Click()
a = 15: b = 30: c = 2.2: dt = 0.001
x = 20: y = 20: z = 5
While i < 20000
i = i + 1
x = x + a * (-x + y) * dt
y = y + (b * x - y - x * z) * dt
z = z + (-c * z + x * y) * dt
If i Mod 10 = 0 Then Cells(i / 10, 1)
If i Mod 10 = 0 Then Cells(i / 10, 2)
If i Mod 10 = 0 Then Cells(i / 10, 3)
If i Mod 10 = 0 Then Cells(i / 10, 4)
Wend
End Sub

=
=
=
=

i * dt
x
y
z

Рис. 9.6. Аттрактор Лоренца: хаотическое движение
9.7. Исследуйте странный аттрактор Ресслера, задаваемый системой уравнений: x& = −( y + z ) , y& = x + y / 5 , z& = 1 / 5 + z ( x − µ ) .
Программа 9.7, решающая эту систему дифференциальных
уравнений, состоит из цикла по времени, в котором рассчитываются
значения x , y , z в следующий момент τ + ∆τ :

υ x = −( y + z ) , υ y = x + y / 5 , υ z = 1 / 5 + z ( x − µ ) ,

76

Майер Р. В. Решение физических задач в электронных таблицах Excel

x t +1 = x t + υ xt +1∆τ , y t +1 = y t + υ ty+1∆τ , z t +1 = z t + υ zt +1∆τ .
Результаты вычислений записываются в ячейки таблицы, на их основе могут быть построены графики x(τ ) , y (τ ) , z (τ ) , а также фазовая
кривая, которая является фракталом (рис. 9.7). Программа позволяет
изучить странный аттрактор Ресслера при различных значениях параметра µ .
Программа 9.7
Private Sub CommandButton1_Click()
a = 0.9: b = 0.3: c = -0.3: dt = 0.005
x = 20: y = 10: z = 5: m = 20
While i < 20000
i = i + 1
x = x - (y + z) * dt
y = y + (x + y / 5) * dt
z = z + (1 / 5 + z * (x - m)) * dt
If i Mod 10 = 0 Then Cells(i / 10, 1) =
If i Mod 10 = 0 Then Cells(i / 10, 2) =
If i Mod 10 = 0 Then Cells(i / 10, 3) =
If i Mod 10 = 0 Then Cells(i / 10, 4) =
Wend
End Sub

i * dt
x
y
-(y + z)

Рис. 9.7. Моделирование аттрактора Ресслера

77

Майер Р. В. Решение физических задач в электронных таблицах Excel
9. Динамический хаос

Содержание

11. Термодинамика

10. МОЛЕКУЛЯРНАЯ ФИЗИКА
10.1. Промоделируйте хаотическое движение броуновской частицы, которая вследствие столкновений с молекулами случайно изменяет
направление и скорость движения.
Используется программа 10.1. Она содержит цикл, в котором
случайно выбирается угол между направлением скорости и осью Ox , а
также смещение L за время ∆τ = 0,05. При этом рассчитываются
координаты броуновской частицы в следующий момент τ + ∆τ . Результаты выводятся в таблицу, на их основе строится траектория
частицы (следует выбрать опцию “Точечная диаграмма, на которой
точки соединены отрезками”). Траектория является фракталом.
Программа 10.1
Private Sub CommandButton1_Click()
While t < 20
t = t + 0.05: i = i + 1
a = 2 * 3.14 * Rnd: L = Rnd
x = x + L * Cos(a): y = y + L * Sin(a)
Cells(i, 1) = x: Cells(i, 2) = y
Wend
End Sub

Рис. 10.1. Моделирование броуновского движения
10.2. Молекула из-за хаотических столкновений с другими молекулами случайным образом движется вдоль оси Ox . При этом ее координата через равные промежутки времени ∆τ принимает значения x , x ,
1

2

x 3 , …, x N , где N – число шагов. Изучите распределение конечной координаты x

N

блуждающей точки при различном числе шагов N .

78

Майер Р. В. Решение физических задач в электронных таблицах Excel

Рассмотрим задачу Пирсона о случайных блужданиях: необходимо определить смещение пешехода, сделавшего N шагов равной длины в случайных направлениях. Применим метод статистических испытаний. Используемая программа 10.2, содержит цикл, в котором
случайным образом определяется направление и величина смещения
молекулы на следующем временном шаге и осуществляется N итераций, в результате чего определяется конечная координата молекулы x N . За один шаг точка смещается на случайную величину из интервала [ -0,5; 0,5 ] , для этого применяется оператор: x = x + Rnd 0.5. Программа воспроизводит N isp = 1000–10000 реализаций N случайных смещений молекулы из начала координат, вычисляет дисперсию конечной координаты x N , подсчитывает число n j попаданий
частицы в интервалы [ j; j + 1] ( j = -15, -14, ..., 14, 15) и выводит их в
таблицу.
Программа 10.2
Private Sub CommandButton1_Click()
Dim n(30): Kol_shag = Cells(1, 7): Kol_isp = 5000
For i = 1 To Kol_isp
For j = 1 To Kol_shag
x = x + Rnd - 0.5: Next
Cells(i, 1) = x: S = x * x + S
For j = 1 To 30
If (x >= j - 15) And (x < j - 14) Then n(j) = n(j) + 1
Next: x = 0: Next
Cells(2, 7) = S / Kol_isp
For j = 1 To 30
Cells(j, 2) = j - 15: Cells(j, 3) = n(j)
Next
End Sub

На рис. 10.2 показаны значения конечной координаты x N после

N = 300 шагов и распределение случайной величины x N , полученное
после 10000 испытаний. Видно, что оно похоже на нормальное распределение Гаусса. Чтобы определить вероятность p j попадания
частицы в интервал [ j; j + 1] ( j = − 15, − 14,..., 14, 15 ), необходимо n j
разделить на количество испытаний N isp . Повторите моделирование
при N = 100 или 500. Убедитесь, что с ростом числа шагов N об-

79

Майер Р. В. Решение физических задач в электронных таблицах Excel

ласть пространства,
увеличивается.

в

которой

может

оказаться

частица,

Рис. 10.2. Моделирование случайного блуждания частицы
10.3. Молекула из-за столкновений с другими молекулами случайным образом перемещается вдоль оси Ox . Рассчитайте средний квадрат ее координаты x N при заданном числе шагов N . Постройте график
зависимости среднего квадрата смещения “блуждающей” частицы

S 2 =< ( x N ) 2 > от N , убедитесь, что это прямая пропорциональность.
Используется программа 10.2, в которой задается число шагов
N = 10, 20, ..., 500 и число испытаний N = 5000 . При этом многократно моделируется случайное одномерное движение частицы из начала
координат O , накапливается сумма квадратов смещений ( x N ) 2 , которая затем делится на N isp . В ячейку G2 записывается значение

S 2 =< ( x N ) 2 > , усредненное по всем испытаниям:
N

1 isp N 2
S =< ( x ) >=
∑ ( xi ) .
N isp i =1
2

N 2

80

Майер Р. В. Решение физических задач в электронных таблицах Excel

Рис. 10.3. Средний квадрат смещения S 2 пропорционален N
Если построить график зависимости среднего квадрата смещения “блуждающей” частицы S 2 =< ( x N ) 2 > от количества шагов N ,
то получится прямая пропорциональность (рис. 10.3). При этом моделируется одномерная диффузия молекул газа, которые сначала находились в небольшом объеме. Результаты можно интерпретировать
так: 1) с течением времени увеличивается область, в которой может
находится хаотически движущаяся молекула; 2) если из небольшого
объема выходит 10000 молекул некоторого вещества, то через некоторое время (пропорциональное количеству шагов N ) молекулы
вследствие диффузии “расползутся” и займут большую область.
10.4. В длинном сосуде находится N молекул газа, которые движутся хаотически. Каждая молекула с равными вероятностями 0,5 оказывается то в правой, то в левой половине сосуда. Постройте график
зависимости вероятности нахождения в левой половине сосуда n молекул, при различных N .
Если в левой половине оказалось n частиц, то в правой – ( N − n)
частиц. Рассчитаем вероятности различных распределений молекул
в левой и правой половинах сосуда и построим график p (n) . Для этого будем случайным образом задавать координаты N молекул и подсчитывать их число n в левой половине сосуда. Многократно (более
1000 раз) повторяя эту процедуру, определим число реализаций n1 ,

n2 , n3 , …, n N −1 , n N , при которых в левой части окажется 1, 2 , 3 , …,

81

Майер Р. В. Решение физических задач в электронных таблицах Excel

N − 1, N молекул. Соответствующие вероятности pi рассчитываются по формуле: pi = n / N isp . Используется программа 10.3, результаты для 30 молекул представлены на рис. 10.4. Максимальную вероятность имеет состояние с наибольшей энтропией, при котором в
левой половине находится N / 2 молекул.
Программа 10.3
Private Sub CommandButton1_Click()
Dim n(50): N_mol = 30: N_isp = 5000
While k < N_isp
k = k + 1: s = 0
For i = 1 To N_mol
xx = Rnd(100): If xx < 0.5 Then s = s + 1
Next: Cells(k, 1) = s
For i = 1 To N_mol
If s = i Then n(i) = n(i) + 1
Next: Wend
For i = 1 To N_mol
Cells(i, 3) = i: Cells(i, 4) = n(i) / N_isp
Next
End Sub

Рис. 10.4. Зависимость вероятности p от числа n
молекул в левой половине сосуда
10.5. Имеется узкий сосуд длиной L , частично заполненный газом.
Сначала молекулы заполняют 1/8 часть сосуда, а затем открывается

82

Майер Р. В. Решение физических задач в электронных таблицах Excel

перегородка, и молекулы, двигаясь хаотически, постепенно заполняют
весь сосуд. Изучите изменение энтропии системы с течением времени.
Вдоль сосуда направим ось Ox . Молекулы газа движутся хаотически, энтропия системы практически не зависит от их координат

yi . Координаты xi в течение малого промежутка времени изменяются на случайные величины ∆x , лежащие в интервале [−2; 2] . Сначала
все молекул N находятся в левой части сосуда, xi < L = 10. В момент
t = 500 перегородка исчезает, и молекулы начинают постепенно заполнять весь сосуд. Для нахождения энтропии разобьем сосуд вертикальными плоскостями на 40 элементарных объемов толщиной 2.
Число молекул в j -м объеме обозначим m j . Энтропия системы:
40

40

S = − ∑ p j ln p j = − ∑

mj

j =1 N

j =1

ln

mj
N

.

Для моделирования этой системы и расчета энтропии в последовательные моменты времени используется программа 10.4. Она создает таблицу, в которой представлены значения энтропии S в различные моменты времени, а также значения m j при j = 1, 2, …, 40. По
этим результатам строят графики m j ( j ) и S (τ ) (рис. 10.5). Видно,
что процесс расширения газа необратимый, энтропия с течением
времени возрастает. График m j ( j ) показывает распределение молекул вдоль оси Ox . При t → ∞ происходит выравнивание концентраций,
нарушаемое флуктуациями xi .
Программа 10.4
Private Sub CommandButton1_Click()
n = 500: Dim x(500): Dim m(40)
For i = 1 To n
x(i) = 10 * Rnd: Next
While t < 3000
t = t + 1: For i = 1 To n
dx = 4 * Rnd - 2: x(i) = x(i) + dx
If x(i) < 0 Then x(i) = 0
If t < 500 Then L = 10 Else L = 80
If x(i) > L Then x(i) = L
Next
If t Mod 50 = 0 Then
enrop = 0: For j = 1 To 40
m(j) = 0.01: For i = 1 To n
If (x(i) > 2 * j - 2) And (x(i) < 2 * j) Then m(j) = m(j) + 1

83

Майер Р. В. Решение физических задач в электронных таблицах Excel

Next: enrop = enrop - m(j) / n * Log(m(j) / n): Next
Cells(t / 50, 1) = t / 50: Cells(t / 50, 2) = enrop: End If
Wend
For j = 1 To 40
Cells(j, 3) = m(j): Next
End Sub

Рис. 10.5. Возрастание энтропии при расширении газа
10.6. Найдите распределение молекул кислорода по скоростям при
температуре 293 К. Молярная масса кислорода M = 32 г/моль. Определите: 1) наиболее вероятную скорость молекул; 2) среднюю квадратическую скорость молекул; 3) доли от общего числа молекул, скорости которых лежат в интервалах от 0 до 500 м/c и от 200 до 600 м/c.
Масса молекулы m = M / N A , где N A = 6,02 ⋅ 10

23

моль

−1

. Распре-

деление молекул газа по скоростям задается уравнением Максвелла:

f (υ ) =

∆N
= Aυ 2 exp(− Bυ 2 ) ,
N ⋅ ∆υ

A = 4π (m / 2πkT )3 / 2 = 4 B 3 / π ,
где k = 1,38 ⋅ 10

− 23

B = m / 2kT ,

Дж/К. Функция распределения f (υ ) показывает от-

носительное число молекул, приходящихся на единицу интервала ско-

84

Майер Р. В. Решение физических задач в электронных таблицах Excel

ростей. Доля молекул, скорость которых меньше заданного значения
υ ′ , определяется интегралом:

∆N (υ < υ ' )
S=
=
N

υ'

n

0

i =1

∫ f (υ )dυ ≈ ∑ f (υi )∆υ ,

где n соответствует υ ′ и равен n = υ ′ / ∆υ . Средняя квадратическая
и наиболее вероятная скорости молекул соответственно равны:

υср.кв. = 3kT / m , υв = 2 RT / M , где R = 8,31 Дж/(моль·К). С помощью
программы 10.5 можно получить таблицу, содержащую столбцы υ ,
f , S , и решить задачу. Получающиеся графики f (υ ) и S (υ ) представлены на рис. 10.6. Из таблицы можно определить, что количества молекул, скорости которых лежат в интервалах от 0 до 500 м/c и
от 200 до 600 м/c, составляют 66 и 72 % от их общего числа.

Рис. 10.6. Распределение молекул по скоростям
Программа 10.5
Private Sub CommandButton1_Click()
m0 = 16 * 0.002 / 6.022E+23
B = m0 / 2 / 1.38E-23 / 293
A = 4 * Sqr(B * B * B / 3.1415): dv = 10
For i = 1 To 100
v = i * dv: f = A * v * v * Exp(-B * v * v)
S = S + f * dv: Cells(i, 1) = v
Cells(i, 2) = f: Cells(i, 3) = S
Next
End Sub

85

Майер Р. В. Решение физических задач в электронных таблицах Excel
10. Молекулярная физика

Содержание

12. Электрическое
и магнитное поле

11. ТЕРМОДИНАМИКА
11.1. В калориметре находится 1 кг льда при –15 0 С. В него опускают нагреватель мощностью P = 250 Вт, который включают на 10 мин.
Как изменяются масса льда и образовавшейся воды, а также температура воды с течением времени?
За время ∆τ выделяется количество теплоты ∆Q = P∆τ , которое идет на нагревание льда, его плавление и нагревание воды.
При нагревании температура увеличивается на ∆T = ∆Q / cm , где c –
удельная теплоемкость вещества, m – его масса. При плавлении
температура остается постоянной, причем за время ∆τ плавится
лед массой ∆m = P∆τ / λ , где λ – удельная теплота плавления. Используется программа 11.1, результаты приведены на рис. 11.1.

Рис. 11.1. Зависимость температуры воды от времени
Программа 11.1
Private Sub CommandButton1_Click()
dt = 0.2: Tem = -20: m = 0.3: c1 = 2100
c2 = 4200: P = 250: Lambda = 330000
While Tem < 0
i = i + 1: t = t + dt: dQ = P * dt
Tem = Tem + dQ / c1 / m
Cells(i, 1) = t: Cells(i, 2) = Tem
Wend: Tem = 0: mL = m
While mL > 0
i = i + 1: t = t + dt: mL = mL - P * dt / Lambda
Cells(i, 1) = t: Cells(i, 2) = Tem
Wend:
While Tem < 24

86

Майер Р. В. Решение физических задач в электронных таблицах Excel

i = i + 1: t = t + dt: dQ = P * dt
Tem = Tem + dQ / c2 / m
Cells(i, 1) = t: Cells(i, 2) = Tem
Wend
End Sub

11.2. Металлическое тело массой 1 кг при температуре T1 = 10 o С
опустили в жидкость массой 2 кг при температуре T2 = 95 o С. Как изменяются средние температуры тел с течением времени? Коэффициент
теплопроводности k = 60 Вт / м ⋅ К , удельные теплоемкости тела и
жидкости равны 450 Дж / кг ⋅ К и 4200 Дж / кг ⋅ К .
В процессе теплообмена за время ∆τ тело отдает жидкости

количество теплоты ∆Q = k (T2 − T1 ) ∆τ . При этом температура воды T2 уменьшается, а температура тела T1 растет:

T1t +1 = T1t + ∆Q / c1m1 , T2t +1 = T2t − ∆Q / c 2 m2 .
Используемая программа 11.2 содержит цикл по времени, в котором
вычисляются значения T1 и T2 в последовательные моменты времени, а результаты выводятся во второй и третий столбцы таблицы.
Затем они используются для построения графиков T1 (τ ) и T2 ( x) (рис.
11.2). Видно, что разность температур с течением времени уменьшается до нуля.

Рис. 11.2. Моделирование теплообмена между двумя телами
Программа 11.2
Private Sub CommandButton1_Click()
k = 60: dt = 0.01: Temp1 = 10: Temp2 = 95
m1 = 1: m2 = 2: c1 = 4200: c2 = 450

87

Майер Р. В. Решение физических задач в электронных таблицах Excel

While t < 30
t = t + dt: dQ = k * (Temp2 - Temp1) * dt
Temp1 = Temp1 + dQ / c1 / m1
Temp2 = Temp2 - dQ / c2 / m2
i = i + 1: Cells(i, 1) = t
Cells(i, 2) = Temp1: Cells(i, 3) = Temp2
Wend
End Sub

11.3. Газ, взятый в количестве 2 моль, занимает объем 1 л при
температуре 320 К. Сначала газ изобарно расширяется до объема 2,5 л,
а затем расширяется изотермически до объема 7 л. Определите давление и температуру в конце процесса, а также совершенную газом работу. Постройте графики в осях p и V ; V и T ; p и T .
Для идеального газа p

t +1

= D /(V t +1 ) γ , pV / T = νR = B . Рассчитав

B = νR , можно определить начальное давление p0 = BT / V . Используемая программа 11.3 содержит два цикла. В первом цикле моделируется изобарное расширение газа до объема 2,5 л, вычисляется объем,
температура и работа:

V t +1 = V t + ∆V ,

T t +1 = p0V t +1 / B ,

At +1 = At + p0 ∆V .

Во втором цикле рассчитывается изотермический процесс при температуре, которую газ достиг в конце первого процесса. Применяются формулы: V

t +1

= V t + ∆V , p t +1 = BT / V t +1 , At +1 = At + p0 ∆V . Ре-

зультаты моделирования представлены на рис. 11.3. В столбцы A, B,
C, D таблицы выводятся объем, давление, температура и работа
газа. Суммарная работа, совершенная газом, равна 22 кДж.

Рис. 11.3. Изменение состояния газа. Графики p (V ) , V (T ) , p (T )
Программа 11.3
Private Sub CommandButton1_Click()
dv = 0.0001: V = 0.001: nu = 2: R = 8.31: T = 320

88

Майер Р. В. Решение физических задач в электронных таблицах Excel

B = nu * R: p0 = B * T / V: p = p0: i = 1
While V < 0.0025
Cells(i, 1) = V: Cells(i, 2) = p: Cells(i, 3) = T
V = V + dv: i = i + 1: p = p0: T = p * V / B
A = A + p * dv: Cells(i, 4) = A
Wend
While V < 0.007
Cells(i, 1) = V: Cells(i, 2) = p: Cells(i, 3) = T
V = V + dv: i = i + 1: p = B * T / V: A = A + p * dv
Cells(i, 4) = A
Wend
End Sub

11.4. Два моля газа имеют объем 1 л при температуре 300 К. В ходе цикла Карно газ: 1) изотермически расширяется до 3 л; 2) адиабатически расширяется до 7,5 л; 3) изотермически сжимается до 2 л; 4)
адиабатически сжимается до начального объема. Постройте график в
осях p и V , вычислите работу газа. Показатель политропы γ = 1,67 .
Для изотермического процесса pV = νRT = C . Рассчитав C ,
можно определить начальное давление p0 = C /V0 . При изотермическом расширении до 3 л: V t +1 = V t + ∆V , p t +1 = С / V t +1 . При адиабатном расширении от 3 до 7,5 л: D = pV γ = const , E = TV γ −1 = const ,

V t +1 = V t + ∆V , p t +1 = D /(V t +1 ) γ , T t +1 = E /(V t +1 ) γ −1 . При изотермическом сжатии от 7,5 л до 2 л: C = νRT , V t +1 = V t − ∆V , p t +1 = С / V t +1 .
При адиабатном сжатии от 2 до 1 л: D = pV γ = const , E = TV γ −1 ,

V t +1 = V t − ∆V ,

p t +1 = D /(V t +1 ) γ , T t +1 = E /(V t +1 ) γ −1 . Используется

программа 11.4, результаты представлены на рис. 11.4. При расширении газ совершает положительную работу A1→3 > 0, а при сжатии –
отрицательную работу A3→1 < 0; суммарная работа равна 2,5 кДж.

89

Майер Р. В. Решение физических задач в электронных таблицах Excel

Рис. 11.4. Цикл Карно в осях p и V . Зависимость работы от объема
Программа 11.4
Private Sub CommandButton1_Click()
dv = 0.0001: V = 0.001: nu = 2: R = 8.31: T = 300
gamma = 1.67: C = nu * R * T: p0 = C / V: p = p0: i = 1
While V < 0.003
Cells(i, 1) = V: Cells(i, 2) = p: Cells(i, 3) = T
V = V + dv / 10: i = i + 1: p = C / V
Wend: D = p * V ^ gamma: E = T * (V ^ (gamma - 1))
While V < 0.0075
Cells(i, 1) = V: Cells(i, 2) = p: Cells(i, 3) = T
V = V + dv: i = i + 1
p = D / (V ^ gamma): T = E / (V ^ (gamma - 1))
Wend: C = nu * R * T
While V > 0.002
Cells(i, 1) = V: Cells(i, 2) = p: Cells(i, 3) = T
V = V - dv: i = i + 1: p = C / V
Wend: D = p * V ^ gamma: E = T * (V ^ (gamma - 1))
While (p < p0)
Cells(i, 1) = V: Cells(i, 2) = p: Cells(i, 3) = T
V = V - dv / 10: i = i + 1
p = D / (V ^ gamma): T = E / (V ^ (gamma - 1))
Wend:
End Sub

11.5. Имеется однородный стержень длиной L с коэффициентом
теплопроводности α . Задано начальное распределение температуры
T (x) и мощности источников тепла q(x) . Необходимо рассчитать температуру различных точек стержня в произвольный момент времени τ ' .
Уравнение теплопроводности для стержня:

∂T
∂ 2T q
=α 2 +
.
∂τ
c
ρ
∂x

90

Майер Р. В. Решение физических задач в электронных таблицах Excel

Построим одномерную сетку с шагом h = ∆x = 1, шаг по времени
∆τ = 0,01. В конечных разностях получаем:

Ti

t +1

Tit−1 − 2Tit + Tit+1
qi
= Ti + α
∆
τ
+
∆τ .
cρ
∆x 2
t

Используется программа 11.5. Число оборотов цикла t ' , соответствующее времени τ ' = t '⋅∆τ , следует записать в ячейку E1. При запуске макрос создает таблицу из двух столбцов x и T (x) , а содержимое
ячейки E1 увеличивает на 500. На основе этой таблицы можно построить график T = T (x ) . При каждом следующем нажатии на кнопку,
выполняются очередные 500 итераций и рисуется новый график для
момента τ 'k +1 = (t 'k +500)∆τ . На рис. 11.5.1 приведен график T = T (x )
для случая: левый конец стержня поддерживается при постоянной
температуре, правый конец – теплоизолирован, имеется источник
тепла ( q > 0 ). Граничные условия: T1 = 4, T100 = T99 .
Программа 11.5
Private Sub CommandButton1_Click()
Dim t(100) As Single
Dim t1(100) As Single
dx = 1: dt = 0.01
For i = 1 To 100: Cells(i, 1) = i * dx
If i < 20 Then t(i) = 4 Else: t(i) = 0
Next i: Vremya = Cells(1, 5)
For k = 1 To Vremya: For i = 2 To 99
If (i > 75) And (i < 80) Then q = 0.2 Else q = 0
t1(i) = t(i) + 1.5 * (t(i - 1) - 2 * t(i) + t(i + 1)) * dt / dx /
dx + q * dt
Next i
For i = 1 To 100: t(i) = t1(i): Next i
t(1) = 4: t(100) = t(99): Next k
For i = 1 To 100: Cells(i, 2) = t(i): Next i
Cells(1, 5) = Cells(1, 5) + 500
End Sub

Рис. 11.5. Распределение температуры вдоль стержня

91

Майер Р. В. Решение физических задач в электронных таблицах Excel

11.6. Используя решение предыдущей задачи, рассчитайте T (x)
для однородного стержня в различные моменты времени, если левый
конец теплоизолирован, правый поддерживается при постоянной температуре, имеется источник тепла ( q > 0 ) и источник холода ( q < 0 ). Пример получающегося распределения температуры показан на рис. 11.5.
11.7. Имеется однородная пластина размером Lx х L y с коэффициентом теплопроводности α . Задано начальное распределение температуры T ( x, y ) и мощности источников тепла q( x, y ) . Необходимо рассчитать температуру различных точек пластины в произвольный момент τ ' .
Уравнение теплопроводности для пластины:

 ∂ 2T ∂ 2T  q
∂T
= α  2 + 2  +
.
∂τ
c
ρ
∂
x
∂
y


Задача решается аналогично. Дискретизируют двумерную область,
переходя к сетке N x M с шагом h = ∆x = ∆y , и записывают конечноразностное уравнение:

 Tit−1, j − 2Tit, j + Tit+1, j Tit, j −1 − 2Tit, j + Tit, j +1 
q
∆τ + i , j ∆τ .
+α
+


cρ
∆x 2
∆y 2


Используется программа 11.6. Время τ ' записывают в ячейку

Tit, +j 1

= Tit, j

A33 и нажимают на кнопку запуска макроса. На экране появляется
двумерная матрица значений температуры в узлах сетки в момент
τ ' (рис. 11.6). При повторном запуске программы выполняются следующие 25 итераций, вычисляются значения T (i, j ) в узлах двумерной
сетки для момента τ ' k +1 = (t ' k +25)∆τ . При необходимости можно отрегулировать ширину столбцов и вручную изменить цвет ячеек, температура которых находится в различных диапазонах. На рис. 11.6.1
показано распределение температуры T (i, j ) при наличии одного источника тепла; на рис. 11.6.2 – при наличии источника тепла ( q > 0 ) и
источника холода ( q < 0 ).
Программа 11.6
Private Sub CommandButton1_Click()
N = 30: M = 30: h = 1: dt = 0.01
Dim t(30, 30) As Single
Dim t1(30, 30) As Single
For i = 1 To N: For j = 1 To M

92

Майер Р. В. Решение физических задач в электронных таблицах Excel

If i < 20 Then t(i, j) = 10 Else t(i, j) = 0
Next j: Next i: Vremya = Cells(33, 1)
For k = 1 To Vremya
For i = 2 To N - 1: For j = 2 To M - 1
If (Abs(i - 12) < 5) And (Abs(j - 15) < 6)
Then q = 10 Else q = 0
AA = t(i - 1, j) - 4 * t(i, j) + t(i + 1, j) +
t(i, j - 1) + t(i, j + 1): t1(i, j) = t(i, j) +
0.5 * (AA) * dt / h / h + q * dt
Next j: Next i
For i = 2 To N: For j = 2 To M
t(i, j) = t1(i, j): Next j: Next i: Next k
For i = 2 To N: For j = 2 To M
Cells(i, j) = t(i, j): Next j: Next i
Cells(33, 1) = Cells(33, 1) + 25
End Sub

Рис. 11.6. К задаче о теплопроводности пластины

93

Майер Р. В. Решение физических задач в электронных таблицах Excel
11. Термодинамика

Содержание

13. Электрические цепи
постоянного тока

12. ЭЛЕКТРИЧЕСКОЕ И МАГНИТНОЕ ПОЛЕ
12.1. На оси Ox в точках с координатами x1 и x2 расположены два
точечных заряда q1 и q2 . Рассчитайте потенциал ϕ и проекцию напряженности E x электрического поля на ось Ox в различных ее точках и
постройте графики зависимости ϕ = ϕ (x ) и E x = E x (x) .
Потенциал и проекция напряженности электрического поля на
ось Ox вычисляются так:

ϕ=k

q1
q2
+k
,
| x − x1 |
| x − x2 |

E1 = k

| q1 |
,
2
( x − x1 )

E2 = k

| q2 |
.
2
( x − x2 )

Если заряды q1 и q2 положительны, то:

 E , если x > x1,
 E , если x > x2 ,
E1x =  1
E2 x =  2
Ex = E1x + E2 x .
−
E
,
−
E
,
x
<
x
;
x
<
x
;
если
если
 1
 2
1
2
Используется программа 12.1. Она содержит цикл по x , в котором с
шагом 0,01 перебираются все точки заданного интервала [0,1; 0,7],
вычисляются значения ϕ (x) и E x (x) . Результаты заносятся в таблицу, на их основе строятся графики. Рис. 12.1.1 и 12.1.3 соответствуют положительным q1 и q2 , а рис. 12.1.2 и 12.1.4 – отрицательному

q1 и положительному q2 .

Рис. 12.1. Получающиеся графики зависимостей ϕ (x) и E x (x)

94

Майер Р. В. Решение физических задач в электронных таблицах Excel

Программа 12.1
Private Sub CommandButton1_Click()
q1 = 5: q2 = 15: k = 0.009
x1 = 0.3: x2 = 0.5: x = 0.1: dx = 0.01
While x < 0.7
x = x + dx: i = i + 1
r1 = Abs(x - x1): r2 = Abs(x - x2)
fi = k * q1 / r1 + k * q2 / r2
If fi > 10 Then fi = 10
If fi < -10 Then fi = -10
E1 = k * q1 / r1 / r1: E2 = k * q2 / r2 / r2
If x < x1 Then E1 = -E1
If x < x2 Then E2 = -E2
E = E1 + E2: If E > 200 Then E = 200
If E < -200 Then E = -200
Cells(i, 1) = x: Cells(i, 2) = fi: Cells(i, 3) = E
Wend
End Sub

12.2. Два точечных заряда 10 нКл и –15 нКл расположены на расстоянии 1 м в точках с координатами (0 ; 0) и (1; 0) (в метрах). Рассчи-

r

тайте вектор напряженности E электрического поля и потенциал ϕ в
точке C, имеющей координаты (0,7; 0,5). Диэлектрическая проницаемость среды ε = 2,6 .

Рис. 12.2. К расчету напряженности и потенциала поля двух зарядов
По принципу суперпозиции

r r r
E = E1 + E2 . Введем обозначения:

L =| AB | , r1 =| AC | , r2 =| BC | . Чтобы рассчитать напряженность
электростатического поля двух зарядов, используются формулы:

E1 = k

q1
q2
,
E
=
k
, r1 = x 2 + y 2 , r2 = ( L − x ) 2 + y 2 ,
2
2
2
ε ⋅ r1
ε ⋅ r2

cos α = x / r1 , sin α = y / r1 , cos β = ( L − x) / r2 , sin β = y / r2 ,
E x = E1x + E2 x , E y = E1 y + E2 y , E = Ex2 + E y2 , γ = arctan( E y / E x ) .

95

Майер Р. В. Решение физических задач в электронных таблицах Excel

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

r

r

напряженности E в точке С и угол

γ , между вектором E и осью Ox .

Результаты должны выводиться в другие ячейки таблицы. Для расчета потенциала используется формула ϕ = k (q1 / r1 + q2 / r2 ) / ε .
12.3. По двум параллельным направляющим скользит металлический стержень, скорость которого изменяется по закону υ x (τ ) = 0,2 +

+ 0,3 ⋅τ 2 (м/c). Направляющие замкнуты через резистор сопротивлением
14 Ом, расстояние между ними 16 см. Система находится в магнитном

r

поле с индукцией B ( x ) = 0,6 + 0,2 x (Тл), причем вектор B перпендикулярен плоскости движения стержня. Как изменяется ЭДС индукции и сила
тока в промежутке времени от 0 до 2,4 с ?

Рис. 12.3. К расчету ЭДС и силы тока в движущемся стержне
При движении стержня в нем возникает ЭДС индукции ei = Bυl и
течет ток i = ei / R . Программа должна содержать цикл по времени, в
котором вычисляются скорость υ (ti ) , координата x(ti ) , индукция
магнитного поля Bx , ЭДС индукции e(ti ) и сила тока i (ti ) в последовательные моменты времени ti = i ⋅ ∆τ . Используются формулы:

υ t = 0,2 + 0,3 ⋅τ 2 , xt = xt −1 + υ t ∆τ , eit = υ t B ( x t ) L , i t = eit / R .
Задача решается с помощью программы 12.2; получающиеся графики
представлены на 12.4.
Программа 12.2
Private Sub CommandButton1_Click()
L = 0.16: R = 14: dt = 0.01
While t < 2.4
t = t + dt: i = i + 1: v = 0.2 + 0.3 * t * t
x = x + v * dt: B = 0.6 + 0.2 * x: eds = v * B * L: tok = eds / R

96

Майер Р. В. Решение физических задач в электронных таблицах Excel

If i Mod 10 = 0 Then
Cells(i / 10, 1) = t: Cells(i / 10, 2) = eds
Cells(i / 10, 3) = tok
End If
Wend
End Sub

Рис. 12.4. Зависимость ЭДС индукции и силы тока от времени
12.4. Три параллельных проводника с токами I1 , I 2 и I 3 расположены перпендикулярно плоскости xOy и пересекают ее в точках с координатами ( xi ; yi ), i = 1,2,3 . Постройте силовую линию магнитного поля,
которая проходит через заданную точку С0 ( x0 ; y0 ).
Сначала решим более простую задачу. Пусть два проводника с
токами I1 и I 2 пересекают плоскость xOy в точках A1 и A2 с координатами ( x1 , y1 ) и ( x2 , y 2 ). Если бы I1 и I 2 были направлены от нас, то
индукция магнитного поля в точке C (рис. 12.5.1) вычислялась так:

r r r
B = B1 + B2 , Bx = − B1 sin α1 − B2 sin α 2 , B y = B1 cos α1 + B2 cos α 2 ,

cos α1 = ( x − x1 ) / r1 , sin α1 = ( y − y1 ) / r1 , B1 = kI1 / r1 ,
cos α 2 = ( x − x2 ) / r2 , sin α 2 = ( y − y 2 ) / r2 , B2 = kI 2 / r2 ,
r1 = ( x − x1 ) 2 + ( y − y1 ) 2 , r2 = ( x − x2 ) 2 + ( y − y 2 ) 2 .

97

Майер Р. В. Решение физических задач в электронных таблицах Excel

r
Рис. 12.5. К расчету силовой линии магнитного поля B
Если ток направлен на нас, то I i и Bi берутся со знаком минус. Сило-

r

вые линии B для трех токов показаны на рис. 12.5.2. В используемой
программе 12.3 задаются координаты точек A1 , A2 , A3 плоскости

xOy , через которые проходят три параллельных проводника, и текущие по ним токи I1 , I 2 , I3 . В заданной точке C0 определяются проек-

r

ции вектора индукции B на оси координат, и в его направлении строится отрезок небольшой длины
a =| С0С1 | по формулам:

xi +1 = xi + a cos α , yi +1 = yi + a sin α , cos α = B x / B , sin α = B y / B , где
r
α – угол между вектором B и Ox . Определяются координаты точки
C1 , результат записывается в таблицу, и все повторяется снова.
Выполнив 10000 шагов, программа заканчивает свою работу. По полученным данным строится силовая линия (рис. 12.6).
Программа 12.3
Private Sub CommandButton1_Click()
x1 = 150: y1 = 100: x2 = 400: y2 = 100: x3 =
y3 = 300: I1 = 1: I2 = -2: I3 = 2: x = 40: y
While k < 10000
r1 = Sqr((x - x1) * (x - x1) + (y - y1) * (y
r2 = Sqr((x - x2) * (x - x2) + (y - y2) * (y
r3 = Sqr((x - x3) * (x - x3) + (y - y3) * (y
B1 = I1 / r1: B2 = I2 / r2: B3 = I3 / r3
cosa1 = (x - x1) / r1: sina1 = (y - y1) / r1
cosa2 = (x - x2) / r2: sina2 = (y - y2) / r2
cosa3 = (x - x3) / r3: sina3 = (y - y3) / r3
Bx = -B1 * sina1 - B2 * sina2 - B3 * sina3
By = B1 * cosa1 + B2 * cosa2 + B3 * cosa3
B = Sqr(Bx * Bx + By * By): x = x + 0.2 * Bx
y = y + 0.2 * By / B: k = k + 1

98

300
= 20
- y1))
- y2))
- y3))

/ B

Майер Р. В. Решение физических задач в электронных таблицах Excel

Cells(k, 1) = x: Cells(k, 2) = y: Cells(k, 3) = x1
Cells(k, 4) = y1: Cells(k, 5) = x2: Cells(k, 6) = y2
Wend
End Sub

Рис. 12.6. Результаты расчета силовой линии магнитного поля
12.5. Рассчитайте стационарное распределение потенциала электростатического поля между двумя бесконечно большими пластинами,
расположенными перпендикулярно оси Ox , на расстоянии b = 1 м друг
от друга. Плотность заряда зависит от координаты так:

5( x − 0,10), если | x − 0,25 |< 0,15;

ρ ( x) =  − 0,8,
если | x − 0,70 |< 0,05;

0,
в противном случае.

Потенциалы пластин равны ϕ1 = 120 В, ϕ N = – 80 В.
Требуется решить уравнение Пуассона для одномерной среды:

ρi
ρ ( x) ϕi −1 − 2ϕi + ϕ i +1
ρi ∆x 2 
1  k
k +1
k
=−
,
= − , ϕi =  ϕi −1 + ϕi +1 +
,
2
ε 
ε
ε
∂x 2
∆x 2

∂ 2ϕ

где k – номер итерации. Распределение плотности заряда задается
с помощью одномерного массива rho(i). В цикле перебираются узлы
одномерной сетки и вычисляются новые значения потенциала, которые записываются в массив f1(i). Затем происходит переобозначение
этого массива в f(i) и осуществляется новая итерация. После окончания заданного в ячейке E1 количества итераций результаты вычислений выводятся в столбец B. По ним строится график ϕ (x) , соот-

99

Майер Р. В. Решение физических задач в электронных таблицах Excel

ветствующий данному количеству итераций (рис. 12.7). При следующем нажатии на кнопку CommandButton1 выполняются следующие
2000 итераций. После того как программа выполнит большое количество итераций, получающееся распределение потенциала ϕ (x) перестанет изменяться и будет близко к истинному.
Программа 12.4
Private Sub CommandButton1_Click()
Dim f(100) As Single: Dim f1(100) As Single
Dim rho(100) As Single: dx = 1: eps = 1
For i = 1 To 100: Cells(i, 1) = i * dx
If Abs(i - 25) < 15 Then rho(i) = 0.05 * (i - 10) Else rho(i) = 0
If Abs(i - 70) < 5 Then rho(i) = -0.8
Next i: Iter = Cells(1, 5)
For k = 1 To Iter: For i = 2 To 99
f1(i) = (f(i - 1) + f(i + 1) + rho(i) * dx * dx / eps) / 2
Next i
For i = 1 To 100: f(i) = f1(i): Next i
f(1) = 120: f(100) = -80: Next k
For i = 1 To 100: Cells(i, 2) = f(i)
Cells(i, 3) = rho(i) * 200: Next i
Cells(1, 5) = Cells(1, 5) + 2000
End Sub

Рис. 12.7. Расчет распределения потенциала в одномерной области

100

Майер Р. В. Решение физических задач в электронных таблицах Excel

12.6. Рассчитайте стационарное распределение потенциала электростатического поля между двумя бесконечно большими пластинами,
расположенными перпендикулярно оси Ox , на расстоянии b = 1 м друг
от друга. Распределение плотности заряда задано так:

 0,02, если | x − 0,30 |< 0,08;

ρ (x) = − 0,3, если | x − 0,80 |< 0,04;
 0,
в противном случае.


ϕ1 = –140 В, а правая – ϕ N = 40 В.

Левая пластина имеет потенциал

Среда неоднородная, диэлектрическая проницаемость

4, если x < 0,50;
1, если x ≥ 0,50.

ε (x) = 

Распределение потенциала в одномерной неоднородной среде
удовлетворяет следующему дифференциальному уравнению:

∂
∂ϕ 
 ε ( x)  = − ρ ( x) ,
∂x 
∂x 

∂ε ∂ϕ
∂ 2ϕ
⋅
+ ε 2 = −ρ .
∂x ∂x
∂x

В конечных разностях получаем:

(ε i +1 − ε i −1 ) (ϕi +1 − ϕi −1 )
ϕ − 2ϕi + ϕi +1
⋅
+ ε i i −1
= − ρi ,
2∆x
2∆x
∆x 2

ϕik +1

1  (ε i +1 − ε i −1 ) (ϕik+1 − ϕik−1 )
k
k
2
ε
(
ϕ
ϕ
)
ρ
x
=
⋅
+
+
+
∆

.
i i −1
i +1
i
2ε i 
2
2


Эта формула позволяет рассчитать распределение потенциала на
(k + 1) –й итерации. При нажатии на кнопку CommandButton1 программа осуществляет 2000 итераций, при повторном нажатии – следующие 2000 и т. д. Получающиеся результаты представлены на
рис. 12.8.
Программа 12.5
Private Sub CommandButton1_Click()
Dim f(100) As Single: Dim f1(100) As Single
Dim rho(100) As Single: Dim e(100) As Single: dx = 1
For i = 1 To 100: Cells(i, 1) = i * dx
If Abs(i - 30) < 8 Then rho(i) = 0.2
If Abs(i - 80) < 4 Then rho(i) = -0.3
If i < 50 Then e(i) = 4 Else e(i) = 1
Next i: Iter = Cells(1, 5)
For k = 1 To Iter: For i = 2 To 99
f1(i) = ((e(i + 1) - e(i - 1)) * (f(i + 1) –
f(i - 1)) / 4 + e(i) * (f(i + 1) + f(i - 1))

101

Майер Р. В. Решение физических задач в электронных таблицах Excel

+ rho(i) * dx * dx) / 2 / e(i): Next i
For i = 1 To 100: f(i) = f1(i): Next i
f(1) = -140: f(100) = 40: Next k
For i = 1 To 100: Cells(i, 2) = f(i):
Cells(i, 3) = rho(i) * 200
Next i: Cells(1, 5) = Cells(1, 5) + 2000
End Sub

Рис. 12.8. Расчет распределения потенциала в одномерной области

102

Майер Р. В. Решение физических задач в электронных таблицах Excel
12. Электрическое
и магнитное поле

Содержание

14. Переменный ток.
Нестационарные процессы
в электрических цепях

13. ЭЛЕКТРИЧЕСКИЕ ЦЕПИ ПОСТОЯННОГО ТОКА
13.1. Определите сопротивление бесконечной цепочки резисторов,
изображенной на рис. 13.1, если r1 = 2 Ом и r2 = 40 Ом.

Рис. 13.1. Бесконечная цепочка резисторов
Рассмотрим сначала одну секцию A1B1 из последовательно соединенных резисторов r1 и r2 , затем две секции A2 A1B1B2 , три секции

A3 A2 A1B1B2 B3 , … N секций и т. д. Их сопротивления рассчитываются
по формулам:

R1 = r1 + r2 , R2 = r1 +

r2 R1
rR
rR
, R3 = r1 + 2 2 , …, RN = r1 + 2 N −1 .
r2 + R1
r2 + R2
r2 + RN −1

Для расчета сопротивления цепочки резисторов используется программа 13.1, получающийся график представлен на рис. 13.1. При
N → ∞ сопротивление RN приближается к ≈ 10,1 Ом.
Программа 13.1
Private Sub CommandButton1_Click()
r1 = 2: r2 = 40: RR = r1 + r2: i = 1
While i < 10
Cells(i, 1) = i: Cells(i, 2) = RR
i = i + 1: RR = r1 + r2 * RR / (r2 + RR)
Wend
End Sub

103

Майер Р. В. Решение физических задач в электронных таблицах Excel

Рис. 13.2. Зависимость сопротивления RN цепи от числа секций N
13.2. К источнику постоянного напряжения с ЭДС E = 10 В и внутренним сопротивлением r = 5 Ом подключен переменный резистор. При
каком сопротивлении внешней нагрузки R выделяющаяся на нем мощность P максимальна? Постройте график P = P (R ) .
Сила тока и мощность рассчитываются по формулам:

I = E /( R + r ) ,

P = I 2R .

Решим эту задачу без использования макросов, для этого создадим
таблицу. Столбец A содержит значения сопротивления нагрузки R ,
столбцы B и C – результаты расчетов силы тока и мощности (в
ячейки B2 и С2 записывают “=10/(A2+5)”, “=B2*B2*A2”). Результаты
вычислений представлены на рис. 13.3. Видно, что на нагрузке выделяется максимальная мощность, когда ее сопротивление равно внутреннему сопротивлению источника 5 Ом.

Рис. 13.3. Расчет зависимостей I (R ) и P (R )
104

Майер Р. В. Решение физических задач в электронных таблицах Excel

13.3. К источнику постоянной ЭДС с внутренним сопротивлением
r0 подключен реостат сопротивлением R , использующийся в качестве
делителя напряжения (рис. 13.4). К реостату подключена нагрузка RH .
Изучите зависимость тока I через источник, тока I H и напряжения U H
на нагрузке от координаты x подвижного контакта реостата.
Пусть длина обмотки ab реостата равна L . Сопротивление
участка ac обозначим через r1 , а cb – через r2 . Если k = x / L , то

r2 = kR , r1 = R − r2 . Используются формулы
Roб = r0 + r1 + r2 RH /(r2 + RH ) , I = E / Rоб ,
U H = E − I (r0 + r1 ) ,

I H = U H / RH .

Используется программа 13.2, результаты – на рис. 13.4.
Программа 13.2
Private Sub CommandButton1_Click()
E = 100: R = 100: RN = 20: L = 100
For i = 1 To 100
k = 0.01 * i: r2 = R * k: r1 = R - r2: x = k * L
Rob = r0 + r1 + r2 * RN / (r2 + RN)
tok = E / Rob: UN = E - tok * (r0 + r1)
Inag = UN / RN: Cells(i, 1) = k
Cells(i, 2) = x: Cells(i, 3) = tok
Cells(i, 4) = UN: Cells(i, 5) = Inag
Next
End Sub

Рис. 13.4. Результаты моделирования работы делителя напряжения
105

Майер Р. В. Решение физических задач в электронных таблицах Excel

13.4. К источнику постоянной ЭДС E = 2 В подключены последовательно соединенные резистор R = 50 Ом и нелинейный элемент, вольтамперная характеристика которого может быть аппроксимирована функ2
цией I (U НЭ ) = 0,01 ⋅ U НЭ
. Необходимо рассчитать ток в цепи и напряже-

ние на элементах.
Так как E = U НЭ + IR , то напряжение на нелинейном элементе
равно U НЭ = E − IR . Для нахождения тока I в цепи следует решить
нелинейное алгебраическое уравнение: I = 0,01( E − IR ) . Искомое зна2

чение тока лежит в интервале от 0 до E / R = 0,04 (A). Построим график функции F ( I ) = 0,01( E − IR ) − I и найдем значение тока I , при
2

котором функция F (I ) обращается в ноль. Используется программа
13.3, результаты ее работы – на рис. 13.5.
Программа 13.3
Private Sub
di = 0.001:
While tok <
tok = tok +
U = E - tok
Cells(k, 1)
Cells(k, 3)
Wend
End Sub

CommandButton1_Click()
E = 2: R = 50: tok = 0
E / R
di: k = k + 1
* R: F = 0.01 * U * U - tok
= tok: Cells(k, 2) = F
= U

Рис. 13.5. График функции F ( I ) = 0,01( E − IR ) − I
Из получившихся графика и таблицы следует, что ток в цепи
будет равен 0,015 – 0,016 А, а напряжение на нелинейном элементе
2

106

Майер Р. В. Решение физических задач в электронных таблицах Excel

1,20 – 1,25 В. Повторите вычисления с меньшим шагом ∆I и найдите
более точное значение I . Напряжение на элементах: U R = IR ,

U НЭ = E − U R .
13.5. Используя решение предыдущей задачи, определите ток в
анализируемой цепи в случае, когда ЭДС источника E изменяется от 0
до 5 В. По полученным данным постройте график зависимости I = I (E )
в электронных таблицах Excel.
13.6. К источнику постоянной ЭДС E подключены последовательно
соединенные резистор R = 50 Ом и нелинейный элемент, вольтамперная характеристика которого может быть аппроксимирована многочленом третьей степени I (U ) = 0,03U 2 − 0,01U 3 (А). Определите ток в
цепи и напряжение на элементах при E = 1,8 В; при E = 2,4 В.

107

Майер Р. В. Решение физических задач в электронных таблицах Excel
13. Электрические цепи
постоянного тока

Содержание

15. Явления
геометрической оптики

14. ПЕРЕМЕННЫЙ ТОК.
НЕСТАЦИОНАРНЫЕ ПРОЦЕССЫ В ЭЛЕКТРИЧЕСКИХ ЦЕПЯХ
14.1. К источнику переменного напряжения регулируемой частоты
подключен последовательный колебательный контур, состоящий из резистора R , конденсатора C и катушки индуктивности L . Рассчитайте
емкостное X C , индуктивное X L , полное сопротивления цепи Z и силу
тока I на разных частотах f . Постройте резонансную кривую I = I ( f ) .
Для моделирования резонанса напряжений используют формулы:

ω = 2π ⋅ f , X L = ω L , X C = 1 /(ω ⋅ C ) , Z = R 2 + ( X L − X C ) 2 ,
I = U / Z , U L = X L I , U C = X C I , f 0 = 1 /(2π LC ) .
Здесь f 0 – собственная частота колебательного контура. Для того
чтобы можно было изменять параметры R , C , L , их значения следует поместить в ячейки L1, L2, L3 и использовать абсолютную адресацию. Например, в ячейку С2 следует ввести “=B2*$L$2”, а в ячейку
D2 нужно записать “=1/(B2*$L$3)”. Результаты моделирования представлены на рис. 14.1.

Рис. 14.1. Резонансная кривая в последовательном контуре

108

Майер Р. В. Решение физических задач в электронных таблицах Excel

14.2. Последовательный колебательный контур подключен к источнику переменного напряжения регулируемой частоты. Пользуясь решением предыдущей задачи, рассчитайте частотные зависимости напряжения на конденсаторе U C ( f ) , катушке индуктивности U L ( f ) и
сдвига фаз ϕ ( f ) между колебаниями тока и напряжения на полюсах источника.

Рис. 14.2. Моделирование резонанса напряжений
Таблицу на рис. 14.1 следует дополнить тремя столбцами, в которых рассчитываются U C , U L и ϕ на разных частотах по следующим формулам: U C = IX C , U L =IX L , ϕ = arctg[( X L − X C ) / R ] . Вблизи
собственной частоты контура f 0 ≈ 145 Гц имеет место резонанс
напряжений: X L ( f 0 ) = X C ( f 0 ) , U L ( f 0 ) = U C ( f 0 ) , а ϕ = 0 . Получаются
графики, изображенные на рис 14.2. Кривую ϕ ( f ) постройте самостоятельно.
14.3. К источнику переменного напряжения регулируемой частоты
подключен параллельный колебательный контур, состоящий из резистора R , конденсатора C и катушки индуктивности L . Рассчитайте проводимость конденсатора YC , катушки индуктивности YL , полную проводимость цепи Y и силу тока I через источник на разных частотах f .
Промоделируйте резонанс напряжений, постройте резонансную кривую
I = I ( f ) и график ϕ = ϕ ( f ) .
Для моделирования резонанса токов используются формулы:

109

Майер Р. В. Решение физических задач в электронных таблицах Excel

ω = 2π · f , YR = 1 / R , YL = 1 / ωL , YC = ω ⋅ C , Y = YR2 + (YC − YL ) 2 ,
I R = YRU , I L = YLU , I C = YCU , I = YU , ϕ = arctg[(YC − YL ) / YR ] .
Задача решается путем создания таблицы, в столбцах которой рассчитываются все искомые величины при различных частотах f . Из
рис. 14.3 и 14.4 видно, что на частоте примерно 150 Гц наступает
резонанс токов: проводимости YL , YC и токи I L , I C равны, общий
ток достигает минимума, сдвиг фаз обращается в 0.

Рис. 14.3. Резонанс токов в параллельном колебательном контуре

Рис. 14.4. Зависимость сдвига фаз от частоты

110

Майер Р. В. Решение физических задач в электронных таблицах Excel

14.4. Последовательно соединенные резистор 50 ом, конденсатор
150 мкФ и катушка индуктивности 0,03 Гн подключены к источнику прямоугольных импульсов амплитудой 90 В и частотой ω . Получите кривую
тока i (τ ) и напряжения на конденсаторе uC (τ ) и катушке индуктивности

u L (τ ) при различных частотах ω .

Из второго закона Кирхгофа:

e(τ ) = u R + uC + u L = iR +

 90, при sin ωτ > 0,
1
di
q + L , где e(τ ) = 
C
dτ
− 90, при sin ωτ ≤ 0.

В конечных разностях получаем:

i ' = (et − i t R − q t / C ) / L , i t +1 = i t +i ' ∆τ , q t +1 = q t + i t +1 ∆τ .
Используется программа 14.1; кривая тока i (τ ) представлена на
рис. 14.5. Учитывая, что uC = q / C и u L = Li ' , постройте графики
uC (τ ) и u L (τ ) .
Программа 14.1
Private Sub CommandButton1_Click()
Um = 90: w1 = 100: R = 50
L = 0.03: C = 0.00015: dt = 0.0001
While t < 0.25
t = t + dt: i = i + 1
If Sin(w1 * t) > 0 Then e = Um Else e
ii = (e - R * tok - q / C) / L
tok = tok + ii * dt: q = q + tok * dt
If i Mod 10 = 0 Then Cells(i / 10, 1)
If i Mod 10 = 0 Then Cells(i / 10, 2)
If i Mod 10 = 0 Then Cells(i / 10, 3)
Wend
End Sub

= -Um

= t
= e
= 20 * tok

Рис. 14.5. Получающиеся графики зависимостей e(τ ) и i (τ )

111

Майер Р. В. Решение физических задач в электронных таблицах Excel

14.5. Последовательно соединенные резистор R , конденсатор C и
катушка индуктивности L подключены к генератору импульсов. Промоделируйте работу цепи, рассчитайте кривые тока и напряжения на конденсаторе и катушке индуктивности, если генератор вырабатывает:
1) прямоугольные импульсы амплитудой U 0 длительностью T1 , разделенные промежутком T2 ; 2) пилообразные импульсы периодом T и амплитудой U 0 ; 3) импульсы, получающиеся при однополупериодном выпрямлении, имеющие период T и амплитуду U 0 .
Используется программа 14.1, в которую необходимо внести незначительные изменения. Чтобы задать пилообразные импульсы,
применяются операторы: U = U + 0.2: If U > Um Then U = Um. Результаты моделирования представлены на рис. 14.6 и 14.7.

Рис. 14.6. Прямоугольные импульсы e(τ ) . Графики u L (τ ) и uC (τ )

Рис. 14.7. Пилообразные импульсы e(τ ) . Графики q (τ ) и i (τ )

112

Майер Р. В. Решение физических задач в электронных таблицах Excel

14.6. Колебательный контур состоит из последовательно соединенных резистора R , конденсатора C и катушки индуктивности L . Конденсатор заряжают до напряжения U 0 и замыкают цепь. Определите
заряд конденсатора q0 и силу тока i в контуре в произвольный момент
времени t , постройте графики q (τ ) , i (τ ) .
В контуре происходит переходный процесс, характеризующийся
уравнением: iR + q / С + Ldi / dt = 0 . Для решения задачи необходимо в
программе 14.1 правильно задать напряжение питания e(τ ) = 0 и заряд конденсатора в начальный момент времени q(τ ) = q0 . Если R мало, то после замыкания цепи конденсатор периодически заряжается и
перезаряжается через катушку индуктивности, возникают затухающие колебания. При больших R сила тока через индуктивность и заряд конденсатора изменяются по апериодическому закону.
14.7. Колебательный контур состоит из последовательно соединенных резистора R , конденсатора C и катушки индуктивности L . Контур подключают к источнику ЭДС (рис. 14.8.1). Как изменяется заряд
конденсатора и ток через катушку индуктивности, если контур подключить к источнику постоянной ЭДС? К источнику ЭДС, вырабатывающему
прямоугольные импульсы напряжения?

Рис. 14.8. Переходный процесс; источник создает постоянную ЭДС
Используется программа 14.1. Если индуктивность L мала, а
сопротивление R и емкость C велики, то при замыкании ключа
(рис. 14.8.1) происходит заряд конденсатора: кривые q (τ ) и i (τ ) показаны на рис. 14.8.2. Если L и C велики, R мало, а источник вырабатывает прямоугольные импульсы напряжения, то в колебательном

113

Майер Р. В. Решение физических задач в электронных таблицах Excel

контуре возникает нестационарные процесс; кривые q (τ ) и i (τ ) показаны на рис. 14.9.

Рис. 14.9. Отклик цепи на прямоугольный импульс
14.8. Выпрямитель состоит из трансформатора, к вторичной обмотке которого подключен полупроводниковый диод и резистор нагрузки
R (рис. 14.10.1). Параллельно резистору включен конденсатор C , сглаживающий пульсации. Рассчитайте работу выпрямителя, постройте кривые тока iD (τ ) через диод и напряжения U N (τ ) на нагрузке.

Рис. 14.10. Выпрямитель с емкостным фильтром (1);
индуктивный фильтр (2)
По первичной обмотке трансформатора течет ток i1 (τ ) =

= I m sin(ωτ ) . ЭДС индукции во вторичной обмотке e2 (τ ) = Mdi1 / dτ .
Будем считать, что потенциал ϕO = 0 , тогда ϕ a = e2 (τ ) . Диод имеет
одностороннюю проводимость, его сопротивление равно:

500 /(ϕ a − ϕb ), если ϕ a > ϕb ,
RD = 
10000,
если ϕ a ≤ ϕb .

Через диод течет ток iD = (ϕ a − ϕb ) / RD . Токи через нагрузку и конденсатор равны i N = ϕ b / R и iC = iD − iN . За время ∆t конденсатор увели-

114

Майер Р. В. Решение физических задач в электронных таблицах Excel

чивает свой заряд на величину ∆qC = iC ∆τ . Потенциал точки b равен
напряжению на резисторе нагрузки: ϕ b = q / C = i N R = u N .
Для компьютерного моделирования работы выпрямителя с емкостным фильтром используется программа 14.2. Она содержит цикл
по времени, в котором вычисляются значения токов iD , iC , i N и напряжений u2 = e2 , u N в последовательные моменты времени, что
позволяет построить соответствующие графики (рис. 14.11).

Рис. 14.11. Результат моделирования выпрямителя с фильтром
Программа 14.2
Private Sub CommandButton1_Click()
N1 = 1000: dt = 0.0001: S = 0.001: R = 10: N2 = 1000:
M = 0.6: r2 = 5: Rn = 350: C = 0.002
While t < 3
t = t + dt: i = i + 1
i11 = i1: i1 = 10 * (1 – 0 * Exp(-t)) * Sin(t * 10)
e2 = M * (i1 - i11) / dt
If e2 > fia Then Rd = 500 / (e2 - fia) Else Rd = 10000
ID = (e2 - fia) / Rd: id3 = id2: id2 = id1: id1 = ID
inag = fia / Rn: ic = ID - inag
q = q + ic * dt: fia = q / C: un = inag * Rn
If i Mod 50 = 0 Then
Cells(i / 50, 1) = t: Cells(i / 50, 2) = 60 * ID
Cells(i / 50, 3) = 3 * un: Cells(i / 50, 4) = 3 * e2
End If
Wend
End Sub

115

Майер Р. В. Решение физических задач в электронных таблицах Excel

14.9. Напряжение с выпрямителя подается на резистор нагрузки R ,
последовательно с которым включена катушка индуктивности L
(рис. 14.10.2). Рассчитайте кривые тока i (τ ) и напряжения u N (τ ) на нагрузке в случаях: 1) однополупериодного выпрямления; 2) двухполупериодного выпрямления.
Катушка индуктивности играет роль индуктивного фильтра;
она сглаживает пульсации выпрямленного тока. Из второго закона
Кирхгофа следует

u (τ ) = u N + u L = iR + L

10 sin ωτ , при sin ωτ > 0,
di
, где u (τ ) = 
0,
при sin ωτ ≤ 0.
dτ


где u (t ) – пульсирующее напряжение, получающееся в результате
однополупериодного выпрямления. В случае двухполупериодного выпрямления u (τ ) = U m | sin ωτ | . Перейдем к конечным разностям:

i t +1 − i t
u =i R+L
,
∆τ
t

t

i

t +1

ut − it R
=i +
∆τ .
L
t

Используется программа 14.3, результаты – на рис.14.12.

Рис. 14.12. Кривые i (τ ) и u N (τ ) : двухполупериодное выпрямление
Программа 14.3
Private Sub CommandButton1_Click()
dt = 0.0001: Rn = 10: L = 0.4
While t < 0.3
t = t + dt: i = i + 1: u = 20 * Abs(Sin(100 * t))
tok = tok + (u - tok * Rn) * dt / L: un = tok * Rn
If i Mod 10 = 0 Then
Cells(i / 10, 1) = t: Cells(i / 10, 2) = 10 * u
Cells(i / 10, 3) = 200 * tok: Cells(i / 10, 4) = 10 * un
End If
Wend
End Sub

116

Майер Р. В. Решение физических задач в электронных таблицах Excel

14.10. Промоделируйте работу однополупериодного выпрямителя.
Изучите зависимость амплитуды пульсаций выпрямленного напряжения
от частоты f переменного напряжения, емкости C сглаживающего конденсатора и сопротивления R нагрузки.
Для решения задачи используется программа 14.2. В ней необходимо правильно задать u (t ) . Изменяя f , C и R , можно установить,
что при увеличении частоты f , емкости C и сопротивления R амплитуда пульсаций уменьшается.
14.11. Создайте компьютерную модель двухполупериодного выпрямителя с индуктивным фильтром, состоящим из катушки индуктивности, соединенной последовательно с нагрузкой. Изучите зависимость
амплитуды пульсаций от частоты напряжения, индуктивности катушки и
сопротивления нагрузки.
Для моделирования выпрямителя с индуктивным фильтром используется программа 14.3. Проводя серию вычислительных экспериментов при различных частоте f , индуктивности L , сопротивлении

R , можно убедиться в том, что при увеличении f и L и уменьшении
R амплитуда пульсаций уменьшается.

117

Майер Р. В. Решение физических задач в электронных таблицах Excel
14. Переменный ток.
Нестационарные процессы
в электрических цепях

Содержание

16. Волновые явления

15. ЯВЛЕНИЯ ГЕОМЕТРИЧЕСКОЙ ОПТИКИ
15.1. Луч света падает на границу раздела стекло-воздух. Постройте график зависимости угла преломления от угла падения, если показатель преломления стекла 1,5.
15.2. Решите предыдущую задачу для случая, когда луч света проходит из воздуха в слой воды, а затем в стекло. Показатели преломления воды и стекла соответственно равны 1,33 и 1,5.
15.3. Луч света падает под углом α на плоскопараллельную стеклянную пластину толщиной d = 4 см. Постройте график зависимости
смещения луча от угла падения. Показатель преломления стекла 1,5.
15.4. Согласно принципу Ферма свет распространяется по пути,
требующему экстремальное время. Рассчитайте время прохождения
светом пути АCB при преломлении на плоской границе раздела двух
сред MN (рис. 15.1.1) при различных координатах точки падения C. Убедитесь в том, что выполняется закон преломления sin α / sin β = n2 / n1 .

Рис. 15.1. Принцип Ферма при преломлении и отражении
В первой среде свет распространяется со скоростью υ1 = c / n1 , а
во второй – со скоростью υ 2 = c / n2 . Пусть | MN |= L , | AM |= a ,

| BN |= b (рис. 15.1.1). Время прохождения светом пути АСB равно:

t=

AC

υ1

+

CB

υ2

=

a2 + x2

υ1
118

+

( L − x) 2 + b 2

υ2

.

Майер Р. В. Решение физических задач в электронных таблицах Excel

Истинному ходу луча соответствует закон преломления, из которого следует: n1 sin α = n2 sin β . Ниже представлена программа, вычисвремя t распространения света, а также функции
f1 ( x) = n1 sin α и f 2 ( x) = n2 sin β . На основе полученных результатов

ляющая

строятся графики (рис. 15.2). Видно, что истинному ходу луча, при
котором f1 ( x) = f 2 ( x) (то есть выполняется закон преломления), соответствует минимальное время.
Программа 15.1
Private Sub CommandButton1_Click()
n1 = 1: n2 = 1.8: c = 3000000000#
L = 100: a = 70: b = 70: v1 = c / n1: v2 = c / n2
For x = 1 To 100
s1 = Sqr(a * a + x * x): s2 = Sqr(b * b + (L - x) * (L - x))
sinalfa = x / s1: sinbeta = (L - x) / s2
Cells(x, 1) = x
Cells(x, 2) = (s1 / v1 + s2 / v2) * 30000000
Cells(x, 3) = n1 * sinalfa: Cells(x, 4) = n2 * sinbeta
Next
End Sub

Рис. 15.2. Проверка принципа Ферма для преломления света

119

Майер Р. В. Решение физических задач в электронных таблицах Excel

15.5. Докажите, что при отражении света от плоской границы выполняется принцип Ферма, то есть свет затрачивает минимальное время. Рассчитайте время прохождения светом пути АCB при отражении от
зеркала MN (рис. 15.1.2) при различных координатах точки отражения С.
Убедитесь, что при этом выполняется закон отражения α = β .
15.6. На стеклянный полуцилиндр радиусом R = 10 см падает луч,
идущий параллельно оси Ox на расстоянии a от нее (рис. 15.3.1). Определите координату точки x пересечения луча и оси Ox после выхода
из стекла при различных a . На листе бумаги постройте ход лучей при
a = 0, 1, 3, 5, 7 см. Показатель преломления стекла 1,5.

Рис. 15.3. Преломление и отражение на сферической поверхности
Когда a не очень велико, из закона преломления для точки A следует: sin β = n sin α . Из геометрических соображений: sin α = a / R ,

γ = β − α . Координата точки B равна x =| OC | + | CB |= R cos α + a ⋅ ctgγ .
В столбец A электронной таблицы помещают значения параметра a , изменяющиеся от 0 до 10 см с шагом 0,1–0,5. В ячейки
B2–G2 записывают формулы:
B2: “=A2/10”; C2: “=1,5*B2”; D2: “=ASIN(B2)”; E2: “=ASIN(C2)”;
F2: “=E2-D2”; G2: “=10*COS(D2)+A2/TAN(F2)”.
Эти формулы копируют вниз, получая таблицу, похожую на таблицу
15.1. При достаточно больших значениях смещений луча a от оси Ox
в точке А происходит полное внутреннее отражение и луч света не
выходит из стекла. При этом n sin α оказывается больше 1.

120

Майер Р. В. Решение физических задач в электронных таблицах Excel

Таблица 15.1

15.7. На цилиндрическое (или сферическое) вогнутое зеркало радиусом R = 10 см падает луч, идущий параллельно главной оптической
оси Ox на расстоянии a от нее (рис. 15.3.2). Определите координату x
точки B пересечения луча и оси Ox после отражения. На листе бумаги
постройте ход лучей при a = 0, 1, 3, 5, 7 и 9 см.
Из геометрических соображений: sin α = a / R , α = arcsin(a / R ) . Если 2α ≥ π / 2 , то xB =| OC | + | CB |= R cos α + atg (2α − π / 2) . Случай при

2α < π / 2 проанализируйте самостоятельно. В столбец A помещают
значения параметра a , изменяющиеся от 0 до 10 с шагом 0,1–0,5. В
ячейки B2 и C2 записывают формулы: B2: “=ASIN(A2/10)”; C2: “=10*
COS(B2)+A2*TAN(2*B2-3,14/2)”. Эти формулы копируют вниз, получая
таблицу. Из результатов вычислений видно, что при увеличении a
от 0 до R ≈ 8,6 см точка B сначала медленно, а потом быстро смещается вправо, ее координата увеличивается от 0 до R . При R > 8,8 см
световой луч дважды отражается от зеркала и только потом достигает оси Ox . На листе бумаги следует нарисовать цилиндрическое
зеркало, падающие и отраженные лучи, соответствующие различным
значениям смещения луча a .
15.8. Кусок стекла ограничен полусферической поверхностью радиусом R = 10 см (рис. 15.4.1). На эту поверхность падает луч света,
идущий параллельно оси симметрии Ox на расстоянии a от нее. Най-

121

Майер Р. В. Решение физических задач в электронных таблицах Excel

дите координату точки B (x ) пересечения луча и оси Ox после преломления. На листе бумаги постройте ход лучей при различных a .
Для решения задачи следует создать таблицу, содержащую
столбцы a , sin α , α , sin β , β , x . При этом используются формулы:

sin α = a / R , sin β = sin α / n , | OC |= R cos α , ∠CAB = (π / 2 − α ) + β = γ ,
| CB |= atgγ , x =| CB | − | OC |= a ⋅ tg (π / 2 − α + β ) − R cos α .

Рис. 15.4. К расчету хода лучей в задачах 15.8 и 15.9
15.9. На стеклянный шар радиусом R падает пучок света
(рис. 15.4.2). Рассчитайте ход луча, смещенного от оси симметрии шара
на расстояние a . Найдите точку пересечения F луча с оптической осью
шаровой линзы. Постройте ход лучей, для которых смещение a принимает значения 0,2 R , 0,4 R , 0,6 R , 0,8 R .
Точка A вхождения светового луча в шар имеет координаты
x A = − R cos α , y A = a , где α – угол падения, его синус равен

sin α = a / R . По закону преломления sin α / sin β = n , sin β = sin α / n , где
β – угол преломления в точке A, равный углу падения в точке B. Координаты точки B выхода светового луча из шара находятся так:

∠CAO = 90 0 − α , ∠CAB = ∠CAO + β , ∠ABD = 360 0 − ∠CAB − 90 0 − 90 0 ,
∠OBD = ∠ABD − β , x B =| OD |= R sin(∠OBD) , y B =| BD |= R cos(∠OBD) .
При этом учитывается, что сумма углов прямоугольника ABDC равна
360 0 . Из закона преломления, записанного для точки B, следует, что
угол преломления луча света при выходе из шара равен α . Так как
угол является внешним по отношению к треугольнику OBF, то
γ + ∠BOD = α . Отсюда определяем угол γ и вычисляем координату
точки F по формуле: x F = R sin(∠OBD ) + y B ctgγ . Задавая различные

122

Майер Р. В. Решение физических задач в электронных таблицах Excel

смещения a , можно рассчитать и построить ход лучей внутри шара
и после выхода из него.
15.10. На прямую треугольную стеклянную призму с преломляющим углом ϕ = 0,54 рад под углом α1 подает световой луч (рис. 15.5).
Рассчитайте его угол отклонения ψ при различных значениях α1 , постройте график зависимости ψ = ψ (α1 ) . При каком угле падения α1 угол
отклонения ψ минимален?

Рис. 15.5. Отклонение луча света призмой на угол ψ
Для решения задачи используются формулы:

sin α1
sin α 2 1
= n,
= , β1 + α 2 = ϕ , ψ = (α1 − β1 ) + ( β 2 − α 2 ) .
sin β1
sin β 2 n
При этом учитывается, что: 1) угол ϕ для треугольника ABM является внешним и поэтому равен сумме двух других углов несмежных с
ним; 2) угол отклонения луча ψ призмой складывается из углов NAB и
NBA, которые равны α1 − β1 и β 2 − α 2 соответственно. Угол ψ вы-

β1 = arcsin(sin α1 / n) ,
α 2 = ϕ − β1 ,
β 2 = arcsin(n sin α 2 ) , ψ = (α1 − β1 ) + ( β 2 − α 2 ) .
Угол падения α1 изменяется от 0 до π / 2 ≈ 1,57 рад с шагом
∆α1 =0,05 рад; его записывают в первый столбец таблицы. В остальных столбцах находятся результаты вычислений β1 , α 2 , β 2 , ψ в радианах при различных α1 . Значения показателя преломления n и преломляющего угла ϕ задаются в ячейках H1 и H2 (рис. 15.6).
числяется так:

123

Майер Р. В. Решение физических задач в электронных таблицах Excel

Рис. 15.6. Получающиеся графики зависимостей β1 , α 2 , β 2 , ψ от α1
Из результатов вычислений видно, что при α1 =0,95 рад луч света падает на вторую преломляющую поверхность нормально
( α 2 = β 2 = 0 ). При любых α1 из интервала [0; π / 2 ] угол отклонения ψ
положительный и достигает минимума при α1 =0,40-0,45 рад. Чтобы
определить положение минимума точнее, необходимо повторить вычисления при меньшем шаге ∆α1 .
15.11. Решите задачу 15.10 для преломляющего угла ϕ = 0,83 рад.

При каком угле падения происходит полное внутреннее отражение?
15.12. Две линзы L1 и L2 с оптической силой D1 и D2 установлены
на расстоянии l друг от друга так, что их главные оптические оси совпадают. На расстоянии a1 от линзы L1 расположен предмет высотой h .
Определите высоту H и расположение его изображения.
Линза L1 создает изображение h' , которое является предметом
для линзы L2 (рис. 15.7). Из формулы тонкой линзы следует:

1 1 1
1 1
1
+ = , D2 =
+ =
.
a1 b1 F1
a2 b2 F2
aF
a F
Найдем b2 и H : b1 = 1 1 , a2 = l − b1 , b2 = 2 2 ,
a1 − F1
a2 − F2
D1 =

124

Майер Р. В. Решение физических задач в электронных таблицах Excel

h' b1
= ,
h a1

H b2
bb
= , H = 1 2 h.
h ' a2
a1a2

Для решения задачи можно использовать программу 15.2.

Рис. 15.7. Ход лучей в оптической системе из двух линз

Рис. 15.8. Таблица в документе Excel к задаче 15.12
Программа 15.2
Private Sub CommandButton1_Click()
a1 = Cells(1, 2): F1 = Cells(2, 2)
F2 = Cells(3, 2): L = Cells(4, 2)
h = Cells(5, 2)
b1 = a1 * F1 / (a1 - F1)
a2 = L - b1
b2 = a2 * F2 / (a2 - F2)
HH = Abs(h * b1 * b2 / a1 / a2)
Cells(6, 2) = b1: Cells(7, 2) = a2
Cells(8, 2) = b2: Cells(9, 2) = HH
End Sub

15.13. Пользуясь решением предыдущей задачи, докажите, что оптическая сила двух близко расположенных тонких линз ( l = 1 см) равна
сумме оптических сил этих линз D1 + D2 . Выполните проверку для случаев: 1) обе линзы собирающие; 2) обе линзы рассеивающие; 3) одна
линза собирающая, а другая рассеивающая.

125

Майер Р. В. Решение физических задач в электронных таблицах Excel
15. Явления
геометрической оптики

Содержание

17. Физика микромира

16. ВОЛНОВЫЕ ЯВЛЕНИЯ
16.1. Два когерентных источника с мощностью P1 и P2 излучают

гармонические звуковые волны с длиной λ синфазно (рис. 16.1). Источники отстоят друг от друга на расстоянии d и удалены от оси Oy на
расстояние L . Рассчитайте распределение интенсивности I ( y ) вдоль
оси Oy .

Рис. 16.1. Интерференция от двух точечных источников
Результирующая интенсивность в точке наблюдения A ( y ) :

I = I1 + I 2 + 2 I1I 2 cos(2π ⋅ ∆x / λ ) ,
где I1 и I 2 – интенсивности волн в точке A ( y ) , ∆x = x2 − x1 – их разность хода, x1 = ( L2 + ( y − d / 2) 2 )0,5 и x2 = ( L2 + ( y + d / 2) 2 )0,5 – расстояния, проходимые волнами. Так как вся излучаемая источником
энергия равномерно распределяется по сферической волновой поверхности, то I1 = P1 /(4π ⋅ x12 ) и I 2 = P2 /(4π ⋅ x22 ) . Для решения задачи следует создать такую же таблицу, как на рис. 16.2 и построить график
I ( y ) . Изменяя длину волны, мощности источников, расстояние между
ними и до экрана, можно изучить зависимость интерференционного
распределения интенсивности вдоль оси Oy от перечисленных
величин.

126

Майер Р. В. Решение физических задач в электронных таблицах Excel

Рис. 16.2. Распределение интенсивности при интерференции
16.2. Два источника периодически и одновременно излучают цуги
волн с равными частотами f = 1 / T и длительностью τ к , которые интерферируют. Изучите зависимость интенсивности в точке наблюдения
от разности хода x и длины когерентности Lк = υ ⋅ τ к .
Результирующие колебания находятся как сумма складываемых
колебаний. Средняя интенсивность пропорциональна интегралу от
квадрата смещения частиц среды в точке наблюдения за время
∆τ ' = τ 2 '−τ1 ' , много большее периода испускания цугов, деленное на

∆τ ' . Длина когерентности Lк равна произведению длительности цуга (то есть времени когерентности τ к ) на скорость волны υ . Если
число длин волн в цуге обозначить через n , то получим следующее:

 A sin(ω ⋅τ ),
0,

 A sin(ω ⋅ (τ + ∆τ )),
ξ 2 (τ ) = 
0,


0 ≤ τ ≤ nT ,
τ < 0, τ > nT ;
∆τ ≤ τ ≤ nT + ∆τ ,
τ < ∆τ , τ > nT + ∆τ ;
где ∆τ – время отставания второй волны от первой, x = υ∆τ – их

ξ1 (τ ) = 

если
если
если
если

разность хода. При наложении волн смещения ξ1 и ξ 2 складываются;

127

Майер Р. В. Решение физических задач в электронных таблицах Excel

чтобы получить интенсивность, результат следует возвести в
квадрат, проинтегрировать и разделить на время ∆τ ' = τ 2 '−τ1 ' :
τ2'

k
I=
(ξ1 + ξ 2 ) 2 dτ .
τ 2 '−τ 1 '

∫

τ1 '

Алгоритм построения графика I = I (x ) зависимости интенсивности при интерференции цугов волн от разности хода состоит в
следующем: 1) задается разность хода x ; 2) вычисляется средняя
интенсивность I , для этого складываются колебания, создаваемые
источниками в точке наблюдения при заданной разности хода x , результат возводится в квадрат и интегрируется по времени; 3) на
координатной плоскости “интенсивность–разность хода” ставится
точка; 4) разность хода x увеличивается на небольшую величину, после чего компьютер снова повторяет операции 2, 3, 4. Численное
интегрирование
осуществляется
методом
прямоугольников
(рис. 16.3.1).
Программа 16.1
Private Sub CommandButton1_Click()
dt = 0.002: pi = 3.1415926: n = 6
For i = 1 To 300
tau = i * 6 * pi / 150: t = 0: S = 0
While t < 40
t = t + dt
If (t > 0) And (t < 2 * pi * n) Then x1 = Sin(t) Else x1 = 0
If (t + tau > 0) And (t + tau < 2 * pi * n) Then x2 = Sin(t + tau)
Else x2 = 0
S = S + (x1 + x2) * (x1 + x2) * dt
Wend: Cells(i, 1) = tau: Cells(i, 2) = S: Next
End Sub

Рис. 16.3. Расчет интенсивности при интерференции цугов

128

Майер Р. В. Решение физических задач в электронных таблицах Excel

Используется программа 16.1. Получающиеся графики зависимости интенсивности I от разности хода x для цугов, длина Lк которых равна четырем и шести длинам волн, изображены на рис. 16.3.2.
Из них видно, что число минимумов интенсивности, наблюдаемых при
интерференции цугов, равно удвоенному числу длин волн в цуге (так
как x может быть и положительным, и отрицательным).
16.3. На прямой край непрозрачной пластины нормально падает
плоская звуковая волна с длиной λ = 0,2 м. За пластиной на расстоянии

d = 5 м установлена направляющая Oy ' , по которой перемещается датчик (рис. 16.4.1). Постройте спирали Корню для различных положений
датчика.

Рис. 16.4. Дифракция на краю пластины. Построение спирали Корню
Мысленно разрежем волновую поверхность на N элементарных
полосок шириной ∆x , параллельных краю пластины; каждая из них является источником Si элементарной волны (рис. 16.4.1). Рассчитаем
распределение интенсивности вдоль оси y . Результат сложения всех
волн в точке наблюдения M ( y ) можно представить с помощью векторной диаграммы, которую удобно построить на комплексной плос-

r

кости (рис. 16.4.2). От точки O отложим вектор a1 , длина которого
пропорциональна амплитуде, а угол

ϕ1 , образованный с осью дейст-

вительных чисел, равен фазе колебаний, создаваемых в точке наблю-

r

r

дения источником S1 . От конца вектора a1 отложим вектор a2 , соответствующий колебаниям, создаваемым в точке наблюдения источником S 2 и т. д. Из принципа суперпозиции следует, что комплекс

129

Майер Р. В. Решение физических задач в электронных таблицах Excel

& = a& + a& + ... + a& . Его
амплитуды результирующих колебаний равен: A
1
2
N
действительная и мнимая части:
N

ARE = ∑ ai cos ϕi ,
i =1

N

AIM = ∑ ai sin ϕi .
i =1

Амплитуда результирующих колебаний в точке наблюдения M ( y )
равна

2
2
A = ARE
+ AIM
, интенсивность пропорциональна квадрату

2
2
+ AIM
).
амплитуды: I = kA2 = k ( ARE

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

r

ной плоскости последовательно откладываются вектора ai . На рис.
16.4.3 и 16.5.1 представлена спираль Корню для точки O ′ , расположенной за прямым краем непрозрачной пластины ( y = 0 ), на который
падает плоская волна. На рис. 16.5.2 показана спираль Корню, получающаяся для точки наблюдения M с координатой y = 1 м.
Программа 16.2
Private Sub CommandButton1_Click()
pi = 3.1415: dx = 0.01: d = 5: lambda = 0.2: y = 0.8
For i = 1 To 6000
x = i * dx: l = Sqr(d * d + (x - y) * (x - y))
Re = Re + Cos(2 * pi * l / lambda) * d / (l * l)
Im = Im + Sin(2 * pi * l / lambda) * d / (l * l)
Cells(i, 1) = Re: Cells(i, 2) = Im
Next
End Sub

Рис. 16.5. Спираль Корню для края пластины (1, 2) и щели (3)
16.4. На щель шириной b = 3 см нормально падает плоская волна
длиной λ = 0,2 см (рис. 16.6.1). За щелью на расстоянии d = 5 см нахо-

130

Майер Р. В. Решение физических задач в электронных таблицах Excel

дится экран. Постройте спираль Корню для произвольной точки M ( y )
экрана.
Задача решается аналогично, необходимо только уменьшить
количество итераций в цикле так, чтобы захватить часть волновой
поверхности, проходящей через щель. Используется программа 16.3,
результаты для y = 1 см представлены на рис. 16.5.3.
Программа 16.3
Private Sub CommandButton1_Click()
pi = 3.1415: dx = 0.01: d = 5: lambda = 0.2: y = 1
For i = 1 To 300
x = i * dx: l = Sqr(d * d + (x - y) * (x - y))
Re = Re + Cos(2 * pi * l / lambda) * d / (l * l)
Im = Im + Sin(2 * pi * l / lambda) * d / (l * l)
Cells(i, 1) = Re: Cells(i, 2) = Im
Next
End Sub

16.5. Плоская волна нормально падает на пластину с прямым краем (рис. 16.4.1). Рассчитаем дифракционную картину, получающуюся в
плоскости, отстоящей от пластины на расстоянии d = 5 м, если длина
волны λ = 0,2 м.

Рис. 16.6. Дифракция на краю пластины и на щели
Волновую поверхность разрежем на N узких полосок шириной
∆x . Будем последовательно перебирать точки оси y с шагом ∆y и
для каждой точки вычислять результирующую амплитуду колебаний.
Фаза колебаний в точке M ( y ) , создаваемых элементарным источником с координатой xi , зависит от длины пути li = (d 2 + ( y − xi ) 2 )0,5 и
равна

r

ϕi = 2π (li / λ ) . Длина каждого вектора ai пропорциональна

cos α i = d / li ( α i – угол между нормалью к волновой поверхности и на-

131

Майер Р. В. Решение физических задач в электронных таблицах Excel

правлением на точку наблюдения) и обратно пропорциональна длине
пути li . Колебаниям, создаваемым i -й полоской, соответствует число a&i = (cos 2π (li / λ ) + j sin 2π (li / λ )) cos α i / li . Результирующие колебания в M ( y ) на комплексной плоскости изображаются вектором:
N

A& = ∑ (cos 2π (li / λ ) + j sin 2π (li / λ )) cos α i / li = ARE + jAIM .
i =1

Интенсивность волны пропорциональна квадрату модуля амплитуды: I = kA = k ( ARE + AIM ) . Программа 16.4, рассчитывающая
2

2

2

дифракционную картину, должна содержать цикл, в котором перебираются точки оси y с шагом ∆y , для каждой отдельно суммируются
действительная и мнимая части чисел a&i , определяется модуль амплитуды A( y ) и интенсивность I ( y ) . Результаты моделирования
представлены на рис. 16.6.2. По таблице можно определить координаты y1 , y2 , y3 , … максимумов дифракционной картины, координаты
минимумов, значения интенсивностей в минимумах, максимумах
и т. д.
Программа 16.4
Private Sub CommandButton1_Click()
pi = 3.1415: dx = 0.005: d = 5: lambda =
For j = 1 To 200
y = j * 0.05 - 0.5: Re = 0: Im = 0
For i = 1 To 4000
x = i * dx: l = Sqr(d * d + (x - y) * (x
Re = Re + Cos(2 * pi * l / lambda) * d /
Im = Im + Sin(2 * pi * l / lambda) * d /
S = Re * Re + Im * Im: Next
Cells(j, 1) = y: Cells(j, 2) = S: Next
End Sub

0.2: y = 1

- y))
(l * l)
(l * l)

16.6. Плоская волна нормально падает на пластину с узкой и длинной щелью (рис. 16.6.1). Рассчитаем дифракционную картину, получающуюся в плоскости, отстоящей от пластины на расстоянии d = 5 cм. Ширина щели b = 10 см, длина волны λ = 0,5 см.
Задача решается тем же способом (программа 16.5). Волновая
поверхность, ограниченная щелью, разрезается на N узких полосок
шириной ∆x = b / N и вычисляется сумма всех комплексных амплитуд,
создаваемых в точке наблюдения M ( y ) каждой такой полоской. Перебирая в цикле точки оси y с заданным шагом, программа строит
132

Майер Р. В. Решение физических задач в электронных таблицах Excel

график зависимости I ( y ) (рис. 16.6.3). В аналогичных опытах с монохроматическим светом на экране получается система светлых и
темных полос, параллельных щели. При уменьшении длины волны ширина максимумов становится меньше, их количество больше. От
этой задачи легко перейти к задаче о дифракции на круглом отверстии, при которой дифракционные максимумы образуют систему
концентрических колец.
Программа 16.5
Private Sub CommandButton1_Click()
pi = 3.1415: dx = 0.005: d = 5: lambda =
For j = 1 To 400
y = j * 0.05 - 2.5: Re = 0: Im = 0
For i = 1 To 2000
x = i * dx: l = Sqr(d * d + (x - y) * (x
Re = Re + Cos(2 * pi * l / lambda) * d /
Im = Im + Sin(2 * pi * l / lambda) * d /
S = Re * Re + Im * Im: Next
Cells(j, 1) = y: Cells(j, 2) = S / 1000:
End Sub

0.5

- y))
(l * l)
(l * l)
Next

16.7. Левый конец струны совершил одно колебание и остановился. Рассчитайте смещение ξ (x) различных точек струны в произвольный
момент времени τ ' . Правый конец нити закреплен.
Запишем уравнение одномерной волны в конечных разностях:
2
∂ 2ξ
2∂ ξ
=υ
,
∂τ 2
∂x 2

ξ it +1

=

2ξ it

− ξ it −1

+υ

2

ξ it−1 − 2ξ it + ξ it+1
∆x

2

∆τ 2 ,

где t – целочисленная переменная, которой соответствуют моменты τ = t ⋅ ∆τ . Используется программа 16.6. При запуске она считывает значения t ' из ячейки E1, вычисляет смещение ξ (x) узлов сетки в
момент t '⋅∆τ и создает таблицу, на основе которой строится график
ξ (x) . После этого t ' увеличивается на 100. При повторном запуске
макрос рассчитывает ξ (x) в момент (t '+100)∆τ и строит новый график (рис. 16.7.1). Если левый конец струны резко сместить из положения ξ = 0 , а затем через некоторое время вернуть обратно, то по
струне побежит волна, показанная на рис. 16.7.2. Смещение левого
конца задается так: If t < 8 Then xi1(1) = 20 Else xi1(1) = 0.
Программа 16.6
Private Sub CommandButton1_Click()
Dim xi(100) As Single: Dim xi1(100) As Single

133

Майер Р. В. Решение физических задач в электронных таблицах Excel

Dim xi2(100) As Single
dx = 1: dt = 0.02: v = 8: Vremya = Cells(1, 5)
For i = 1 To 100: Cells(i, 1) = i * dx: Next i
For k = 1 To Vremya: t = t + dt
If 2 * t < 6.28 Then xi1(1) = 20 * Sin(2 * t) Else xi1(1) = 0
For i = 2 To 99
xi2(i) = 2 * xi1(i) - xi(i) + v * (xi1(i - 1) - 2 * xi1(i)
+ xi1(i + 1)) * dt * dt / dx / dx
Next i
For i = 2 To 99: xi(i) = xi1(i): xi1(i) = xi2(i)
Next i: xi(100) = 0: Next k
For i = 1 To 100: Cells(i, 2) = xi(i): Next i
Cells(1, 5) = Cells(1, 5) + 100
End Sub

Рис. 16.7. Результат решения волнового уравнения
16.8. Левый конец струны совершает одно колебание, и по ней бежит импульс. Струна неоднородная: скорость распространения волны по
ее левой половине υ1 не равна скорости распространения υ 2 по правой
половине. Рассчитайте смещение ξ (x) различных точек струны в последовательные моменты времени, промоделируйте отражение одномерной волны от границы раздела двух сред и ее прохождение во вторую
среду, когда υ1 > υ 2 и υ1 < υ 2 .
Задача решается тем же способом. Используется программа
16.7. Результат – на рис. 16.8: волна слева подходит к “границе раздела сред”, частично отражается и частично проходит во вторую
среду.
Программа 16.7
Private Sub CommandButton1_Click()
Dim xi(200) As Single: Dim xi1(200) As Single
Dim xi2(200) As Single
dx = 1: dt = 0.02: Vremya = Cells(1, 5)
For i = 1 To 200: Cells(i, 1) = i * dx: Next i
For k = 1 To Vremya: t = t + dt
If 2 * t < 6.28 Then xi1(1) = 20 * Sin(2 * t) Else xi1(1) = 0

134

Майер Р. В. Решение физических задач в электронных таблицах Excel

For i = 2 To 199
If i > 100 Then v = 20 Else v = 8
xi2(i) = 2 * xi1(i) - xi(i) + v * (xi1(i - 1) - 2 * xi1(i) + xi1(i
+ 1)) * dt * dt / dx / dx
Next i
For i = 2 To 199: xi(i) = xi1(i): xi1(i) = xi2(i)
Next i: xi(200) = 0: Next k
For i = 1 To 200: Cells(i, 2) = xi(i): Next i
Cells(1, 5) = Cells(1, 5) + 100
End Sub

Рис. 16.8. Отражение и прохождение волны через границу двух сред
16.9. Левый конец струны совершает гармонические колебания,
правый конец закреплен. Образующаяся волна доходит до закрепленного конца, отражается и интерферирует с падающей. Промоделируйте
образование стоячей волны.

135

Майер Р. В. Решение физических задач в электронных таблицах Excel

Задача решается аналогично (программа 16.8). Сначала волна
распространяется вправо, затем отражается (рис. 16.9.1–16.9.3). В
области наложения падающей и отраженной волн возникает стоячая
волна (рис. 16.9.4–16.9.6). Некоторые точки струны колеблются
с максимальной амплитудой (пучности), а некоторые – не колеблются
вообще (узлы).

Рис 16.9. Интерференция падающей и отраженной волн
Программа 16.8
Private Sub CommandButton1_Click()
Dim xi(100) As Single: Dim xi1(100) As Single
Dim xi2(100) As Single
dx = 1: dt = 0.02: v = 8: Vremya = Cells(1, 5)
For i = 1 To 100: Cells(i, 1) = i * dx: Next i
For k = 1 To Vremya: t = t + dt: xi1(1) = 20 * Sin(t)
For i = 2 To 99
xi2(i) = 2 * xi1(i) - xi(i) + v * (xi1(i - 1) - 2 * xi1(i) +
xi1(i + 1)) * dt * dt / dx / dx: Next i
For i = 1 To 99: xi(i) = xi1(i): xi1(i) = xi2(i)
Next i: xi(100) = 0: Next k
For i = 1 To 100: Cells(i, 2) = xi(i): Next i
Cells(1, 5) = Cells(1, 5) + 30
End Sub

136

Майер Р.В. Решение физических задач в электронных таблицах Excel
16. Волновые явления

Содержание

18. Другие задачи по физике

17. ФИЗИКА МИКРОМИРА
17.1. Постройте график зависимости спектральной светимости r
абсолютно черного тела от частоты ν излучения при различных температурах T .
Зависимость спектральной светимости r абсолютно черного
тела от частоты излучения ν и температуры T выражается формулой Планка:

2π ⋅ hv 3

r ( v, T ) =

.

c 2 (exp(hv / kT ) − 1)
От абсолютной температуры T и частоты излучения v перейдем к
параметрам
1/ 3

k  2π ⋅ h 
T '= T ⋅

h  c2 

1/ 3

 2π ⋅ h 
и v ' = v

 c2 

.

Тогда формула Планка приобретет простой вид:

r (v' , T ' ) = v'3 /(exp(v' / T ' ) − 1) .
Используемая программа 17.1 строит графики r = r (v' ) при T ' = 500,
1000, 1500, 2000, 2500 (рис. 17.1). Видно, что при увеличении температуры спектральная светимость черного тела на всех частотах
возрастает, максимум смещается в сторону высоких частот.

Рис. 17.1. Графики зависимостей r = r (v' ) при различных T '

137

Майер Р. В. Решение физических задач в электронных таблицах Excel

Программа 17.1
Private Sub CommandButton1_Click()
For i = 1 To 5
T = 500 * i
For j = 1 To 1000
nu = 20 * j: F = nu * nu * nu / (Exp(nu / T) - 1)
Cells(j, 1) = nu: Cells(j, 1 + i) = F
Next: Next
End Sub

17.2. Методом численного интегрирования определите интегральную светимость R абсолютно черного тела для данной температуры.
Подтвердите, что R прямо пропорциональна четвертой степени T :

R = σT 4 (закон Стефана-Больцмана).
Для расчета интегральной светимости черного тела R можно
применить метод прямоугольников:
∞

ν '3i ∆ν '
R = ∫ r (ν ' )dν ' ≈∑ r (ν 'i )∆ν ' = ∑
,
exp(
ν
'
/
T
'
)
−
1
i
i =1
i =1
0
n

n

где ν 'i = i∆ν ' . Используется программа 17.2, перед ее запуском значение T ' следует записать в ячейку D1. После нажатия на кнопку запуска программа создает таблицу из двух столбцов ν ' и r , содержащую
20000 строк, а в ячейке E1 появляется значение интегральной светимости R . Изменяя T ' , заполняют таблицу 17.1, из которой следует, что σ = R / T 4 = const , то есть R увеличивается пропорционально

T 4 . Все величины измеряются в условных единицах.
Таблица 17.1

Программа 17.2
Private Sub CommandButton1_Click()
T = Cells(1, 4): h = 1: Max = 0
For nu = 1 To 20000
F = nu * nu * nu / (Exp(nu / T) - 1)
Cells(nu, 1) = nu: Cells(nu, 2) = F

138

Майер Р. В. Решение физических задач в электронных таблицах Excel

R = R + F * h
If F > Max Then wm = nu: Max = F
Next: Cells(1, 5) = R: Cells(2, 5) = 1 / wm
End Sub

17.3. Найдите длину волны λm , соответствующую максимуму спектральной светимости r = r (v ' ) абсолютно черного тела, имеющего температуру T ' . Подтвердите, что эта длина волны λm обратно пропорциональна его температуре: λm = b / T (закон смещения Вина).
Используется та же программа 17.2. Температуру T ' записывают в ячейку D1, после запуска макроса получающиеся значения
λm = c /ν (в условных единицах) появляются в ячейке E2. С помощью
программы можно рассчитать b = λmT и заполнить таблицу 17.1. Из
нее видно, что по мере увеличения T длина волны λm , соответствующая максимуму спектральной светимости, уменьшается так, что
величина b = λmT остается постоянной.
17.4. Изотоп A распадается с периодом полураспада 2,9 мин, образующийся при этом изотоп B имеет период полураспада 24 мин. Постройте графики зависимости количества атомов A и B от времени, график зависимости активности смеси от времени. Определите момент
времени, когда количество изотопа B будет максимально.
Запишем закон радиоактивного распада:

N A = N A0 e − λ A t , N B = N B 0 e − λ B t , λ A = ln 2 / TA , λB = ln 2 / TB .
Пусть в момент t имеется N A атомов А и N B атомов B. Через время dt число нераспавшихся ядер А равно N A (t + dt ) = N A (t ) exp(−λ ⋅ dt ) .
При этом распадается

dN ' A = N A (t ) − N A (t + dt ) = N A (1 − exp(−λ A ⋅ dt )) ,
dN ' B = N B (t ) − N B (t + dt ) = N B (1 − exp(−λ B ⋅ dt ))
ядер А и B соответственно. После распада из ядер А получаются ядра В, поэтому
N A (t + dt ) = N A (t ) − dN ' A , N B (t + dt ) = N B (t ) − dN ' B + dN ' A .
Активность равна A(t ) = (dN ' A + dN ' B ) / dt . Используется программа
17.3, результаты приведены на рис. 17.2.
Программа 17.3
Private Sub CommandButton1_Click()
TA = 2.9 * 60: TB = 24 * 60: NA = 10000: dt = 10

139

Майер Р. В. Решение физических задач в электронных таблицах Excel

lambdaA = 0.693 / TA: lambdaB = 0.693 / TB
While t < 2000
t = t + dt: i = i + 1
dNA = NA * (1 - Exp(-lambdaA * dt))
NA = NA - dNA: NB = NB + dNA - dNB
dNB = NB * (1 - Exp(-lambdaB * dt))
Cells(i, 1) = t: Cells(i, 2) = NA: Cells(i, 3) = NB
Cells(i, 4) = (dNA + dNB) * 200 / dt
Wend
End Sub

Рис. 17.2. К решению задачи 17.4 о радиоактивном распаде
17.5. Ядро атома магния состоит из 23 нуклонов из которых 12 протонов. Масса атома равна M = 22,99414 а.е.м. Напишите программу,
которая вычисляет удельную энергию связи ядра в электрон-вольтах.
Для расчета удельной энергии связи используются формулы:

Eсв = ( Zm p + ( A − Z )mn − M Я )с 2 ,

M Я = M − Zme , E уд = Eсв / A .

Если массы в а.е.м, то Eсв = 931,5( Zm p + ( A − Z ) mn − M Я ) МэВ. Следует

учесть,

что

m p = 1,00728

а.е.м.,

mn = 1,00866

а.е.м.,

me = 0,0005486 а.е.м. Используется программа 17.4.
Программа 17.4
Private Sub CommandButton1_Click()
A = Cells(1, 2): Z = Cells(2, 2): M = Cells(3, 2)

140

Майер Р. В. Решение физических задач в электронных таблицах Excel

c2 = 931.5: mp = 1.00728: mn = 1.00866
m_e = 0.0005486
M_yadra = M - Z * m_e
E = (Z * mp + (A - Z) * mn - M_yadra) * c2
Cells(4, 2) = E / A
End Sub

Рис. 17.3. Расчет удельной энергии связи в таблицах Excel
17.6. Используя решение предыдущей задачи, вычислите удельную энергию связи для ядер лития, железа, урана и т. д. Постройте график зависимости удельной энергии связи от массового числа (атомного
номера) А.
Используется программа 17.4. С ее помощью можно вычислить
удельную энергию связи для различных ядер и построить график зависимости E уд от атомного номера A (рис 17.4).

Рис. 17.4. Зависимость удельной энергии связи от массового числа
17.7. При столкновении двух ядер Я1 и Я2 происходит ядерная реакция и получаются ядра Я3 и Я4. Напишите программу в Excel, которая
вычисляла бы энергетический выход реакции.

141

Майер Р. В. Решение физических задач в электронных таблицах Excel
17. Физика микромира

Содержание

Список литературы

18. ДРУГИЕ ЗАДАЧИ ПО ФИЗИКЕ
18.1. Заяц бежит по полю вдоль прямой; его координаты изменяются так: xZ (t ) = 0,9 ⋅ τ (м), y Z (t ) = 5 + 0,4 ⋅ τ (м). Из точки O (0; 0) начинает двигаться волк со скоростью υV = 1,15 м/с, все время направленной к
зайцу. Определите траектории движения зайца и волка, время и место
их встречи.
Угол между направлением движения волка и осью Ox обозначим
через α . Запишем используемые формулы:

r = ( xZ − xV ) 2 + ( yZ − yV ) 2 , cos α = ( xZ − xV ) / r , sin α = ( yZ − yV ) / r ,
xVt +1 = xVt + υVt +1 cos α ⋅ ∆τ ,

yVt +1 = yVt + υVt +1 sin α ⋅ ∆τ .

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

Рис. 18.1. Траектории движения зайца и волка;
зависимость расстояния между ними от времени
Программа 18.1
Private Sub CommandButton1_Click()
dt = 0.2: v = 0.15: r = 0.2
While (t < 150) And (r > 0.1)
t = t + dt: i = i + 1: xz = 0.1 * t: yz = 0.05 * t + 5
r = Sqr((xz - xv) * (xz - xv) + (yz - yv) * (yz - yv))
cosa = (xz - xv) / r: sina = (yz - yv) / r
xv = xv + v * cosa * dt: yv = yv + v * sina * dt
Cells(i, 1) = t: Cells(i, 2) = xz
Cells(i, 3) = yz: Cells(i, 4) = xv
Cells(i, 5) = yv: Cells(i, 6) = r
Wend
End Sub

142

Майер Р. В. Решение физических задач в электронных таблицах Excel

18.2. Заяц бежит по окружности радиусом R со скоростью

υ Z . Из

центра окружности O начинает двигаться волк. Его скорость υV остается
постоянной и всегда направлена в сторону зайца. Определите траектории движения зайца и волка, время и место встречи.
18.3. Изучите движение шарика в поле тяжести по поверхности, задаваемой уравнением z = 0,25( x 2 + y 2 ) 2 − 0,5( x 2 + y 2 ) и представляющей замкнутую канавку вокруг оси Oz . Постройте траекторию шарика
при различных начальных координатах и скорости.
Потенциальная энергия шарика массой m , движущегося по рассматриваемой поверхности, равна:

U ( x, y ) = mg (0,25( x 2 + y 2 ) 2 − 0,5( x 2 + y 2 )) .
Горизонтальные проекции силы реакции, действующей на шарик:

Fx = −∂U / ∂x = − mg ( x 3 + y 2 x − x),

Fy = −∂U / ∂y = − mg ( y 3 + x 2 y − y ) .

Самостоятельно рассчитайте траекторию в плоскости XOY , используя решения задачи 5.8.
18.4. Из пушки вылетает снаряд с известной скоростью

υ 0 . Опре-

делите угол α , под которым необходимо произвести выстрел, чтобы
снаряд попал в цель с координатами ( x' , y ' ) . Сила сопротивления воз-

r

r

духа, действующая на снаряд, пропорциональна его скорости FC = − rυ .
Используется метод стрельбы. Запишем второй закон Ньютона

r
r
r
mg − rυ = ma в конечных разностях (см. задачу 5.1):
a tx+1 = − rυ xt / m ,

a ty+1 = − g − rυ ty / m ,

υ xt +1 = υ xt + a tx+1∆τ ,

xt +1 = xt + υ xt +1∆τ ,

υ ty+1 = υ ty + a ty+1∆τ ,

y t +1 = y t + υ ty+1∆τ .

Самостоятельно напишите программу, рассчитывающую траекторию движения снаряда при заданной начальной скорости и произвольном угле α между стволом орудия и горизонталью. Если снаряд не
поразил цель, то угол α увеличивают на величину шага ∆α и повторяют расчет, затем еще раз и еще раз до тех пор, пока снаряд не попадет в цель. Искомое значение угла в радианах выводится на экран
компьютера. Задача может иметь два решения, одно или не иметь
решения.

143

Майер Р. В. Решение физических задач в электронных таблицах Excel

18.5. Напряжение изменяется по закону u (τ ) =| 10 sin ωτ | . Разложите импульсы напряжения в ряд Фурье. Восстановите сигнал по 3, 7, 15 и
25 первым гармоникам и получите график u ' (τ ) .
По теореме Фурье периодическая функция u (τ ) с периодом T
представима в виде суммы постоянной составляющей a0 и бесконечного числа гармонических колебаний (гармоник), с частотами кратными основной частоте ω = 2π / T (частоты ω , 2ω , 3ω , …):
∞

u ' (τ ) = a0 + ∑ (an cos nωτ + bn sin nωτ ) , где
n =1

1T
2T
2T
a0 = ∫ u (τ )dτ , an = ∫ u (τ ) cos(nωτ )dτ , bn = ∫ u (τ ) sin(nωτ ) dτ .
T0
T0
T0
Используется программа 18.2, результаты разложения исходной
функции u (t ) =|10 sin(ωτ ) | в ряд Фурье и ее восстановления представлены на рис. 18.2.

Рис. 18.2. Гармонический анализ и синтез функции u (τ ) =| 10 sin(ω ⋅ τ ) |
Программа 18.2
Private
Ch_garm
w = 2 *
For i =
t = i *
U0 = U0
For n =
aa = 0:
For i =

Sub CommandButton1_Click()
= 25: Pi = 3.1415926: dt = 2 * Pi / 1000: Period = 2 * Pi
Pi / Period: 'Deffnn U(t) = Abs(10 * Sin(w * t))
1 To 1000
dt: U0 = U0 + Abs(10 * Sin(w * t)) * dt: Next i
/ Period: Cells(1, 5) = U0
1 To Ch_garm
bb = 0
1 To 1000

144

Майер Р. В. Решение физических задач в электронных таблицах Excel

t = i * dt: U = Abs(10 * Sin(w * t))
aa = aa + U * Cos(n * w * t) * dt
Next i: a = 2 * aa / Period
For i = 1 To 1000
w = 2 * Pi / Period: U = Abs(10 * Sin(w * t))
'{If t>Period/2 then U:=0 else U:=10;}
t = i * dt: bb = bb + U * Sin(n * w * t) * dt: Next i
b = 2 * bb / Period: Cells(n, 1) = a: Cells(n, 2) = b
Next n
End Sub
Private Sub CommandButton2_Click()
Ch_garm = 25: Pi = 3.1415926
dt = 2 * Pi / 1000: Period = 2 * Pi
w = 2 * Pi / Period: U0 = Cells(1, 5)
For i = 1 To 500
t = i * dt * 10: UV = U0
For n = 1 To Ch_garm
a = Cells(n, 1): b = Cells(n, 2)
UV = UV + a * Cos(n * w * t) + b * Sin(n * w * t)
Next n: Cells(i, 3) = t: Cells(i, 4) = UV: Next i
End Sub

18.6. Цепь переброшена через невесомый блок так, что с одной
стороны свисает 19 м, а с другой 20 метров цепи (рис. 18.3.1). Как будет
двигаться цепь, если ее отпустить? Постройте графики зависимости ускорения, координаты и скорости конца цепи B от времени.
Все точки цепи движутся с одинаковым по модулю ускорением,
которое можно найти из второго закона Ньютона:

(m1 + m2 )a = (m2 − m1 ) g , a =

m2 − m1
L −L
g = 2 1g.
m1 + m2
L1 + L2

Допустим, при τ = 0 координата точки B (правого конца цепи) x = 0 .
Ее скорость и координата в моменты времени ti = i ⋅ ∆τ равны:

υ t +1 = υ t + a t +1∆τ ,

x t +1 = x t + υ t +1∆τ .

Длина левой части цепи уменьшается, а правой части – увеличивается: L1t +1 = L1t − υ t +1∆τ , Lt2+1 = Lt2 + υ t +1∆τ . Так продолжается до тех пор,
пока L1 не уменьшится до 0. Ускорение цепи возрастает от 0 до

g = 9,81м/c 2 , когда она полностью соскочила с блока (рис. 18.3.2).

145

Майер Р. В. Решение физических задач в электронных таблицах Excel

Рис. 18.3. Движение цепи, переброшенной через неподвижный блок
Программа 18.3
Private Sub CommandButton1_Click()
L1 = 19: L2 = 20: dt = 0.01
While t < 7
i = i + 1: t = i * dt
a = 9.81 * (L2 - L1) / (L1 + L2)
v = v + a * dt: x = x + v * dt
If L1 > 0 Then
L1 = L1 - v * dt: L2 = L2 + v * dt
End If
If i Mod 10 = 0 Then
k = i / 10: Cells(k, 1) = t: Cells(k, 2) = a
Cells(k, 3) = v: Cells(k, 4) = x
End If
Wend
End Sub

18.7. В инерциальной системе отсчета OXYZ в точке с координатами ( x0 , y0 , z0 ) в момент времени t0 происходит событие A. Рассчитайте пространственные и временные координаты (t ' , x' , y ' , z ' ) события в
подвижной инерциальной системе отсчета O ' X ' Y ' Z ' , которая движется
относительно ИСО OXYZ со скоростью υ в направлении Ox . Найдите
квадрат радиус-вектора события в пространстве Минковского в обеих
системах отсчета.
Радиус вектор события в четырехмерном пространстве равен
пространственно-временному интервалу между событием O (0,0,0,0)
и событием A. В системе отсчета OXYZ квадрат радиус-вектора ра-

146

Майер Р. В. Решение физических задач в электронных таблицах Excel

вен: S 2 = −c 2 t 2 + x 2 + y 2 + z 2 . Найдем пространственно-временные
координаты события A в движущейся системе отсчета O ' X ' Y ' Z ' .
Запишем преобразования Лоренца:

x' = γ ( x − υ ⋅ t ) , y ' = y , z ' = z , t ' = γ (t − (υ / c 2 ) x) ,
где

β = υ / c , γ = 1 / 1 − β 2 . В системе отсчета O' X ' Y ' Z ' квадрат

радиус–вектора рассчитывается так: S '2 = −c 2t '2 + x'2 + y '2 + z '2 . Это
инвариант преобразований Лоренца, поэтому S = S ' . Используя записанные выше формулы, вычислите в таблицах Excel S 2 , x ' , y ' , z ' , t ' и

S '2 . Докажите, что S 2 ≈ S '2 .

Рис. 18.4. Расчет S 2 в различных системах отсчета.

147

Майер Р. В. Решение физических задач в электронных таблицах Excel
Содержание

18. Другие задачи по физике

Список литературы
1. Булавин, Л. А., Выгорницкий, Н. В., Лебовка, Н. И. Компьютерное
моделирование физических систем. Долгопрудный: Интеллект, 2011.
352 c.
2. Васильев, А. Н. Научные вычисления в Microsoft Excel. М.: Вильямс, 2004. 512 с.
3. Гельман, В. Я. Решение математических задач средствами
Excel: Практикум. СПб.: питер, 2003. 240 с.
4. Гулд, Х., Тобочник, Я. Компьютерное моделирование в физике. В
2 ч. Ч. 2. М.: Мир, 1990. 400 с.
5. Компьютеры, модели, вычислительный эксперимент. Введение в
информатику с позиций математического моделирования. М.: Наука,
1988. 176 c.
6. Кунин, С. Вычислительная физика. М.: Мир, 1992. 518 с.
7. Майер, Р. В. Задачи, алгоритмы, программы [Электронный ресурс]. Глазов: Глазов. гос. пед. ин-т, 2012. URL: http://maier-rv.glazov.net
8. Майер, Р. В. Как стать компьютерным гением, или Книга об информационных системах и технологиях. Глазов: Глазов. гос. пед. ин-т,
2008. 204 c.
9. Майер, Р. В. Компьютерное моделирование физических явлений. Глазов: Глазов. гос. пед. ин-т, 2009. 112 с.
10. Майер, Р. В. Компьютерное моделирование: учебно-методическое пособие для студентов педагогических вузов [Электронное учебное
издание на компакт-диске]. – Глазов: Глазов. гос. пед. ин-т, 2015. – 24,3 Мб.
URL: http://maier-rv.glazov.net
11. Майер, Р. В. Решение физических задач с помощью электронных таблиц MS Excel // International Journal of Open Information Technologies. 2014, vol. 2, no. 9. P. 18–23.
12. Могилев, А. В., Пак, Н. И., Хеннер, Е. К. Информатика: учеб.
пособие для студ. пед. вузов. М.: Академия, 2003. 816 с.
13. Поттер, Д. Вычислительные методы в физике. М.: Мир, 1975.
392 с.

148

Майер Р. В. Решение физических задач в электронных таблицах Excel

14. Ращиков, В. И., Рошаль, А. С. Численные методы решения физических задач: учеб. пособие. СПб.: Лань, 2005. 208 с.
15. Рябенький, В. С. Введение в вычислительную математику:
учеб. пособие. М.: Физматлит, 2000. 296 с.
16. Самарский, А. А., Вабищевич, П. Н., Самарская, Е. А. Задачи и
упражнения по численным методам: учеб. пособие. М.: Эдиторал УРСС,
2000. 208 с.
17. Самарский, А. А., Михайлов, А. П. Математическое моделирование: Идеи. Методы. Примеры. М.: Физматлит, 2005. 320 с.
18. Советов, Б. Я., Яковлев, С. А. Моделирование систем: учебник
для вузов. М.: Высш. шк., 2001. 343 с.
19. Старовиков, М. И. Введение в экспериментальную физику:
учеб. пособие. СПб.: Лань, 2008. 240 с.
20. Угринович, Н. Д. Исследование информационных моделей.
Элективный курс: учеб. пособие. М.: БИНОМ. Лаборатория знаний, 2004.
183 с.
21. Уокенбах, Д. Профессиональное программирование на VBA в
Excel 2002. М., СПб., Киев, 2003. 781 с.
22. Федоренко, Р. П. Введение в вычислительную физику: учеб. пособие для вузов. М.: Изд-во Моск. физ.-техн. ин-та, 1994. 528 с.
23. Эксперимент на дисплее: Первые шаги вычислительной физики. М.: Наука, 1989. 175 с.
24. Giordano, N. J. Computational Physics. New Jersey, Prentice Hall,
1997. 419 p.
25. Phillipson, P. E., Schuster, P. Modeling by Nonlinear Differential
Equations: Dissipative and Conservative Processes. World Scientific Publishing, 2009. 225 p.
26. Woolfson, M. M., Pert, G. J. An Introduction to Computer Simulation. Oxford University Press, 1999. 311 p.

149

Майер Р. В. Решение физических задач в электронных таблицах Excel
Титульный экран

СОДЕРЖАНИЕ
Введение
1. Кинематика точки: одномерное движение
2. Кинематика точки: двумерное движение
3. Динамика: одномерное движение
4. Механические колебания
5. Динамика: двумерное движение в силовом поле
6. Задача двух тел
7. Динамика твердого тела
8. Аналитическая механика
9. Динамический хаос
10. Молекулярная физика
11. Термодинамика
12. Электрическое и магнитное поле
13. Электрические цепи постоянного тока
14. Переменный ток. Нестационарные процессы
в электрических цепях
15. Явления геометрической оптики
16. Волновые явления
17. Физика микромира
18. Другие задачи по физике
Список литературы

150