36424

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

Лабораторная работа

Информатика, кибернетика и программирование

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

Русский

2013-09-21

161.5 KB

35 чел.

Лабораторная работа 7

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

Материалы для заданий 1, 2

Потери пучка при прохождении через вещество

В этой работе можно познакомиться с основным методом моделирования, применяемым при исследовании прохождения пучков частиц через вещество — методом статистического моделирования (называемым методом Монте-Карло). При этом «судьба» каждой частицы «разыгрывается» с помощью случайного выбора, а полученные для множества частиц результаты подвергаются статистической обработке. Метод применяется, например, при проектировании ядерных реакторов, детекторов частиц на ускорителях и обработке получаемых результатов (а также во многих других случаях, скажем, при исследовании распространения мутаций в среде живых организмов).

Мы будем изучать, естественно, очень простой вариант задачи — прохождение пучка тяжелых частиц (A) через слой газа, состоящего из легких частиц (O). К тому же будем полагать скорости частиц газа до столкновения с частицами пучка равными нулю. Не будем также учитывать изменений, происходящих в газе вследствие прохождения пучка. Это не ограничивает существенно общности задачи.

1. Эффективные сечения

Частицы можно представлять шариками радиусов Ra и Ro. Пусть частица A летит так, что центр ее должен пролететь на расстоянии ρ от центра частицы O. Столкновение произойдет, если ρ < R, где R = Ra + RO. Фактически частица О образует для частиц A «преграду», площадь которой σ = πR2 . Эта величина называется эффективным сечением столкновения . В зависимости от величины р (называемой прицельным параметром) определяются угол отклонения, передача энергии от частицы A частице O и т.п. Соответственно можно определить величины эффективных сечений, например, для потери частицей A определенной энергии от ε1 до ε2, отклонения на угол, больший данного и т.п. Определяются и дифференциальные эффективные сечения, например, дифференциальное эффективное сечение потери энергии частицей A

определяется так, что величина площадки  = f(ε) отвечает потере энергии в интервале от ε до ε +  (при достаточно малой величине ).

Характерные размеры пучка частиц — сантиметры или даже микроны — во много раз превосходят характерные прицельные параметры. Именно поэтому здесь уместен статистический подход. Пусть концентрация газа-мишени составляет nO частиц O на 1 см3 . Тогда в очень тонком слое (толщины dx ) на площадь S приходится nOS dx частиц O, а перекрываемая ими площадь равна dS = σnOSdx. Можно сказать также, что dW = dS/S = σnOdx представляет долю площади, перекрытую частицами O в слое dx. Иначе говоря, это вероятность столкновения частицы в данном слое. Имея в виду возможность выбрать достаточно тонкий слой, мы можем пренебречь возможностью, что какая-то из частиц O будет «затенена» другими. Фактически для этого достаточно условия dW <C 1.

(Для столкновений атомов характерные величины эффективных сечений имеют порядок 10 см2 , для нейтронов и атомных ядер — 10 см2 , на современных ускорителях изучаются сечения еще на 10 - 12 порядков меньшие.)

2. Потери частиц пучка при прохождении слоя

Начнем с простейшей задачи о выбывании частиц из пучка. Пусть любое столкновение ведет к выбыванию частицы — сечение выбывания равно σ . Задача об ослаблении пучка легко решается в таком случае аналитически. Пусть количество частиц A в исходном пучке равно Na0 . После прохождения слоя dx это число уменьшится на dNA = NAdW. Введя число частиц, оставшихся в пучке, Na(x), имеем

dNA(x) = −NA(x)nOσdx.                  (1)

Решение этого дифференциального уравнения, удовлетворяющее условию , т.е.

                               (2)

где b = nOσ. Величина 1/b определяет толщину слоя, в котором пучок потеряет значительное число частиц (ослабнет в e =2.7 раза). Приведенное решение задачи относится к средним значениям числа частиц.

Теперь покажем, как решается эта задача методом Монте-Карло. Будем «пропускать» частицу через тонкие слои толщиной dx, пока она не пройдет слой x (или поглотится). «Розыгрыш» для слоя dx состоит в том, что мы выбираем случайное число с равномерным распределением от 0 до 1 (функция rand в MATLAB) и сравниваем с величиной dW; если rand< dW, то частица выбывает, в противном случае перемещается в следующий слой. Пропустив Na частиц через слой толщины x, мы можем получить число оставшихся частиц Na(x).

Приведем отрезок программы, в котором через слой толщины Xmax «пропускается» Na частиц. При реализации этой части программы предполагается, что для всех частиц заведены два вектора - вектор координат, достигнутых каждой из частицей к текущему моменту времени, и вектор состоящий из 0 и 1. Каждая 1 соответствует еще не выбывшей частице, а 0 - уже выбывшей из потока.

Следующая программа описывает моделирование выбывания частиц из пучка (листинг А).

% Программа расчета потерь пучка %

%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear;

Na=250; k=(1:Na)'; no=100; SecA=0.05; dx=0.01;

x=zeros(size(k)); % Начальный вектор координат

sc=ones(size(k)); % Начальный вектор-счетчик

dW=no*SecA*dx;          % Вероятность поглощения на шаге dx
hl=line(x,k); % Подготовка Na точек для рисования

set(hl, 'Marker', 'o', 'MarkerSize',3);
axis([0 1 0 Na+1]);          % Масштабирование осей
pause; % Пауза перед запуском основного цикла

% Далее выполняется цикл, на каждом шаге которого

% координата x возрастает на dx до тех пор, пока

% частица не пройдет путь Xmax или не поглотится

while (any(sc)> 0)

 ra=rand(size(k)); % Расчет вероятностей с помощью

% датчика случайных чисел для

% всех частиц

 k1=find(ra-dW< 0);       % Определение номеров частиц

% подлежащих отбраковке

 sc(k1)=0; % Зануление счетчика отбракованной частицы

 x=x+sc*dx; % Продвижение на dx остальных частиц

 set(hl,'XData',x); % Рисование

end; % Конец цикла while

%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

Число частиц Nf, прошедших через достаточно толстый слой, оказывается малым и сильно флуктуирует (изменяется от запуска к запуску). Распределение вероятностей для доли прошедших частиц p = Nf/Naэто распределение Бернулли.

3. Превращение частиц пучка при прохождении слоя

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

Приведем участок программы, моделирующей «судьбу» частиц в этом случае (листинг B).

% Моделирующая часть программы с поглощением A   %
% и превращением A -> B %

% SecA - сечение поглощения A

% SecAB - сечение превращения A в B

% scA - счетчик числа частиц типа A

% scB - счетчик числа частиц типа B

clear; % Очистка рабочей области

Na=25; % Начальное число частиц A

L=5; k=(1:Na)'; no=100;dx=0.01; SecA=0.1; SecAB=0.15;

x=zeros(size(k)); % Начальный вектор координат

scA=ones(size(k)); % Вектор-счетчик частиц A

scB=zeros(size(k)); % Вектор-счетчик частиц B

dWa=no*dx*SecA; % Вероятность поглощения A на шаге dX

dWab=no*dx*SecAB;        % Вероятность A -> B на шаге dX

% Имитация поглощения и превращения частиц

% Пока есть частицы A и не пройден путь L

while (any(scA)> 0 & all(x)< L) ra=rand(size(k));

 ka=find(ra-dWa< 0);       % Номера поглощенных частиц A
 
scA(ka)=0; % Выбывание частиц типа A

 kb=find(dWa< ra & ra< dWa+dWab)

 scB(kb)=scA(kb);          % Превращение A в B
 
scA(kb)=0; % Исчезновение этих A

 x=x+scA*dx+scB*dx;

 pause(1);

end;

Для подобной модели несложно также составить и решить уравнения, описывающие изменения среднего числа частиц A и B, подобно уравнениям (1), (2).

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

Материалы для задания 3

Метод Монте-Карло при нахождении значения определенного интеграла

Во многих исследованиях для построения модели некоторого физического процесса необходимо решать задачи нахождения значений определенных интегралов. Рассмотрим метод статистических испытаний (метод Монте-Карло) для решения указанной задачи.

Рассмотрим метод Монте-Карло применительно к задаче нахождения значения определенного интеграла

.

Используем известный факт, что для нахождения значения определенного интеграла необходимо найти площадь криволинейной трапеции, ограниченной отрезком [a,b] оси OX и графиком функции f(x) на этом отрезке.

Для простоты изложения пусть известна аналитическая запись подынтегральной функции f(x) и задан интервал интегрирования [a,b].

Решение задачи состоит в следующем:

на основании равномерного распределения  случайных величин (например, функция rand в MatLab) случайным образом выбираем координаты точки (x,y),

где

Если

,

то данная точка (x,y) вносит вклад в значение определенного интеграла (такие точки образуют множество M1).

Если

,

то точка (x,y) не засчитывается при нахождении значения определенного интеграла (такие точки образуют множество M2).

Следовательно, множество M1 образуют точки, находящиеся ниже графика функции f(x), множество M2 образуют точки, находящиеся выше графика функции f(x).

Обозначим, S1 – площадь криволинейной трапеции, т.е. площадь фигуры, расположенной ниже графика f(x), S2 – площадь фигуры, расположенной над графиком функции f(x).

Рис. 1. Фигуры S1 и S2, ограниченные графиком функции f(x)

Пусть после проведения серии N испытаний в множестве M1 имеем N1 точек, а в множестве M2 имеем N2 точек.

Тогда, отношение площадей фигур S1 и S2 и отношение количества точек N1 и N2 связаны следующим равенством

.

Сумма площадей S1 и S2 равняется площади прямоугольника S0 со сторонами (b-a) и  ymax.

Следовательно, так как

и

то

,

,

.     (3)

Формула (3) позволяет находить значение определенного интеграла методом Монте-Карло.

Материалы для задания 4

Маятник

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

Такое движение маятника хорошо известно — это гармонические колебания. Закон движения можно записать в виде x = acos(ω0t +0), где ω0 – частота колебаний, a — амплитуда, 0  — начальная фаза. (Угол отклонения маятника мы обозначили здесь x). Малые колебания описываются уравнением

линейным относительно функции x, поэтому их обычно называют линейными.

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

Свободные колебания

Результаты исследования движения маятника удобно представить в виде набора кривых на плоскости (x,p), где x - угол отклонения маятника, p = скорость изменения угла.

Плоскость (x,p) называется фазовой плоскостью, переменная p — импульсом, а кривые, определяемые параметрически законом движения как x = x(t), p = p(t), — фазовыми траекториями.

Фазовая траектория определяется, например, начальными значениями координаты x(0) и импульса p(0).

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

Угол отклонения маятника достаточно задавать в некоторых конечных пределах, например, принимая − π ≤ x < π.

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

При некоторой энергии колебания сменяются вращением.

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

       (4)

Имеется аналитическое решение этого уравнения (довольно сложное), однако мы будем исследовать движение маятника численно.

Запишем это дифференциальное уравнение в виде системы уравнений первого порядка

,

.    (5)

Основная идея численного расчета чрезвычайно проста (метод Эйлера): зная значения координаты и скорости в момент времени t, можно приближенно найти их значения через малый промежуток ∆t

,

,

взяв значения величин х и р из уравнений (5). Многократно повторяя такие вычисления, можно найти зависимости x(t) и p(t). Большей точности можно достигнуть, практически не усложняя расчеты, если использовать так называемую вычислительную схему с перешагиванием. В этой схеме вычисляются значения координат в моменты времени

t−1/2∆t,     t + 1/2∆t,     t + 3/2∆t,     t + 5/2∆t,...,

а значения скоростей в моменты времени

t, t + ∆t,      t + 2∆t,      t + 3∆t,...

Данная вычислительная схема и использована в предлагаемой ниже программе.

Чтобы отклонение от точного решения уравнения (1) было небольшим, должен быть достаточно малым шаг по времени ∆t. Точность вычислений можно проверить, применяя метод повторного счета с уменьшенным шагом ∆t. Если при этом решение остается прежним, значит, шаг был выбран достаточно малым и можно быть уверенным в правильности результата.

Далее приводится текст простейшей программы на языке MATLAB.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Фазовая траектория математического маятника     %

clear; % Очистка рабочей области

% Задание начальных значений

x1=2.2; % Начальная координата

p1=0.0; % Начальный импульс

dt=0.025; % Шаг по времени

axis([-pi pi -pi pi]); % Задание диапазона осей

hl=line(x1,p1); % Задание дескриптора линии

% Задание параметров выводимой линии

% 'EraseMode',’none' - это режим вывода без стирания

% 'LineStyle',':'  - вывод линии в виде пунктира

% по умолчанию выводится сплошная линия

% 'Color',’r’      - задание цвета линии

set(hl,'EraseMode', 'none', 'LineStyle',':','Color', 'r');

grid on; % Задание вывода координатной сетки

pause; % Пауза, обеспечивает немедленный вывод

% рисунка на экран для продолжения программы
% необходимо нажать любую клавишу
while 1 % Бесконечный цикл

% Основной алгоритм расчета

x2=x1+p1*dt;

p2=p1-sin(x2)*dt;

% Склейка граничных условий

if x2> pi

x2=x2-2*pi; end; if x2< -pi

x2=x2+2*pi; end;

% Вывод очередного участка фазовой траектории

set(hl,'XData',x2,'YData',p2);

% Переприсвоение начальных значений

x1=x2;p1=p2;

end;

(Здесь шаг по времени ∆t обозначен dt, pi=π.)

Следует сделать несколько замечаний. После запуска программа откроет графическое окно и нарисует в нем одно окно с системой координат, осями и координатной сеткой, после чего перейдет в режим ожидания (команда pause). Для продолжения расчета необходимо нажать любую клавишу, и программа начнет вывод фазовой траектории. Поскольку используется бесконечный цикл, выполнение программы можно прервать клавишей Ctr+C. Для повторного прогона с другими параметрами необходимо ввести соответствующие изменения в текст программы и перезапустить ее.

Задания лабораторной работы

Написать файл-функции для решения поставленных далее задач.

Задачи:

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

2. Построить модель (используя метод Монте-Карло) прохождения частицами пучка некоторого слоя заданной ширины и имеющим заданное эффективное сечение поглощения с учетом превращения частиц А пучка в частицы В (эффективное сечение превращения задано). Условия выполнения приведены далее.

В задачах 1 и 2 для каждого варианта определены следующие значения:

- ширина слоя Xmax,

- количество частиц в исходном пучке Na,

- ширина элементарного слоя dx,

- концентрация частиц слоя (преграды) no,

- эффективное сечение поглощения (A) SecA,

- эффективное сечение превращения (AB) SecB.

Используя текст листингов А и В, необходимо выполнить следующие задания:

- В задачах 1 и 2 необходимо визуализировать процесс прохождения частицами слоя и превращения частиц А пучка в частицы В.

- В задаче 2 необходимо визуализировать процесс превращения частиц А пучка в частицы В.

- Ширина слоя в задачах 1 и 2 не должна превышать величины Xmax.

3. Найти значение определенного интеграла методом Монте-Карло.

Построить график функции f(x) и отметить на графике «случайные» точки, образованные при решении.

4. Построить фазовую траекторию математического маятника. Визуализируйте процесс движения самого математического маятника.

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

Варианты:

Значения переменных для задания 1, 2.

Вариант

Xmax

Na

dx

no

A

AB

1

4

30

0.01

100

0.05

0.1

2

5

40

0.01

100

0.07

0.08

3

6

45

0.01

100

0.08

0.17

4

7

35

0.01

100

0.09

0.15

5

8

30

0.01

100

0.1

0.2

6

7

35

0.01

100

0.05

0.15

7

6

40

0.01

100

0.1

0.25

8

5

30

0.01

100

0.07

0.3

9

4

45

0.01

100

0.01

0.15

10

7

35

0.01

100

0.15

0.2

Значения переменных для задания 3.

Для всех вариантов – количество испытаний N=2000.

Вариант

f(x)

a

b

ymax

1

sin(x)+2

0

3

max f(x) на [a,b]

2

cos(x)+3

1

5

max f(x) на [a,b]

3

x2-3x+10 

2

6

max f(x) на [a,b]

4

sin(2x)+1

3

7

max f(x) на [a,b]

5

cos(x/2)+5

0

4

max f(x) на [a,b]

6

5x2-8x+14 

1

3

max f(x) на [a,b]

7

sin(4x)+2

2

5

max f(x) на [a,b]

8

cos(x/3)+15

3

7

max f(x) на [a,b]

9

7x2-x+2 

1

5

max f(x) на [a,b]

10

sin(3x)+5

2

8

max f(x) на [a,b]

Значения переменных для задания 4.

Вариант

x1

Импульс p1

dt

1

3

5

0.02

2

2

1

0.02

3

2.5

1

0.02

4

3.1

0.1

0.02

5

3

0

0.02

6

2.5

7

0.02

7

0

3

0.02

8

1

4

0.02

9

0.5

5

0.02

10

3

0

0.02

PAGE  1


 

А также другие работы, которые могут Вас заинтересовать

19424. Понятие алгоритма. Исполнитель алгоритма. Система команд исполнителя (на примере учебного исполнителя) 70 KB
  Понятие алгоритма. Исполнитель алгоритма. Система команд исполнителя на примере учебного исполнителя. Свойства алгоритма. Способы записи алгоритмов; блоксхемы. Появление алгоритмов связывают с зарождением математики. Более 1000 лет назад в 825 году ученый из города Хор
19425. Основные алгоритмические структуры: следование, ветвление, цикл; изображение на блок-схемах 87.5 KB
  Основные алгоритмические структуры: следование ветвление цикл; изображение на блоксхемах. Разбиение задачи на подзадачи. Вспомогательные алгоритмы. Основные виды алгоритмов алгоритмических структур: 1. Линейный алгоритм еще называют следование; 2. Циклический а
19426. Величины: константы, переменные, типы величин. Присваивание, ввод и вывод величин. Линейные алгоритмы работы с величинами 62.5 KB
  Величины: константы переменные типы величин. Присваивание ввод и вывод величин. Линейные алгоритмы работы с величинами. Вам уже известно что всякий алгоритм составляется для конкретного исполнителя. Сейчас в качестве исполнителя мы будем рассматривать компьютер осн...
19427. Логические величины, операции, выражения. Логические выражения в качестве условий в ветвящихся и циклических алгоритмах 44 KB
  Логические величины операции выражения. Логические выражения в качестве условий в ветвящихся и циклических алгоритмах. Для того чтобы понять работу ветвящихся и циклических алгоритмов рассмотрим понятие логического выражения. В некоторых случаях выбор варианта де...
19428. Представление о программировании: язык программирования (на примере одного из языков высокого уровня) 32 KB
  Представление о программировании: язык программирования на примере одного из языков высокого уровня; примеры несложных программ с линейной ветвящейся и циклической структурой. Для представления алгоритма в виде понятном компьютеру служат языки программирования. С
19429. Основные компоненты компьютера, их функциональное назначение и принципы работы. Программный принцип работы компьютера 306 KB
  Основные компоненты компьютера их функциональное назначение и принципы работы. Программный принцип работы компьютера. С давних времен люди стремились облегчить свой труд. С этой целью создавались различные машины и механизмы усиливающие физические возможности челов...
19430. Программное обеспечение компьютера, состав и структура. Назначение операционной системы. Командное взаимодействие пользователя с компьютером 673 KB
  Программное обеспечение компьютера состав и структура. Назначение операционной системы. Командное взаимодействие пользователя с компьютером. Графический пользовательский интерфейс. В 5060е годы когда компьютер еще назывался ЭВМ электронновычислительная машина он...
19431. Понятие файла и файловой системы организации данных (папка, иерархическая структура, имя файла, тип файла, параметры файла) 76 KB
  Понятие файла и файловой системы организации данных папка иерархическая структура имя файла тип файла параметры файла. Основные операции с файлами и папками выполняемые пользователем. Понятие об архивировании и защите от вирусов. Все программы и данные хранятся в д...
19432. Информационные ресурсы общества. Основы информационной безопасности, этики и права 60.5 KB
  Информационные ресурсы общества. Основы информационной безопасности этики и права. Информационные ресурсы. Ресурс – это запас или источник некоторых средств. Традиционно различают следующие виды общественных ресурсов: материальные энергетические трудовые финанс