КОМПЬЮТЕРНОЕ МОДЕЛИРОВАНИЕ ФИЗИЧЕСКИХ ЯВЛЕНИЙ

Компьютерное моделирование является одним из эффективных методов изучения физических систем. Часто компьютерные модели проще и удобнее исследовать, они позволяют проводить вычислительные эксперименты, реальная постановка которых затруднена или может дать непредсказуемый результат. Логичность и формализованность компьютерных моделей позволяет выявить основные факторы, определяющие свойства изучаемых объектов, исследовать отклик физической системы на изменения ее параметров и начальных условий.

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

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

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

Принципы моделирования состоят в следующем [3]:

1. Принцип информационной достаточности. При полном отсутствии информации об объекте построить модель невозможно. При наличии полной информации моделирование лишено смысла. Существует уровень информационной достаточности, при достижении которого может быть построена модель системы.

2. Принцип осуществимости. Создаваемая модель должна обеспечивать достижение поставленной цели исследования за конечное время.

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

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

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

Компьютерное моделирование систем часто требует решения дифференциальных уравнений [1-10]. Важным методом является метод сеток, включающий в себя метод конечных разностей Эйлера. Он состоит в том, что область непрерывного изменения одного или нескольких аргументов заменяют конечным множеством узлов, образующих одномерную или многомерную сетку, и работают с функцией дискретного аргумента, что позволяет приближенно вычислить производные и интегралы. При этом бесконечно малые приращения функции f = f(x, y, z, t) и приращения ее аргументов заменяются малыми, но конечными разностями.

1. МЕТОДЫ ЧИСЛЕННОГО ИНТЕГРИРОВАНИЯ И ДИФФЕРЕНЦИРОВАНИЯ

1. Задача. Для известной функции y = y(x), определите первую и вторую производные в точке с координатой x, а также интеграл в интервале от a до b.

2. Теория. Пусть задана функция y = y(x). Разобъем интервал от a до b на элементарные отрезки длиной h = Δx, получив конечное множество узлов сетки x_i = x₁ + ih, где i = 1, 2,..., N, а N -- число узлов. В результате функция непрерывного аргумента будет заменена функцией дискретного аргумента y_i = y(x_i ). Тогда левая, правая и центральная разностные производные первого порядка в точке с координатой x_i соответственно равны:

y'(x_i ) _- = (y(x_i ) - y(x_{i - 1} ))/h,

y'(x_i ) ₊ = (y(x_{i + 1} ) - y(x_i ))/h,

y'(x_i ) = (y(x_{i + 1} ) - y(x_{i - 1} ))/(2h).

Чем меньше шаг сетки h, тем выше точность найденных производных.

Вторая производная определяется из выражения:

y''(x_i ) = (y(x_{i + 1} ) - 2y(x_i ) + y(x_{i - 1} ))/h².

Интеграл функции численно равен площади криволинейной трапеции, ограниченной графиком этой функции и пределами интегрирования a, b. Если эту трапецию разбить на N прямоугольных полосок шириной h = (b - a) / N и длиной y_i = y(a + ih), то площадь будет примерно равна:

Чем меньше шаг h и, соответственно, больше N, тем точнее найденное значение интеграла. Этот метод называется методом прямоугольников.

Более точный метод трапеций заключается в том, что каждая i - ая полоска рассматривается как трапеция высотой h с длинами оснований y_i = y(a + ih) и y_{i + 1} = y(a + (i + 1)h), поэтому ее площадь равна ΔS = (y_i + y_{i + 1} )h/2. Интеграл функции равен сумме всех элементарных площадей этих трапецевидных полосок:

Метод Монте - Карло нахождения площади криволинейной трапеции под кривой y = y(x) состоит в следующем. Представим себе прямоугольник, ограниченный пределами интегрирования a и b, осью x и горизонталью y = c, внутри которого находится эта криволинейная трапеция. Площадь прямоугольника равна (b - a)c. Задавая случайным образом координаты x_i , y_i , поместим внутрь прямоугольника N точек. Подсчитаем число n точек, оказавшихся внутри криволинейной трапеции, то есть удовлетворяющих условию y_i <y(x_i ). Площадь криволинейной трапеции будет во столько раз меньше площади выбранного прямоугольника, во сколько раз n меньше N. Поэтому при N стремящемся к бесконечности дробь n(b - a)c / N стремится к пределу, равному искомому интегралу.

3. Компьютерная программа. Самостоятельно составьте алгоритмы нахождения производных и интегралов. Ниже представлены примеры программ. Первая программа позволяет вычислить первую и вторую производные функции y = x³ - x² + 3 в точке с координатой x = 0, а также найти ее интеграл в интервале от a = 1 до b = 3 методом трапеций. Вторая программа определяет интеграл функции y = x² в интервале от 0 до 1 методом Монте - Карло.

program PROGRAMMA1_1;
uses crt;
var  x,y1,y2,y3,a,b,h,S :real;
Function Funct(x:real):real;
begin  {Задание функции}
   Funct:=x*x*x-x*x+3; 
end;
BEGIN  {Основная программа}
   clrscr; x:=3; h:=0.001;
   y1:=Funct(x-h); 
   y2:=Funct(x); 
   y3:=Funct(x+h);
   Writeln('Первая производная  ', (y2-y1)/h:3:3);
   Writeln('Вторая производная  ', (y1-2*y2+y3)/(h*h):3:3);
   a:=1; b:=3; x:=a; S:=0;
   Repeat   {Интеграл}
   S:=S+0.5*(Funct(x)+Funct(x+h))*h; x:=x+h;
   until x>b;
   Writeln('Интеграл  ',S:3:3);
   Repeat until KeyPressed;
END.

program PROGRAMMA1_2;
uses crt;
const NN=10000;
var x,y,xx,yy: real;
    n,i: integer;
function Funct(x:real):real;
begin Funct:=x*x; end;
BEGIN  {Основная программа}
  clrscr; Randomize; n:=0;
  for i:=1 to NN do
    begin
       x:=Random(1000)/1000; 
       yy:=Random(1000)/1000;
       if yy<Funct(x) then n:=n+1;
    end;
  writeln('Интеграл равен ',n/NN);
  Repeat until KeyPressed;
END.

4. Задания для самостоятельного решения.

1. Точка движется прямолинейно со скоростью v(t) = 3 + 2t + 1t². Численными методами и аналитически определите ускорение точки в момент времени 1 с и пройденный путь за время от 1 до 3 с. Повторите расчеты при различных значениях шага и сравните результаты.

2. Имеется пластинка толщиной h ограниченная кривой y = x² и прямой y = 1. Ее плотность есть функция координаты y: ρ(y) = ρ₀ (1 + α y), где α -- произвольный коэффициент пропорциональности. Определите ее площадь и массу методом Монте - Карло.

3. Определите координаты центра масс пластины толщины h, ограниченной прямыми x = 0, y = 0, y = 4 - x², плотность которой равна ρ.

4. Пластина толщиной h имеет форму круга радиуса R. Ее плотность с ростом расстояния r до центра убывает по закону ρ(r) = ρ₀ (1,5 - r/R). Методом численного интегрирования определите момент инерции пластины относительно оси проходящей через ее центр и лежащей в ее плоскости.

8. Постройте кривую зависимости излучательной способности абсолютно черного тела от частоты при постоянной температуре T, выражаемую формулой Планка: f_ω (ω,T) = Aω³/ (e^Bω/T - 1), где A и B -- постоянные коэффициенты. Постройте график при различных T. Методом численного интегрирования найдите интегральную светимость абсолютно черного тела, взяв интеграл от f_ω (ω,T).

---

1. Задача. Имеется физическая система с одной степенью свободы, состоящая из инерционного элемента массой m, упругого элемента жесткостью k и диссипативного элемента с коэффициентом сопротивления r. Определить отклик системы x(t), а также ее первую и вторую производные v(t), a(t), на внешнее воздействие F_x = F_x (t), если известны начальные условия x(0), v(0).

2. Теория. Из второго закона Ньютона следует линейное неоднородное дифференциальное уравнение второго порядка:

Процессы, происходящие в последовательном колебательном контуре, состоящем из последовательно соединенных резистора R, конденсатора C и катушки индуктивности L, на который подано напряжение U(t), описываются уравнением:

где q -- заряд, проходящий по цепи, i = dq/dt -- сила тока. При этом реализуется электромеханическая аналогия "сила - напряжение".

Итак, неоднородное диффуравнение второго порядка описывает широкий класс задач, изучаемых в школьном и вузовском курсе физики: движение связанного с пружиной тела в вязкой среде под действием произвольной силы F_x (t), протекание тока через последовательно соединенные резистор, конденсатор и катушку индуктивности, подключенные к источнику ЭДС e(t) или источнику тока i(t).

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

• внешняя сила отсутствует;
• внешняя сила постоянна;
• внешняя сила изменяется по гармоническому закону;
• внешняя сила изменяется по произвольному периодическому закону;
• внешняя сила изменяется по произвольному непериодическому закону.

Кроме того, физические явления, возникающие в системе, зависят от ее параметров m, r, k и начальных условий, к которым относятся координата x(0) и скорость v_x (0) в начальный момент времени.

3. Алгоритм. Дифференциальное уравнение второго порядка может быть решено методом конечных разностей Эйлера. Он состоит в том, что бесконечно малые приращения функции x(t) и ее первых двух производных v_x , a_x заменяются конечными разностями, а бесконечно малые приращения dt -- малыми, но конечными приращениями Δ t. Сначала, исходя из параметров системы m, r, k, ее координаты x и скорости v_x в момент времени t, расчитывается ее координата и скорость в следующий момент t + Δ t:

a_x (t + Δ t) = (F_x (t) - r v_x (t) - kx(t))/m,

v_x (t + Δ t) = v_x (t) + a_x (t + Δ t)Δ t.

x(t + Δ t) = x(t) + v_x (t + Δ t)Δ t.

Затем это состояние рассматривается как исходное и процедура расчета повторяется для момента времени t + 2Δ t и так далее. Одновременно с вычислениями строятся графики x = x(t), v_x = v_x (t), a_x = a_x (t).

Рассмотрим алгоритм численного решения уравнения.

1. Задают параметры физической системы m, k, r, зависимость внешнего воздействия от времени F_x (t), а также начальные условия x(0), v_x (0) и шаг по времени Δ t.

2. Начало цикла по t. Дают приращение по времени: переменной t присваивают значение t + Δ t.

3. Определяют ускорение, скорость и координату тела в момент t + Δ t:

a_x (t + Δ t) = (F_x (t) - r v_x (t) - kx(t))/m,

v_x (t + Δ t) = v_x (t) + a_x (t + Δ t)Δ t,

x(t + Δ t) = x(t) + v_x (t + Δ t)Δ t.

4. Результаты вычислений x(t + Δ t), v_x (t + Δ t), a_x (t + Δ t) выводят на экран в числовом виде либо строят соответствующие точки на координатной плоскости.

5. Возвращение к операции 2. Если цикл по t закончился, -- выход из цикла.

4. Компьютерная программа. При запуске программа рисует графики зависимостей координаты x = x(t), проекции скорости v_x = v_x (t) и проекции ускорения a_x = a_x (t).

Некоторые строчки программы заключены в скобки "(*" и "*)". Убрав скобки и активизировав соответствующие операторы, можно промоделировать различные явления.

program PROGRAMMA2;
uses  dos, crt, graph;
Const Fm=10;w=5;m=2;r=0;k=0;
Mx=20; Mv=40; Ma=8; Mf=2; Mt=100;
dt=0.00006;
Var x,v,a,F,t : Real;
j,xx,vv,aa,FF,tt,Gd,Gm : Integer;
BEGIN
  Gd:= Detect;
  InitGraph(Gd, Gm, 'c:\bp\bgi');
  if GraphResult<>grOk then Halt(1);
  t:=0; v:=0; x:=-3;
  line(30,300,650,300);
  line(31,500,31,10);
  OutTextXY(50,20,'X, V, A');
Repeat
  begin {Задание функции F=F(t)}
  t:=t+dt; (* F:=Fm*sin(w*t); *)
(*If sin(w*t)<0 then F:=0;
  If sin(w*t)>0 then F:=Fm;*)
  F:=0; If t<1 then F:= Fm;
  If t>3 then F:=-Fm;
  a:=(F-r*v-k*x)/m; x:=x+v*dt;  v:=v+a*dt; tt:=round(t*Mt);
  xx:=round(x*Mx); vv:=round(v*Mv); aa:=round(a*Ma); FF:=round(F*Mf);
  circle(30+tt,300-xx,1); circle(30+tt,300-vv,1); circle(30+tt,300-aa,2);
  end;
until KeyPressed;
CloseGraph;
END.

5. Задания для самостоятельного решения.

1. На точку массы m действует скачкообразно изменяющаяся сила. Исследуйте движение точки, проанализируйте получившиеся графики зависимостей x = x(t), v_x = v_x (t), a_x = a_x (t).

2. Промоделируйте движение материальной точки, движущейся в вязкой среде под действием постоянной силы, направленной вдоль оси x: F_x = const>0 при начальных условиях x₀ = 0, v_x0 <0. Проанализируйте получающиеся графики x(t), v_x (t), a_x (t). Докажите, что время подъема камня, брошенного вертикально вверх, меньше времени спуска.

3. Создайте модель переходного процесса в цепи, содержащей резистор R и катушку индуктивности L, подключенные к источнику постоянного напряжения, при условии, что i₀ не равно 0. Исследуйте аналогичный переходный процесс в цепи, содержащей последовательно соединенные резистор и конденсатор.

4. Изучите движение колебательной системы в случае слабого затухания, когда r/2m<ω₀ = (k/m)^1/2. Убедитесь в том, что ускорение изменяется в противофазе с координатой, а скорость опережает координату на π/2, причем амплитуды колебаний x(t), v(t), a(t) уменьшаются по экспоненте. Проведите серию вычислительных экспериментов при различных начальных условиях системы.

5. Промоделируйте движение осциллятора в случае сильного затухания при r/2m>ω₀ = (k/m)^1/2. Убедитесь, что в этом случае движение будет апериодическим.

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

7. Промоделируйте работу сглаживающего RL - фильтра при подаче на него пульсирующего напряжения, получающиегося в результате однополупериодного выпрямления. Убедитесь в том, что с ростом индуктивности уменьшается коэффициент пульсаций тока и напряжения на резисторе. Изучите зависимость амплитуды пульсаций от индуктивности L, сопротивления нагрузки R и частоты импульсов ω.

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

Так как в программе осуществляется деление на m (аналог индуктивности L ), то значение этого параметра должно быть очень малым, но не равным нулю.

9. Промоделируйте движение тела в вязкой среде (m, r не равны 0, k = 0), на которое в момент времени t₀ = 0 начинает действовать внешняя гармоническая сила F_x = F_m sin(ω t). Эта ситуация соответствует переходному процессу, происходящему при подключении активно - индуктивной нагрузки к источнику переменного напряжения. При t стремящемся к бесконечности переходный ток стремится к принужденному току, изменяющемуся с той же частотой, что и приложенная ЭДС и отстающему от нее на некоторую фазу.

10. Создайте программу, моделирующую процессы, происходящие в колебательной системе в случае, если на нее действует периодически изменяющаяся сила, частота которой пропорциональна времени: F_x (t) = F_msin(ω(1 + α t)t), где α >0. Значения ω и α подберите так, чтобы резонансная частота колебательной системы находилась в середине рабочего диапазона частот. На рисунке показан получающийся график зависимости x = x(t). Так как частота колебаний прямопропорциональна времени, то огибающая графика является амплитудо - частотной характеристикой колебательной системы, и называется резонансной кривой.

1. Задача. Материальная точка массой m движется в силовом поле F_x = F_x(x,y), F_y = F_y(x,y), при этом на нее действует сила вязкого трения с проекциями F_тx = - r v_x, F_тy = - r v_y, направленная противоположно скорости. Необходимо, зная начальные условия x₀ , y₀ , v_0x , v_0y , построить траекторию движения точки.

2. Теория. Примерами подобного движения являются движение точки, в однородном силовом поле, в центральном силовом поле сил притяжения или отталкивания, в центральном поле сил упругости и т.д. При этом могут быть учтены силы вязкого трения.

Проанализируем основные ситуации.

1. Движение в однородном поле. Во всех точках пространства вектор силы имеет постоянные проекции на оси координат. При отсутствии силы трения точка движется по параболе, а при ее наличии -- по более сложной кривой.

2. Движение в центрально - симметричном поле, действующем по закону обратных квадратов. На точку с координатами x, y действует сила F = GmM/r², r² = x² + y² Ее проекции на оси координат:

F_x = - Fcosα = - Fx/r,

F_y = - Fsinα = - Fy/r.

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

3. Движение в магнитном поле. Движение заряженной частицы в магнитном поле будет двумерным, если начальная скорость частицы перпендикулярна силовым линиям магнитного поля. При этом со стороны поля действует сила Лоренца F = qvB, лежащая в плоскости экрана и направленная перпендикулярно вектору скорости. Введем угол β, который образует вектор скорости с осью x. Проекции силы Лоренца на координатные оси:

F_x = - Fsinβ = Fv_y /v,

F_y = - Fcosβ = - Fv_x /v.

Заряженная частица описывает окружность. При наличии тормозящей силы радиус окружности уменьшается.

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

Если заряд частицы положительный, то на него со стороны электрического поля действует постоянная сила, направленная вверх. Чтобы учесть ее влияние необходимо к вертикальной проекции силы Лоренца прибавить постоянное слагаемое qE:

F_x = Fv_y /v, F_y = qE - Fv_x /v.

Если начальная скорость частицы равна нулю, то траекторией ее движения является циклоида.

3. Алгоритм. Пусть в момент времени t материальная точка имеет координаты x, y и проекции скорости v_x , v_y . Запишем второй закон Ньютона в проекциях:

F_x(x,y) -r v_x = ma_x, F_y(x,y) -r v_y = ma_y.

Отсюда следует, что проекции ускорения точки в момент времени t + Δ t равны:

a_x (t + Δ t) = (F_x (t) - r v_x (t))/m, a_y (t + Δ t) = (F_y (t) - r v_y (t))/m.

Определив координаты и проекции скорости точки в момент времени t + Δ t, можно повторить процедуру вычисления требуемое количество раз и построить траекторию движения точки.

Построим алгоритм модели.

1. Задают массу материальной точки m, коэффициент вязкости r, начальные координаты x₀ , y₀ и проекции скорости v_0x , v_0y , силовое поле F_x = F_x (x,y,z), F_y = F_y (x,y,z), а также шаг по времени Δ t.

2. Начало цикла по t. Дают приращение по времени: переменной t присваивают значение t + Δ t.

3. Определяют ускорение, скорость и координату тела в следующий момент времени:

a_x (t + Δ t) = (F_x (t) - r v_x (t))/m,

a_y (t + Δ t) = (F_y (t) - r v_y (t))/m,

v_x (t + Δt) = v_x (t) + a_x (t + Δt)Δ t,

v_y (t + Δ t) = v_y (t) + a_y (t + Δ t)Δ t,

x(t + Δ t) = x(t) + v_x (t + Δ t)Δ t,

y(t + Δ t) = y(t) + v_y (t + Δ t)Δ t.

4. Результаты вычислений x(t + Δ t), y(t + Δ t) выводят на экран в числовом виде либо строят соответствующие точки на координатной плоскости.

5. Возвращение к операции 2. Если цикл по t закончился, -- выход из цикла.

4. Компьютерная программа. Предлагаемая компьютерная программа позволяет изучить движение материальной точки в различных силовых полях с учетом действующей на точку силы трения.

program PROGRAMMA3;
uses crt, graph;
var v, B, q, F, Fx, Fy : real;
  r, x, y, vx,vy,ax,ay : real; Gd, Gm, i: integer;
const M=500; mm=100; dt=0.005; rr=0.1; k=2;
Begin
 Gd:= Detect; InitGraph(Gd, Gm, 'c:\bp\bgi');
 if GraphResult <> grOk then Halt(1);
 line(320,240,640,240); line(320,240,320,0); circle(320,240,5);
 x:=100; y:=120; vx:=1; vy:=-2;
Repeat
  begin
{--Заданние силового поля--}
(*   Fy:=3;        Fx:=0;           *)
(*   Fx:=-k*x;     Fy:=-k*y;        *)
(*   r:=sqrt(x*x+y*y); F:=M*mm/(r*r);
     Fx:=-F*x/r;       Fy:=-F*y/r; *)
     B:=2; q:=1; F:=B*v*q; v:=sqrt(vx*vx+vy*vy);
     Fx:=F*vy/v; Fy:=-F*vx/v;
(*   B:=2; q:=1; F:=B*v*q; v:=sqrt(vx*vx+vy*vy);
     Fx:=F*vy/v; Fy:=-0.5-F*vx/v;   *)
{--Расчет скоростей и ускорений--}
ax:=(Fx-rr*vx)/mm; ay:=(Fy-rr*vy)/mm;
vx:=vx+ax*dt; vy:=vy+ay*dt; x:=x+vx*dt; y:=y+vy*dt;
circle(round(x)+320,240-round(y),2);   setcolor(12);
circle(round(x)+320,240-round(y),1);   setcolor(15);
end;
until KeyPressed;
CloseGraph;
END.

5. Задания для самостоятельного решения.

1. Материальная точка массы m брошена под углом к горизонту в однородном поле тяжести. Изучите траекторию ее движения при отсутствии силы вязкого трения и при ее наличии.

2. Убедитесь в том, что время подъема материальной точки, движущейся под действием силы тяжести в вязкой среде, меньше времени спуска. Воспользуйтесь тем, что в наивысшей точке подъема проекция скорости на ось y меняет свой знак на противоположный.

3. Промоделируйте движение точки в поле центральных сил упругости F_x = - kx, F_y = - ky, в случае, когда на точку действует сила вязкого трения и когда она равна нулю. По какой траектории движется точка?

4. Исследуйте движение точки в поле сил притяжения, действующих по закону обратных квадратов F = GmM / r². Промоделируйте ситуации, в которых точка движется по гиперболе, параболе, эллипсу. Изучите характер движения искусственного спутника Земли, входящего в верхние слои атмосферы, на который действует сила вязкого трения.

5. Изучите движение точки в поле сил отталкивания. Промоделируйте опыт Резерфорда по отклонению альфа - частиц ядрами атомов золота. Меняя прицельный параметр, проведите серию компьютерных экспериментов.

6. Промоделируйте движение заряженной частицы в камере Вильсона, помещенной в однородное магнитное поле. Учтите, что по мере своего движения частица теряет кинетическую энергию. Магнитное поле направлено перпендикулярно плоскости экрана.

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

8. Изучите движение материальной точки в гравитационном поле двух массивных тел. Проведите компьютерные эксперименты при заданных начальных координатах и скорости точки.

9. Задайте произвольное силовое поле F_x = F_x (x,y), F_y = F_y (x,y) и промоделируйте движение частицы в нем.

1. Задача. Имеется система N материальных точек с массами m_i , i = 1, 2, ..., N, взаимодействующих друг с другом с внутренними силами, на каждую из которых действует внешняя сила. Исходя из начальных координат x_i, y_i, и скоростей v_xi, v_yi, определите координаты и скорости материальных точек в последующие моменты времени.

2. Теория. Рассмотрим механическую систему из N материальных точек. Основной закон динамики:

где F'_ij -- внутренняя сила, действующая на i-ую материальную точку со стороны j -ой материальной точки, F_i -- равнодействующая внешних сил, действующих на i - ую материальную точку со стороны тел не входящих в систему.

Дифференциальное уравнение второго порядка может быть представлено двумя дифференциальными уравнениями первого порядка. Имеем:

Зная внешние и внутренние силы, действующие на каждую материальную точку, можно определить их ускорения. Исходя из координат и скоростей точки в момент времени t, можно расчитать координаты и скорости точки в следующий момент времени t + Δ t.

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

3. Алгоритм.

1. Задают число материальных точек N, их массы m_i , координаты x_i , y_i и проекции начальных скоростей v_ix , v_iy , силовое поле F_x = F_x (x,y), F_y = F_y (x,y), а также шаг по времени Δ t.

2. Начало цикла по t. Дают приращение по времени: переменной t присваивают значение t + Δ t.

3. Определяют проекции F_xi , F_yi равнодействующей всех внешних и внутренних сил, действующих на каждую i - ую материальную точку в момент t + Δt, и записывают их в массивы.

4. В цикле переобозначают координаты всех материальных точек, записывая их в массивы xx[i], yy[i].

5. В цикле перебирают все материальные точки и определяют проекции ускорения, скорости и координаты для каждой из них в момент t + Δt:

a_xi (t + Δt) = F_xi (t + Δt)/m_i ,

v_ix (t + Δ t) = v_ix (t) + a_ix (t + Δt)Δt,

x_i (t + Δ t) = x_i (t) + v_ix (t + Δt)Δt.

По аналогичным формулам вычисляют проекции на ось OY. Результаты записывают в массивы x[i], y[i], vx[i], vy[i].

6. Стирают изображения материальных точек в предыдущий момент времени t, координаты которых сохранены в массивах xx[i], yy[i].

7. На экране строят точки в следующий момент t + Δt, либо рисуют графики или выводят результат в числовом виде.

8. Возвращение к операции 2. Если цикл по t закончился, -- выход из цикла.

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

program PROGRAMMA4;
uses dos, crt, graph;
const  N=50; dt=0.01;
var m,Fx,Fy,x,y,vx,vy,xx,yy : array[1..N] of real;
    Gd, Gm, i, j : integer; ax, ay, F, l : real;
label     Metka,  metka1;
Procedure Init_Graph; {Инициализация графики}
begin
  Gd:=Detect; InitGraph(Gd, Gm, 'c:\bp\bgi');
  if GraphResult <> grOk then Halt(1);
end;
Procedure Sila; {Вычисление действующих сил}
label     Metka;
begin
  For i:=1 to N do begin Fx[i]:=0; Fy[i]:=0; end;
  For i:=1 to N do for j:=1 to N do begin
    if j=i then goto Metka;
    l:=sqrt(sqr(x[i]-x[j])+sqr(y[i]-y[j])); if l<2 then l:=2;
    F:=-50000*m[i]*m[j]/sqr(l)+500000*m[i]*m[j]/sqr(l*l);
    Fx[i]:=Fx[i]+F*(x[i]-x[j])/l; 
    Fy[i]:=Fy[i]+F*(y[i]-y[j])/l+m[i]*10;
    Metka:  end;
end;
Procedure Nach_uslov;
begin  Randomize; {Задание случайных координат и скоростей}
  for i:=1 to N do begin m[i]:=2;
  x[i]:=random(280)+60; y[i]:=random(280)+60;
  vy[i]:=random(30)-15; vx[i]:=random(30)-15;
end; end;
BEGIN
Init_Graph; Nach_uslov;
Repeat Sila;
for i:=1 to N do
  begin
    xx[i]:=x[i];           yy[i]:=y[i]; {Запись предыдущих координат}
    ax:=Fx[i]/m[i];        ay:=Fy[i]/m[i]; {Вычисление ускорений,}
    vx[i]:=vx[i]+ax*dt;    vy[i]:=vy[i]+ay*dt;  {скоростей,}
    x[i]:=x[i]+vx[i]*dt;   y[i]:=y[i]+vy[i]*dt; {координат}
    if (x[i]<50)or(x[i]>350) then vx[i]:=-vx[i];{отражение}
    if (y[i]<50)or(y[i]>350) then vy[i]:=-vy[i];{от стенок}
  end;   delay(500); 
setcolor(8);  for i:=1 to N do circle(round(xx[i]),round(yy[i]),2);
setcolor(15); for i:=1 to N do circle(round(x[i]),round(y[i]),2);
until KeyPressed; Repeat until keypressed; CloseGraph;
END.

5. Задания для самостоятельного решения.

1. На основе компьютерной модели изучите движение двух или трех частиц, между которыми действуют силы притяжения.

2. Изучите движение двух материальных точек, между которыми действуют силы отталкивания. Промоделируйте центральный и нецентральный удар.

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

4. Изучите движение материальной точки в гравитационном поле двух массивных тел. Проведите вычислительные эксперименты при различных начальных координатах и скоростях точки.

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

6. Промоделируйте движение молекул газа в прямоугольном замкнутом сосуде. Учтите, что при соударении молекулы с горизонтальной (вертикальной) стенкой сосуда вертикальная (горизонтальная) проекция ее скорости меняет свой знак на противоположный.

7. Промоделируйте диффузию двух газов. Пусть вначале молекулы с массами m заполняют левую половину сосуда, а молекулы с массами M -- правую. Задайте случайные значения скоростей молекул. Как изменяется концентрация молекул газов в сосуде с течением времени?

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

9. Промоделируйте движение молекул газа в гравитационном поле шара большой массы.

1. Задача: Имеется система из N осцилляторов с массами m_i и жесткостью пружин k_i , связанные между собой упругими связями с жесткостью q_i . На каждый осциллятор действует сила вязкого трения - r η_i . Заданы начальные смещения ξ_i (0) и скорости η_i (0) всех осцилляторов. На отдельные осцилляторы действует вынуждающая сила F_iξ = F_iξ (t) некоторые осцилляторы колеблются по закону ξ_i = ξ_i (t). Известно, что крайние осцилляторы закреплены (свободны). Необходимо построить компьютерную модель колебаний связанных осцилляторов.

2. Теория. Рассмотрим цепочку осцилляторов с массами m_i и коэффициентом жесткости k_i , связанных между собой пружинами жесткостью q_i . На каждый i - ый осциллятор со стороны соседних (i - 1) - ого и (i + 1) - ого осцилляторов действует сила:

F_iξ = - q_i (ξ_i - ξ_{i - 1} ) - q_{i + 1} (ξ_i - ξ_{i + 1} ),

Кроме нее на осциллятор действует сила упругости - kξ_i и сила вязкого трения - r η_i , где η_i -- скорость i - ого осциллятора. По второму закону Ньютона, ускорение i - ого осциллятора равно:

θ_i = (F_iξ - r η_i - kξ_i )/m.

Определив ускорение θ_i , можно найти скорость и смещение i - ой массы в следующий момент времени t + Δt:

η_i = η_i + θ_i Δ t,

ξ_i = ξ_i + η_i Δ t.

3. Алгоритм. Эта задача является частным случаем задачи о моделировании движения системы взаимодействующих между собой частиц.

1. Задают параметры системы и начальные условия: число осцилляторов N, их массы m_i , жесткость упругих связей k_i , жесткость пружины внутри осциллятора q_i , шаг по времени Δt, начальные смещение ξ_i и скорости η_i , вынуждающие силы, действующие на различные материальные точки.

2. Дают приращение по времени: переменной t присваивают значение t + Δt.

3. В цикле перебирают все N осцилляторов и для каждого из них опеределяют равнодействующую силу, ускорение, скорость и смещение в следующий момент t + Δ t:

F_iξ (t + Δt) = - q(ξ_i (t) - ξ_{i - 1} (t)) - q(ξ_i (t) - ξ_{i + 1} (t)),

θ_i (t + Δt) = (F_i (t + Δ t) - r η_i (t) - kξ_i (t))/m,

η_i (t + Δt) = η_i (t) + θ_i (t + Δt)Δt,

ξ_i (t + Δt) = ξ_i (t) + η_i (t + Δt)Δt.

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

5. Возвращение к операции 2. Если цикл по t закончился, -- выход из цикла.

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

program PROGRAMMA5;
uses  dos, crt, graph;
Const   m=0.5; r=0.1; k=0.01; 
        dt=0.001; q=100; N=50;
Var     teta, F, t, y  : Real;
i,xx,vv,aa,FF,tt,Gd,Gm  : Integer;
eta, ksi : array [0..N] of real;
Procedure GraphInit;
begin
  Gd:= Detect;
  InitGraph(Gd, Gm, 'c:\bp\bgi');
  if GraphResult <> grOk then Halt(1);
end;
Procedure Oscillator;
  begin
  F:=q*(ksi[i-1]-ksi[i])+q*(ksi[i+1]-ksi[i]);
  teta:=(F-r*eta[i]-k*ksi[i])/m;
  eta[i]:=eta[i]+teta*dt; ksi[i]:=ksi[i]+eta[i]*dt;
end;
BEGIN
  GraphInit;
  Repeat  begin  t:=t+dt;
  For i:=1 to N do  begin
    y:=ksi[i]; Oscillator;
    if 5*t<3.142 then ksi[1]:=20*sin(5*t) else ksi[1]:=0;
    ksi[N]:=0; {закреплен}
    {ksi[N]:=ksi[N-1]; незакреплен}
    setcolor(8); circle(10*i,240-round(y*10),2);
    setcolor(15); circle(10*i,240-round(ksi[i]*10),2);
  end; end;
  until KeyPressed;
CloseGraph;
END.

5. Задания для самостоятельного решения.

1. Промоделируйте колебания двух связанных осцилляторов. Рассмотрите случаи: 1) на один из них действует вынуждающая сила; 2) один из осцилляторов имеет начальное смещение; 3) один из осцилляторов имеет начальную скорость. Выполните компьютерные эксперименты при различных k_i и q_i .

2. Изучите колебания трех связанных осцилляторов, рассмотрев все перечисленные выше случаи, выполнив компьютерные эксперименты при различных k_i и q_i .

3. Промоделируйте колебания 50 осцилляторов, связанных упругими связями, в случае, когда на левый крайний осциллятор подействовала кратковременная сила. Рассмотрите случаи, когда правый крайний осциллятор закреплен и незакреплен.

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

7. Изучите распространение импульса и его отражение от открытого или закрытого конца струны (одномерной упругой среды), которая моделируется 50 связанными осцилляторами.

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

2. Теория. Рассмотрим обобщенную модель Винера - Розенблюта. Мысленно разобъем экран компьютера на элементы, определяемые индексами i, j и образующими двумерную сеть. Пусть состояние каждого элемента описывается фазой y_i,j (t), и концентрацией активатора u_ij (t), где t -- дискретный момент времени.

Если элемент находится в покое, то будем считать, что y_i,j (t) = 0. Если вследствие близости возбужденных элементов концентрация активатора u_ij (t) достигает порогового значения h, то элемент возбуждается и переходит в состояние 1. Затем на следующем шаге он переключается в состояние 2, затем -- в состояние 3 и т.д., оставаясь при этом возбужденным. Достигнув состояния r элемент переходит в состояние рефрактерности. Через (s - r) шагов после возбуждения элемент возвращается в состояние покоя.

Будем считать, что при переходе из состояния s в состояние покоя 0 концентрация активатора становится равной 0. При наличии соседнего элемента, находящегося в возбужденном состоянии, она увеличивается на 1. Если p ближайших соседей возбуждены, то на соответствующем шаге к предыдущему значению концентрации активатора прибавляется число возбужденных соседей: u_ij (t + Δt) = u_ij (t) + p.

Можно ограничиться учетом ближайших восьми соседних элементов.

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

1. Задают число элементов активной среды, ее параметры s, r, h, начальное распределение возбужденных элементов.

2. Начало цикла по t. Дают приращение по времени: переменной t присваивают значение t + Δt.

3. Перебирают все элементы активной среды, определяя их фазы y_i,j (t + Δt) и концентрацию активатора u_i,j (t + Δt) в момент t + Δt.

4. Очищают экран и строят возбужденные элементы активной среды.

7. Возвращение к операции 2. Если цикл по t закончился, -- выход из цикла.

4. Компьютерная программа. Ниже представлена программа, моделирующая активную среду и происходяшие в ней процессы. В программе заданы начальные значения фазы y_i,j (t + Δt) всех элементов активной среды, а также имеется цикл по времени, в котором расчитываются значения y_i,j (t + Δt) в следующий момент t + Δt и осуществляется графический вывод результатов на экран. Параметры среды r = 6, s = 13, h = 5, то есть каждый элемент кроме состояния покоя может находиться в 6 возбужденных состояниях и 7 состояниях рефрактерности. Пороговое значение концентрации активатора равно 5. Программа строит однорукавную волну, осциллятор и препятствие.

Program PROGRAMMA6;
uses dos, crt, graph; 
Const N=110; M=90; s=13; r=6; h=5; 
Var      y, yy, u     : array [1..N,1..M] of integer;
         ii, jj, j, k, Gd, Gm : integer; i : Longint;
Label met;

BEGIN  
  Gd:= Detect;  InitGraph(Gd, Gm, 'c:\bp\bgi');
  If GraphResult <> grOk then  Halt(1);
  setcolor(8); setbkcolor(15);
(*  y[50,50]:=1;              { Одиночная волна } *)
    For j:=1 to 45 do         { Однорукавная волна }
    For i:=1 to 13 do  y[40+i,j]:=i;
(*  For j:=1 to M do          { Двурукавная волна  }
    For i:=1 to 13 do begin    y[40+i,j]:=i;
    If j>40 then y[40+i,j]:=14-i; end;   *)
Repeat 
    If k=round(k/20)*20 then y[30,30]:=1; {Осциллятор 1}
(*  If k=round(k/30)*30 then y[20,50]:=1; {Осциллятор 2} *)
    For i:=2 to N-1 do For j:=2 to M-1 do begin
    If (y[i,j]>0) and (y[i,j]<s) then yy[i,j]:=y[i,j]+1;
    If y[i,j]=s then begin yy[i,j]:=0; u[i,j]:=0; end;
    If y[i,j] <> 0 then goto met;
    For ii:=i-1 to i+1 do For jj:=j-1 to j+1 do begin
    If (y[ii,jj]>0) and (y[ii,jj]<=r) then u[i,j]:=u[i,j]+1;
    If u[i,j]>=h then yy[i,j]:=1;                end;
met:  end;  Delay(2000);               {Задержка}
    cleardevice; 
    For i:=21 to 70 do begin
    yy[i,60]:=0; yy[i,61]:=0;        {Препятствие}
    circle(6*i-10,500-6*60,3); circle(6*i-10,500-6*61,3); end;
    For i:=1 to N do For j:=1 to M do 
      begin y[i,j]:=yy[i,j]; setcolor(12);
      If (y[i,j]>=1) and (y[i,j]<=r) then circle(6*i-10,500-6*j,3);
      setcolor(8);
      If (y[i,j]>6) and (y[i,j]<=s) then circle(6*i-10,500-6*j,2);
      end;
until KeyPressed; 
CloseGraph;
END.

5. Задания для самостоятельного решения.

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

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

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

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

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

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

7. Промоделируйте образование двухрукавных спиральных волн.

8. Изучите зависимость частоты вращения однорукавной спиральной волны от параметров среды (r, s, h). Повторите этот вычислительный эксперимент для двухрукавной волны.

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

10. Исследуйте распространение и анигиляцию одиночного импульса в одномерной активной среде.

11. Изучите распространение автоволн в одномерной активной среде при наличии осциллятора.

12. Промоделируйте распространение одиночного импульса в одномерной активной среде, последний элемент которой контактирует с первым.

13. Создайте компьютерную модель распространения автоволн в трехмерной активной среде.

ЛИТЕРАТУРА

1. Бурсиан Э.В. Физика 100 задач для решения на компьютере: Учебное пособие. -- СПб.: ИД "МиМ", 1997. -- 256 с.
2. Гулд Х., Тобочник Я. Компьютерное моделирование в физике: В 2-х частях. Часть первая.-- М.: Мир, 1990.-- 400 с.
3. Гультяев А. Визуальное моделирование в среде MATLAB: учебный курс -- СПб.: Питер, 2000. -- 432 c.
4. Заславский Г.М., Сагдеев Р.З. Введение в нелинейную физику: От маятника до турбулентности и хаоса.-- М.: Наука, 1988.-- 368 с.
5. Зельдович Я.Б. Высшая математика для начинающих и ее приложения к физике. -- М.: Наука, 1963. -- 560 с.
6. Зельдович Я.Б., Мышкис А.Д. Элементы прикладной математики. -- М.: Наука, 1965. -- 615 с.
7. Лоскутов А.Ю., Михайлов А.С. Введение в синергетику: Учеб. руководство.-- М.: Наука, 1990.-- 272 с.
8. Санин А.Л. Структуры и хаос. -- проблемы физики.-- Л.: Знание, 1985.-- 32 с.
9. Могилев А.В., Пак Н.И., Хеннер Е.К. Информатика: Учеб. послобие для студ. пед. вузов / Под ред. Е.К.Хеннер.--- М.: Изд.центр "Академия", 2001. --- 816 с.
10. Тихонов А.Н., Самарский А.А. Уравнения математической физики. -- М.: Наука, 1966. -- 724 с.