Практикум по матлабу. практикум по матлабу. Физических процессов с использованием

• Название: Физических процессов с использованием
• Анкор: Практикум по матлабу
• Дата: 21.06.2021
• Размер: 1.13 Mb.
• Имя файла: практикум по матлабу.pdf
• Тип: Учебное пособие; #219898
• страница: 5 из 17

Alpha=1; dt=0.025;
x1=2.5; y1=0;
vx1=0; vy1=0.25;
Emax=-1e18;
Emin=1e18;
r=sqrt(x1^2+y1^2);
% Расстояние до центра
E1=(vx1^2+vy1^2)/4 ...
- Alpha/r;
% Полная энергия (почему так?)
hl=line(x1,y1);
% Задание дескриптора линии
a=Alpha/(2*abs(E1));
% Большая полуось
b=r*vy1/sqrt(2*abs(E1));
% Малая полуось
axis([-2*a 2*a -1.2*b 1.2*b]);
% Задание масштаба
set(hl,’EraseMode’,’none’,’Color’,’g’);
ha=gca;
% Задание разметки осей
set(ha,’XTick’,[-a 0 a],’YTick’,[-b 0 b]);
grid on;
% Введение сетки
pause
% Пауза для обязательного вывода
while 1
% Бесконечный цикл
36

x2=x1+vx1*dt;
% Новые координаты
y2=y1+vy1*dt;
%
r=sqrt(x2^2+y2^2);
% Расстояние до центра
A = Alpha/r^2;
% Модуль ускорения
ax=-A*x2/r;
% Компоненты вектора
ay=-A*y2/r;
% ускорения
vx2=vx1+ax*dt;
% Новые скорости
vy2=vy1+ay*dt;
%
E2=((vx2+vx1)^2+(vy2+vy1)^2)/8 ...
-Alpha/r;
% Полная энергия
if E2
> Emax
% Границы изменения
Emax=E2;
% полной энергии
end;
if E2
< Emin
Emin=E2;
end;
% Обновление данных на графике
set(hl,’XData’,[x1 x2],’YData’,[y1 y2]);
% Переприсвоение
x1=x2;
y1=y2;
vx1=vx2;
vy1=vy2;
end;
% конец цикла ’while 1’
Работа программы понятна из приведенных в ее тексте комментариев. Следует только сделать несколько замечаний. Поскольку используется бесконечный цикл,
выполнение программы можно прервать клавишей Ctr+C. Следует помнить, что в
MATLAB все данные хранятся в рабочей области до тех пор, пока Вы не выйдете из командного окна. Поэтому после первого прогона программы можно остано- вить вывод командой Ctr+C, поменять требуемые параметры (например, шаг по времени), задать режим
hold on(рисовать поверх), и снова запустить программу.
График будет выведен на тот же рисунок
17
Прежде чем выполнять задания, обсудим некоторые особенности начальной программы расчета траектории планеты. Прежде чем воспользоваться возможно-
17
Можно, конечно, уже сейчас, используя опыт предыдущей задачи (Биения, см. п.
2.3
) и
разработанные m-файлы, создать графический интерфейс и работать далее с ним. Но можно
это сделать и позже или не делать вовсе. Это зависит от того, что проще – разрабатывать
интерфейс или использовать примитивные приемы остановки и перезапуска программы.
37

стями анимации (см. п.
2.4.2
), необходимо приготовить дескриптор той линии или линий, на который потом можно будет ссылаться при смене координат. Кроме то- го, задается масштаб по осям и для проведения линий сетки только в определенных местах графика (а не по умолчанию) явно определяются места разметки осей (опе- ратор
set(ha,’XTick’,[-a 0 a],’YTick’,[-b 0 b]) и следующий за ним оператор
grid on.
Задание 1.
Исследуйте влияние величины шага
∆t на точность расчета траекто- рии. При слишком большом значении
∆t возникает прецессия орбиты, свя- занная с неточностью расчета.
Для оценки точности вычислений выведите на экран значение полной энергии в зависимости от координаты x. Это лучше сделать, разбив область графиче- ского вывода на две подобласти – для построения орбиты и для построения энергии. Это можно сделать с помощью функции
subplot или явным зада- нием областей вывода графиков.
Пометьте на траектории другим цветом точки наибольшего и наименьшего удаления от центра.
Задание 2.
Исследуйте влияние на траекторию частицы возмущения вида δU
=
β/r
2
Предусмотрите вывод на экран «центра тяжести» витка траектории (виток –
участок от одной точки траектории r
= r
min
до следующей).
Если представить, что движущаяся частица – электрон в атоме, то интересу- ющий нас «центр тяжести» определяет среднее по времени значение диполь- ного момента. Речь идет, таким образом, о центре тяжести отрезка кривой с учетом времени, проводимого частицей на каждом его участке.
Для исследования точности расчетов можно в некоторый момент изменить направление движения частицы на строго противоположное. В таком случае частица должна будет двигаться по уже «проложенной» траектории в обрат- ном направлении. Проследить, насколько это выполняется, можно, изменив цвет следа.
5.2.
Влияние малого возмущения
Интересно исследовать, как влияет на движение частицы малая добавка δU к по- тенциальной энергии U
0
. Если добавка δU зависит только от r
= |r|, то поле
U
= U
0
+ δU остается центральным, а траектория перестает быть замкнутой и ее
38

можно представлять как прецессирующий эллипс, постепенно заполняющий коль- цо r
min
≤ r ≤ r
max
. Если же мы стартуем от поля, уже отличного от кулоновско- го, то та же добавка δU приведет к изменению скорости прецессии орбиты, но не вызовет качественного изменения траектории.
Таким образом, кулоновское поле в сравнении с другими центральными полями обладает той особенностью, что траектория финитного движения в нем оказывает- ся гораздо более чувствительной к «возмущению» δU .
Подобная особенность кулоновского поля становится более яркой, если в ка- честве возмущения δU выбрано нецентральное поле. Например, добавление к ку- лоновскому очень слабого однородного поля, не параллельного плоскости орбиты,
δU
= −fr, приводит не только к существенным деформациям орбиты, но и к зна- чительным поворотам ее плоскости.
То же слабое однородное поле δU исказит траекторию в другом исходном цен- тральном поле (скажем, U
0
= −α/r + β/r
2
) гораздо меньше
18
. В частности, не будет значительного поворота плоскости орбиты!
Такая повышенная чувствительность орбиты в кулоновском поле к возмуще- ниям связана в конечном счете с тем, что орбита в кулоновском поле при любых значениях M и E <
0 замкнутая.
Чтобы понять причину таких различий, проследим за изменением момента им- пульса. Скорость его изменения равна моменту силы, действующей на частицу
(причем силу притяжения к центру можно исключить): ˙
M = [rf]. Приращение вектора
M за один период движения в кулоновском поле мало, но за несколько пе- риодов подобные приращения накапливаются и приводят к большому его измене- нию. В результате сильно изменяются положение плоскости орбиты (определяемое направлением
M) и малая полуось эллипса (определяемая его величиной). Если же исходное поле U
0
(r) заметно отличается от кулоновского, то траектория представ- ляет собой "прецессирующий эллипс", и приращения
M за несколько периодов взаимно компенсируются, так что плоскость орбиты лишь слегка покачивается.
Задание 3.
Исследуйте влияние на орбиту частицы в кулоновском поле слабого однородного поля δU
= −fr. Если сила лежит в плоскости орбиты, то ор- бита искажается, постепенно превращаясь в вырожденный эллипс – отрезок.
При этом происходит нарушение точности счета. Аккуратный расчет требо- вал бы в этом случае уменьшения шага и перехода к другой, более сложной вычислительной схеме.
18
Подразумевается, что возмущение δU намного меньше, чем отличие поля U
0
от кулонов-
ского β/r
2
.
39

Лучше обойти эти трудности, исследуя движение в пространстве и под дей- ствием возмущающей силы
f, не параллельной плоскости траектории. В этом случае орбита искажается и поворачивается, но в отрезок не вырождается.
Для наблюдения орбиты, если она оказывается пространственной кривой, удоб- но изображать сразу три ее проекции на плоскости (X,Y), (Y,Z) и (X,Z),
разместив их в трех расположенных рядом графических окнах. (Необходи- мо согласовать направление осей координат в разных окнах друг с другом.)
Можно также с помощью свойства кривой ’ZData’ строить пространствен- ное изображение траектории (подробнее о построении трехмерной кривой см.
Дополнение, п.
8.3
).
В этом случае становится ясно, что анализ результатов – задача не менее сложная, чем построение модели. (Аналогично обстоит дело и в обычном,
«железном» эксперименте.)
Задание 4.
Повторите исследование влияния возмущения δU
= −fr на финит- ное движение, если исходное поле сильно отличалось от кулоновского. Мож- но положить, например, U
0
= −α/r + β/r
2
5.3.
Движение двух частиц
Движение нескольких взаимодействующих друг с другом тел лишь в исключитель- ных случаях можно рассчитать аналитически. В то же время численный расчет ока- зывается совсем немногим более сложен, чем для одного тела, т.е. вполне доступен.
При исследовании возможных движений в этом случае наиболее сложной частью задачи оказывается представление результатов – перебор и описание различных видов движения, понимание, какие из движений нужно считать качественно раз- личными. Разумеется, такая постановка вопроса не предполагает однозначного от- вета, да и исследование в рамках данного практикума не может быть полным. Это скорее иллюстрация наличия практически неограниченного числа возможностей.
Задание 5.
Изобразите на экране движение одновременно двух точек в поле
U
= −α/r.
Получите траекторию одной точки с точки зрения другой (т.е. траекторию
Венеры, например, с точки зрения земного наблюдателя).
Введите взаимодействие частиц друг с другом: U
12
= β/|r
1
− r
2
|.
"Выключив"поле U
(r) и сохраняя лишь взаимодействие двух тел, получите их траектории.
40

6.
Случайные блуждания и диффузия
Хаотическое движение броуновских частиц в жидкости или газе представляет со- бой пример случайных блужданий. Теория броуновского движения была построена
А.Эйнштейном и М.Смолуховским в 1905 -1906 гг.
Задача о случайных блужданиях является одной из широко исследуемых задач теории вероятностей и находит множество других приложений.
6.1.
Закономерности случайных блужданий
Закономерности случайных блужданий можно понять, используя простую модель,
которая легко реализуется с помощью компьютера.
N частиц (которые в начальный момент для удобства наблюдений распределе- ны на оси y) смещаются последовательными шагами
∆x вдоль оси x. Каждый шаг каждой частицы выбирается случайным и независимым от других шагов. Однако распределение вероятностей при выборе любого шага одно и то же. Примем, что смещения в противоположные стороны равновероятны. Это значит, что среднее значение смещения
∆x = 0.
(1)
Смысл этого равенства в том, что среднее арифметическое смещений
∆x очень большого числа частиц приближается к нулю по мере роста этого числа. Так по- нимается усреднение и далее. Иногда такие средние величины называют априор-
ными
19
. Кроме того, мы будем использовать «наблюдаемые средние» – средние арифметические для заданного числа частиц (как правило, очень большого). «На- блюдаемое среднее» смещения частицы
∆xн мало, но не равно нулю.
После каждого шага частицы будут «расползаться» от оси y. Обозначим x
(k)
координату некоторой частицы через k шагов. Тогда
x
(k + 1) = x(k) + ∆x.
(2)
Усреднив это равенство (вновь по множеству частиц), получаем
x(k + 1) = x(k),
т.е. среднее значение
x(k) не изменяется от шага к шагу и, следовательно, равно
x(0) = 0. Наблюдаемое значение xн для большого числа частиц
x(k) н =
1
N
N

j=1
x
j
(k)
(3)
19
Например, приступая к какой-то жеребьевке с помощью подбрасывания монеты, мы за-
ранее предполагаем, что вероятность выпадения «орла» равна 1/2.
41

окажется близким к нулю (здесь x
j
-координата j-й частицы)
20
Ширину полосы, по которой распределяются частицы после k-го шага, удобно характеризовать величиной
x
2
(k). Чтобы выяснить зависимость этой величины от числа шагов, возведем равенство (
2
) в квадрат и усредним:
x
2
(k + 1) = x
2
(k) + 2x(k)∆x + (∆x)
2
.
(4)
Ввиду независимости x
(k) и ∆x имеем
x(k)∆x = x(k)∆x = 0.
Обозначим
(∆x)
2
 = a
2
. Из (4) следует
x
2
(k + 1) = x
2
(k) + a
2
,
т.е. средний квадрат координаты растет с каждым шагом на величину a
2
. Значит,
x
2
(k) = ka
2
.
(5)
Наблюдаемое значение
x
2

н
=
1
N
N

j=1
x
2
j
(6)
изменяется приблизительно пропорционально числу шагов.
Распределение частиц в занятой ими полосе более детально характеризуется функцией распределения f
(x), определяющей концентрацию частиц; dW = f(x)dx
– вероятность того, что координата j-й частицы после k-го шага окажется в ин- тервале x
≤ x
j
≤ x + dx. Теория случайных блужданий дает для достаточно большого числа шагов k распределение Гаусса
f
(x) =
1
√
2πka
2
exp


−
x
2 2ka
2


.
(7)
Наблюдаемая функция распределения получается путем разбиения оси x на ко- нечные интервалы и подсчета числа частиц в каждом из них. Результат подсчета представляется графически ступенчатой кривой – гистограммой (рис.
7
).
Обратим внимание на одно свойство зависимости (
5
). Если укрупнить шаги по времени в l раз, то средний квадрат смещения за один шаг a
2
следует в соот- ветствии с (
5
) заменить на ˜
a
2
= la
2
, а число шагов k – на ˜
k
= k/l. В итоге
x
2
(k) = la
2
· k/l = ˜
a
2
˜k, т.е. вид зависимости (
5
) не изменяется при укрупнении шагов.
20
При заданном числе частиц N это справедливо для не слишком большого номера шага k.
42

Рис. 7. Распределение частиц при диффузии (гистограмма и теоретическая кривая)
6.2.
Оценка параметров движения броуновской частицы в жидкости
Приведем оценки для реального броуновского движения. Средняя скорость хао- тического движения броуновской частицы v
T
определяется так же, как средняя скорость молекулы, соотношением
mv
2
T
∼ k
Б
T,
(8)
где m – масса частицы, k
Б
– постоянная Больцмана, T – абсолютная температура среды. Если скорость частицы v
 v
T
, то ее движение определяется уравнением
m
˙v = F,
где
F = −αv – сила трения
21
. Частица тормозится, и время τ , за которое ее скорость существенно уменьшится, можно оценить, подставив
˙v ∼ v/τ:
mv/τ
∼ αv,
откуда
τ
∼ m/α.
(9)
Если же скорость частицы близка к тепловой, v
∼ v
T
, то и сила, естественно,
гораздо меньше, а отклонения ее от среднего значения
−αv весьма существенны.
21
Для шарика радиуса R в жидкости с коэффициентом вязкости η согласно закону Стокса
α
= 6πηR.
43

Именно эти отклонения ответственны за безостановочное хаотическое движение частицы. Если речь идет о таком движении, то τ из (
9
) можно понимать как оцен- ку времени, спустя которое частица «забывает» начальное направление движения.
Но та же величина τ дает грубую оценку интервала времени, в течение которо- го частица «помнит» направление движения. (Может быть, для оценки времени
«гарантированного забывания» стоило бы взять
2τ, а для оценки времени гаранти- рованного сохранения направления τ /
2, но нас интересуют не «гарантированные»,
а средние времена, поэтому будем полагать, что коэффициенты 2, 1/2 и т.п. лежат за пределами принятой точности оценок.)
За время τ частица проходит путь, равный по порядку величины
a
∼ v
T
τ.
(10)
Смещения частицы за различные интервалы времени порядка τ мы можем рас- сматривать как случайные, подобные рассматривавшимся ранее
∆x, только на- правленные не вдоль оси x, а в произвольном направлении (например, как три одновременных и независимых смещения вдоль трех осей координат). Движение частицы за время t
 τ можно разбить на k ∼ t/τ таких шагов. Смещение ча- стицы за время t оценивается по аналогии с (
5
):
r
2
(t) ∼ ka
2
∼ (v
T
τ
)
2
t
τ
∼ v
2
T
τ t.
(11)
Этот результат обычно представляют в виде
r
2
(t) = 6D t,
(12)
где D – коэффициент диффузии
22
. С учетом (
8
), (
9
), (
11
)
D
∼
k
Б
T
α
.
(13)
Если сначала частицы были сосредоточены в каком-то малом объеме, то со вре- менем они распространяются все дальше, занимая область размера
∼ r(t).
Соотношения вида (
12
), (
13
), полученные Эйнштейном и Смолуховским, по- служили основой экспериментов Перрена, в ходе которых была определена масса атомов и которые были приняты «научной общественностью» в качестве убеди- тельного доказательства существования атомов.
Описанные выше закономерности следует понимать как предельный случай, от- вечающий наблюдению бесконечного числа частиц. Реализация же случайных блу- жданий конечного числа частиц, совершающих броуновское движение (настоящих или «компьютерных»), демонстрирует лишь приближенное выполнение этих соот- ношений.
22