Практикум по матлабу. практикум по матлабу. Физических процессов с использованием
 Скачать 1.13 Mb.
|
|
Корреляционной функцией случайных величин x (t), y(t) называют функцию (x(t) − x(t))(y(t) − y(t)). Функцию вида ( 11 ) иногда называют автокорреляционной. В ( 11 ) мы учитываем, что x, v x = 0. 52 Функции x(t)x(t + ξ), x(t)v x (t + ξ) и т.п. зависят не только от ξ, но и от t, и мы их не рассматриваем. Для гармонического осциллятора, в отличие от свободной частицы, эти функции не зависят от t и их интересно изучать наряду с корреляционной функцией скоростей. Задание 4. Получите корреляционную функцию скоростей для броуновских ча- стиц. Для усреднения можно рассматривать одновременно движение многих частиц или суммировать вклады, получаемые за разные интервалы времени для од- ной частицы. Задание 5. Получите для гармонического осциллятора корреляционные функ- ции x(t)x(t + ξ), x(t)v(t + ξ), v(t)v(t + ξ). При этом интересно рассмотреть случаи ωτ 1; ωτ ∼ 1; ωτ 1. 8. Шары В этой работе будем изучать «газ», образованный шарами, которые двигаются без трения по плоскому квадратному столу и абсолютно упруго сталкиваются друг с другом и со стенками («идеальный биллиард»). Едва ли можно наблюдать такой газ где-нибудь, кроме экрана дисплея. Тем не менее, это хорошая модель настоя- щего газа. Пожалуй, самое интересное, что можно увидеть в этой работе, – это возник- новение «молекулярного хаоса» (даже при движении всего лишь двух шаров). Кроме того, можно получать функции распределения по скоростям, наблюдая, как с ростом числа шаров распределения приближаются к максвелловскому. Мож- но изучать смесь газов с частицами разной массы, получать распределение по вы- соте в поле тяжести, получать так называемые стохастический нагрев и стохасти- ческое охлаждение и т.д. 8.1. Расчет движения шаров 8.1.1. Алгоритм расчета Движение шаров за время t рассчитывается довольно просто. Определим момен- ты столкновения для каждой пары шаров без учета возможных помех со стороны других шаров и стенок, а также моменты соударения каждого из шаров с каждой из стенок, неограниченно продолженной, и снова без учета всех помех. 53 Рис. 8. Схема «стола» и системы координат Далее выберем самое раннее из числа этих столкновений (отбросив предвари- тельно те, которые происходили в прошлом). Ясно, что именно это столкновение и произойдет, так как никакое другое не успеет ему помешать. Если промежуток времени до этого столкновения t больше, чем заданный интервал t, то в течение времени t не произойдет никаких столкновений, так что координаты в момент t определяются очевидным образом: r i (t) = r i + v i t. Если же t < t, то нужно рассчитать смещение шаров за время t ,– в результате этого смещения пара шаров (скажем, i- й и j-й) достигают соприкосновения. Затем рассчитываем изменение скоростей этой пары шаров, происходящее в результате столкновения. После этого задача сводится к предыдущей (только «контрольное время» стало меньше t → t − t ): нужно вновь найти ближайшее столкновение, сравнить время до него с новым значением t и т.д. – пока заданный интервал вре- мени не будет исчерпан. Момент столкновения шара, например, с правой стенкой (см. рис. 8 ) опреде- ляется из уравнения x + v x t = L − R, где L– размер стола, R – радиус шара. Изменение скорости при таком столкновении v x → −v x Момент столкновения пары шаров находится из квадратного уравнения (r + vt ) 2 = 4R 2 , где r = r i − r j , v = v i − v j 54 Изменение скоростей при столкновении v i → v i − ∆v, v j → v j + ∆v, где ∆v = n(nv), n = r/(2R) 25 Несложно «включить» также поле тяжести (иначе говоря, наклонить билли- ардный стол). Скажем, если компонента ускорения, создаваемого силой тяжести в направлении против оси Y, равна g, то момент столкновения с «потолком» нахо- дится из уравнения y + v y t − gt 2 / 2 = L − R. 8.1.2. Процедура Balls Описанный алгоритм запрограммирован в виде довольно сложной функции Balls. Для желающих глубже разобраться в алгоритме далее приведен текст процедуры Balls2 на языке MATLAB для 2 шаров. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Функция расчета соударения двух шаров balls2 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function[r,v]=balls2(rin,vin,dt); global rball lx ly ; % Радиус шаров и размеры стола a=rball; a2=4*a*a; eps=1E-10; r=rin; v=vin; % Координаты и скорости шаров ldl=[a a; a a]; % Координаты левого нижнего края % области, доступной центрам шаров, % повторенные дважды - для каждого шара lur=[lx-a ly-a ; lx-a ly-a]; % Правый верхний угол t=dt; % Контрольное время - время до выхода из "balls" while (t > 0) % Учтем возможное наличие очень малых компонент скоростей vm=(v.*v < eps); v1=v+vm; tdl=(ldl-r)./v1; % Моменты столкновения со стенками l и d % учтем столкновения в прошлом и в очень нескорые) 25 Если массы шаров различны, то v i → v i − m j m i + m j ∆v, v j → v j + m i m i + m j ∆v. 55 tdl=tdl+1E10*((tdl < =0)+vm); tur=(lur-r)./v1; % Столкновения со стенками r и u tur=tur+1E10*((tur < =0)+vm); % Выбор самого раннего столкновения со стенкой [t1,j1]=min(tdl(:)); [t2,j2]=min(tur(:)); if(t1 < t2) t0=t1; j0=j1; else t0=t2; j0=j2; end; % Находим момент столкновения шаров друг с другом r0=r(:,1)-r(:,2); v0=v(:,1)-v(:,2); rr=r0’*r0; vv=v0’*v0; rv=r0’*v0; d=rv*rv-(rr-a2)*vv; if (d > 0) & (rv< 0) & (vv> 1E-10) tb=-(sqrt(d)+rv)/vv; else tb=inf; end; % if % Выбор самого раннего соударения if(t < t0)&(t< tb) r=r+v.*t; % сдвиг шаров elseif(t0 < tb) r=r+v.*t0; % и изменение скорости при ударе о стенку v(j0)=-v(j0); else t0=tb; r=r+v.*t0; r0=r(:,1)-r(:,2); v0=v(:,1)-v(:,2); dv=r0’*v0/a2*r0; ddv=[-dv,dv]; % Изменение скоростей при столкновении шаров v=v+ddv; end; % if 56 % После столкновения сделаем маленький сдвиг r=r+v*eps; t=t-t0-eps; end; % while %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Подобная функция для произвольного числа шаров, написанная для ускорения сче- та на языке C, преобразована к виду, допускающему подключение ее к систе- ме Matlab. Она вызывается из программы написанной на языке MATLAB как обычная функция. Ниже приведены необходимые для ее использования сведения. Функция [xf,yf,vxf,vyf]= Balls(n,x,y,vx,vy,dt,R,Lx,Ly,m) по известным значениям координат центров n шаров r i = [x(i), y(i)]) и компо- нент их скоростей v i = [vx(i), vy(i)]), (i = 1, 2, ..., n) определяет значения этих же величин через заданный интервал времени dt и присваивает переменным rf i = [xf(i), yf(i)] vf i = [vxf(i), vyf(i)]). Шары движутся по прямоуголь- ному столу размером Lx × Ly , радиусы шаров одинаковы и равны R, значения масс заданы в массиве m и могут быть различными. 26 Начало координат расположено в левом «нижнем» углу стола, оси координат параллельны его сторонам (см. рис. 8 ). Предусмотрено число шаров n ≤ nmax = 25. Переменные R, Lx, Ly, m можно не задавать. Тогда будут приняты значения R = 20, Lx = Ly = 256, а массы выбраны одинаковыми. Начальные значения координат центров шаров должны быть выбраны в пря- моугольнике R ≤ x(i) ≤ Lx − R, R ≤ y(i) ≤ Ly − R и притом так, чтобы расстояние между любыми центрами было не меньше, чем 2R. 27 Вычисления в процедуре Balls производятся по точным формулам. Программа, которая демонстрирует движение дисков на экране, приведена в файле GAS.m. 26 Во избежание недоразумений отметим, что в процедуре Balls не предусмотрена возмож- ность одновременного столкновения более чем двух шаров, поскольку такое столкновение представляется совершенно невероятным. Если даже будут заданы условия, заведомо при- водящие к тройному соударению, то в процедуре Balls такое соударение будет трактоваться как последовательность парных. 27 В противном случае, т.е. при выборе нефизических начальных условий, возможны неправильные результаты или сбой счета. 57 8.2. Динамический хаос Используя уравнения движения, удается рассчитать с высокой точностью движе- ние планет на тысячи лет вперед (а также и в прошлом). Доказаны теоремы, утвер- ждающие, что движение механической системы, описываемой уравнениями дви- жения, полностью определяется значениями координат и скоростей всех ее точек в некоторый момент времени (теоремы единственности). Казалось бы движение системы 5-10 шаров тоже можно рассчитать хотя бы на сотни «периодов» (пробе- гов от стенки до стенки), тем более, что законы их движения и соударений очень просты. Однако оказывается, что рассчитать движение сталкивающихся шаров на сколько-нибудь длительное время практически невозможно. 8.2.1. Почему движение шаров становится непредсказуемым? Чтобы понять это, представим себе, что направление движения одного из шаров (рис. 9 ) отклонилось от «правильного» на угол ϕ 0 ∼ 10 −8 , и исследуем, как будет изменяться угол отклонения при столкновениях. При этом ограничимся самыми грубыми оценками. Будем иметь в виду, что средний путь шара между столкновениями – l – много больше радиуса шара, l a. За время до очередного столкновения шар движется по прямой и центр его сместится от «правильного» положения на расстояние OO ∼ lϕ 0 . Тогда AA ∼ 1/2 ,OO ∼ lϕ 0 – смещение точки соприкосновения шаров при ударе. Участок поверхности второго шара в окрестности точки касания играет роль «зеркала». Это зеркало двигалось до удара, при ударе оно отскакивает, поэтому направление, в котором отскакивает шар, не определяется правилом «угол падения равен углу от- ражения». Тем не менее поворот «зеркала» на малый угол α ведет к изменению Рис. 9. Схема соударения шаров 58 направления движения отскочившего шара на угол ϕ 1 ∼ α (может быть, в 1.5 - 2 раза больше или меньше – такой уровень точности в данном случае нас устра- ивает). Так как α ∼ AA /a, получаем ϕ 1 ∼ (l/a)ϕ 0 ϕ 0 – при ударе угол отклонения скорости резко увеличивается. После k столкновений ϕ k ∼ (l/a) k ϕ 0 Если l/a ∼ 10, то достаточно 8-10 столкновений, чтобы стало ϕ k ∼ 1 и направле- ние движения шара перестало иметь какое бы то ни было отношение к «правильно- му». Становится очевидным, что никакое разумное повышение точности расчетов не может значительно увеличить правильно рассчитываемый интервал движения шаров. Чтобы предвидеть движение шаров после 100 соударений, согласно этой оценке нужно было бы задать их начальные скорости с фантастической точностью до 100 знаков. Легко понять, что вывод о катастрофическом росте неопределенностей коор- динат относится и к системе большого числа шаров. Относится он и к движению молекул настоящего газа. Только для молекул неопределенности возникают из-за всяческих возмущений, которыми во всех других отношениях можно пренебречь, а для нашего газа шаров – из-за ограниченной точности расчетов. Таким образом, движение шаров (и молекул 28 ) является вполне закономерным за относительно малый промежуток времени и случайным – за долгий промежуток. Отметим, что эта случайность реализуется в рамках закона сохранения энергии. 29 8.2.2. Как убедиться в появлении хаоса? Увидеть, что малое отклонение в направлении движения одного из шаров очень быстро вырастает, очень просто, проведя несколько «запусков» шаров с чуть раз- личными начальными условиями. Чтобы убедиться, что движение шаров спустя достаточный интервал времени становится непредсказуемым, нужно надежно предсказать «силой мысли» какое- то нетривиальное движение. Это совсем не трудно сделать. Представим себе, что в некоторый момент мы заменили направления скоростей всех шаров на противо- положные: v i → −v i . После этого шары должны двигаться строго по прежним траекториям и вернуться в начальное положение. 28 Для молекул есть еще одна причина хаотизации: их движение описывается не классиче- ской механикой, а квантовой, непременно включающей в описание понятие вероятности. 29 Тот факт, что движение молекул газа является хаотическим, понимал и использовал в своих исследованиях еще Л. Больцман. Описанный выше механизм «потери памяти» по- нял и подробно изучил Н.С. Крылов, а исследовал это явление с математической строгостью Я.Г. Синай. 59 Проделать именно такое наблюдение над нашим газом совсем просто. В ре- зультате мы обнаружим, что шары после нескольких столкновений друг с другом «сбиваются с пути». Для качественного исследования достаточно зафиксировать пути частиц на экране ( использовав вместо маркера ’o’ или наряду с ним маркер ’.’ и выбрав для него ва- риант вывода без удаления предыдущего изображения: set(hh, ’EraseMode’,’none’)), а при повторном запуске шаров изменить цвет траекторий. Задание 1. Предусмотрите в программе возможность в некоторый момент tm изменить скорости всех частиц на противоположные: v i → −v i и просле- дите, будут ли шары возвращаться по «проложенным» ранее траекториям. Наблюдайте движение при различных значениях tm. 8.3. Функции распределения Более пригодными, чем зависимости r i (t) для описания движения частиц газа – как настоящих молекул, так и наших «шаров»,– являются средние величины и функции распределения. Будем отмечать на плоскости с координатами (v x , v y ) точки, соответствующие скоростям шаров в разные моменты времени. Эти точки образуют «облако», концентрация точек в котором постепенно убывает к краям. Эта концентрация f 0 (v) становится величиной хорошо определенной, если число N 0 отмеченных точек очень велико. Обычно, функцией распределения по скорости v называют функцию f(v) = f 0 (v)/N 0 , которая в нашем случае удовлетворяет условию ∞ −∞ ∞ −∞ f (v)dv x dv y = 1. Вероятность того, что скорость шара окажется в элементе dv x dv y пространства скоростей вблизи точки v, равна dw = f(v)dv x dv y . (1) Теоретический расчет, который приводить здесь было бы неуместно, 30 дает f (v) = m 2π N − 1 E 1 − ε E N−2 , (2) где ε = mv 2 / 2 – энергия шара; E – суммарная энергия всех N шаров. 30 Этот расчет можно найти в пособии: «Лекции по статистической физике» (Г.Л. Кот- кин,НГУ,1996), §5.1, предназначенном для студентов 3-го курса физфака НГУ. 60 Рис. 10. Схема построения гистограммы модуля скорости При ε E можно заменить 1 − ε/E на exp(−ε/E), так что f = m 2π N − 1 E exp − ε (N − 2) E . Если к тому же N 1, то можно ввести температуру газа: k Б T = E/N, где k Б – постоянная Больцмана. (Для трехмерного газа было бы (3/2)k Б T = E/N). Тогда ( 1 ) сводится к распределению Максвелла: f (v x , v y ) = m 2πk Б T exp − m (v 2 x + v 2 y ) 2k Б T . Если разбить плоскость (v x , v y ) на узкие кольца (рис. 10 ) одинаковой ширины ∆v и подсчитать число точек в каждом из колец ∆N, то можно получить наблюдае- мую функцию распределения (1/N)(∆N/∆v). Соответствующая теоретическая зависимость получается из ( 1 ) заменой dv x dv y на 2πvdv. Подобным же образом можно получить и распределение по энергии. Из равен- ства ε = mv 2 / 2 следует dε = mv dv, поэтому теоретическая зависимость по- лучается заменой 2πv dv на (2π/m)dε. Для получения же наблюдаемой функции (1/N)(∆N/∆ε) нужно разбить плоскость (v x , v y ) на кольца равной площади. Функция распределения по компоненте скорости v x определяется с помощью разбиения плоскости (v x , v y ) на полосы шириной ∆v x (рис. 11 ). Вид функции рас- 61 пределения : dw dv x = ∞ −∞ f (v x , v y )dv y . Интеграл с помощью замены переменной m 2E v y = x 1 − mv 2 x 2E приводится к виду dw dv x = 2m E (N − 1) π 1 − mv 2 x 2E N − 3 2 I (N − 2), где 31 I (N) = 1 0 (1 − x 2 ) N dx = (2N)!! (2N + 1)!! , N !! = N(N − 2)(N − 4)... Задание 2. Получите на экране картину распределения частиц в плоскости (v x , v y ). Для этого следует выводить и сохранять точки с координатами (vx,vy). Ин- тервал времени dt в функции Balls лучше выбрать большим, порядка сред- него времени между столкновениями, чтобы каждый раз на экран выводились новые точки. Задание 3. Получите «наблюдаемые функции распределения» (гистограммы) шаров по компоненте скорости v x , по абсолютной величине скорости v, по энергии. Выведите для сравнения также теоретические функции распределения (полу- чаемые на основе формул ( 2 )). Интересно получить также распределение по скорости движения шаров друг относительно друга, по энергии относительного движения ( m 4 (v 1 − v 2 ) 2 ). 31 Для I (N) с помощью интегрирования по частям легко получить рекуррентное соотноше- ние I (N) = 2N 2N + 1 I (N − 1), сразу же приводящее к ( 31 ) При N 1 основной вклад в I(N) дает, очевидно, область x 1. В этой области 1 − x 2 можно заменить на exp (−x 2 ), после чего интервал интегрирования можно расширить до бесконечности. В итоге I (N) ≈ ∞ 0 exp (−Nx 2 )dx = π/N / 2. 62 Рис. 11. Схема построения гистограммы по проекции скорости v x Задание 4. Получите функцию распределения по расстояниям между центрами шаров (усредненную за длительное время). Укажите существенные отличия таких функций распределения при небольшой концентрации шаров и при условиях, когда среднее расстояние между шарами немногим больше их диаметра. 8.4. Стохастический нагрев и стохастическое охлаждение Если биллиардный стол («ящик с газом») начинает двигаться, очень медленно на- ращивая скорость, а затем так же плавно останавливается, то вместе с ним, есте- ственно, ускоряется и замедляется и газ, сохраняя в конечном счете свою энергию. Если движение «ящика» не является медленным и плавным, то энергия газа в результате такого движения может как вырасти, так и уменьшиться. Какого из- менения энергии газа можно ожидать при длительном периодическом движении «ящика»? В предельном случае, когда газа много, от движущихся стенок распро- страняются звуковые волны, затухающие в объеме газа, поэтому газ нагревается. А если частиц мало? Оказывается, и в этом случае энергия газа будет расти (хотя и не монотонно) – это так называемый стохастический нагрев газа. Процедура Balls определяет движение шаров и в ящике, который движется с 63 постоянной скоростью, естественно, в системе отсчета, связанной с ним. Исполь- зуя эту процедуру, можно найти движение шаров, если скорость «ящика» изменя- ется скачками. Задание 5 Рассмотрите периодическое движение ящика вдоль одной из стенок, при котором ящик половину периода движется со скоорстью U , а второю – со скоростью −U. Постройте график зависимости энергии шаров от времени. Будет ли расти средняя энергия (температура) газа шаров? Для исследования столкновений протонов с антипротонами на ускорителях со встречными пучками необходимо сделать летящие навстречу друг другу сгустки ча- стиц очень концентрированными, а для этого нужно ослабить относительное дви- жение частиц внутри каждого из сгустков. Один из способов сделать это – это так называемое стохастическое охлажде- ние. Можно добиться замедления относительного движения частиц, если по-разному воздействовать на сгусток в зависимости от движения частиц в текущий момент. Такое воздействие достигается благодаря тому, что информация о текущем состо- янии частиц, движущихся по дуге, успевает опередить сгусток, двигаясь по более короткому прямому пути. При этом можно фактически оперировать лишь инфор- мацией, относящейся только к сгустку в целом. Управляя движением «ящика», можно смоделировать стохастическое охлажде- ние. Задание 6 Получите стохастическое охлаждение газа шаров, запрограммировав выбор той или иной величины скорости «ящика» в зависимости от скорости центра масс шаров. 9. Потери пучка при прохождении через вещество В этой работе можно познакомиться с основным методом моделирования, приме- няемым при исследовании прохождения пучков частиц через вещество – методом статистического моделирования (называемым методом Монте-Карло). При этом «судьба» каждой частицы «разыгрывается» с помощью случайного выбора, а полу- ченные для множества частиц результаты подвергаются статистической обработке. Метод применяется, например, при проектировании ядерных реакторов, детекто- ров частиц на ускорителях и обработке получаемых результатов (а также во многих других случаях, скажем, при исследовании распространения мутаций в среде жи- вых организмов). 64 Мы будем изучать, естественно, очень простой вариант задачи – прохождение пучка тяжелых частиц (A) через слой газа, состоящего из легких (O). К тому же будем полагать скорости частиц газа до столкновения с частицами пучка равны- ми нулю. Не будем также учитывать изменений, происходящих в газе вследствие прохождения пучка. Это не ограничивает существенно общности задачи. 9.1. Эффективные сечения Частицы можно представлять шариками радиусов R A и R O . Пусть частица A ле- тит так, что центр ее должен пролететь на расстоянии ρ от центра частицы O. Столкновение произойдет, если ρ < R, где R = R A + R O . Фактически частица O образует для частиц A «преграду», площадь которой σ = πR 2 . Эта величи- на называется эффективным сечением столкновения 32 . В зависимости от величины ρ (называемой прицельным параметром) определяются угол отклонения, передача энергии от частицы A частице O и т.п. Соответственно можно определить величи- ны эффективных сечений, например, для потери частицей A определенной энергии от ε 1 до ε 2 , отклонения на угол, больший данного и т.п. Определяются и диф- ференциальные эффективные сечения, например, дифференциальное эффективное сечение потери энергии частицей A dσ dε = f(ε) определяется так, что величина площадки dσ = f(ε)dε отвечает потере энергии в интервале от ε до ε + dε (при достаточно малой величине dε). Подробнее об этом сказано в работе[ 7 , §18], мы будем использовать взятые оттуда формулы. Характерные размеры пучка частиц – сантиметры или даже микроны – во много раз превосходят характерные прицельные параметры. Именно поэтому здесь уместен статистический подход. Пусть концентрация газа-мишени составляет n O частиц O на 1 см 3 . Тогда в очень тонком слое (толщины dx ) на площадь S при- ходится n O S dx частиц O, а перекрываемая ими площадь равна dS = σn O Sdx. Можно сказать также, что dW = dS/S = σn O dx представляет долю площади, перекрытую частицами O в слое dx. Иначе говоря, это вероятность столкновения частицы в данном слое. Имея в виду возможность выбрать достаточно тонкий слой, мы можем пренебречь возможностью, что какая-то из частиц O будет «затенена» другими. Фактически для этого достаточно условия dW 1. 32 Для столкновений атомов характерные величины эффективных сечений имеют порядок 10 −16 см 2 , для нейтронов и атомных ядер – 10 −26 см 2 , на современных ускорителях изуча- ются сечения еще на 10 - 12 порядков меньшие. 65 9.2. Потери частиц пучка при прохождении слоя Начнем с простейшей задачи о выбывании частиц из пучка. Пусть любое столк- новение ведет к выбыванию частицы – сечение выбывания равно σ. Задача об ослаблении пучка легко решается в таком случае аналитически. Пусть количество частиц A в исходном пучке равно N A0 . После прохождения слоя dx это число уменьшится на dN A = N A dW . Введя число частиц, оставшихся в пучке, N A (x), имеем dN A (x) = −N A (x)n O σdx. (1) Решение этого дифференциального уравнения, удовлетворяющее условию N A (0) = N A0 , N A (x) = N A0 exp(−bx), (2) где b = n O σ. Величина 1/b определяет толщину слоя, в котором пучок потеряет значительное число частиц (ослабнет в e =2.7 раза). Приведенное решение задачи относится к средним значениям числа частиц. Теперь покажем, как решается эта задача методом Монте-Карло. Будем «про- пускать» частицу через тонкие слои толщиной dx, пока она не пройдет слой x (или поглотится). «Розыгрыш» для слоя dx состоит в том, что мы выбираем случайное число с равномерным распределением от 0 до 1 (функция rand в MATLAB) и сравниваем с величиной dW ; если rand< dW , то частица выбывает, в противном случае перемещается в следующий слой. Пропустив N a частиц через слой толщины x, мы можем получить число оставшихся частиц N a (x). Приведем отрезок программы, в котором через слой толщины Xmax «пропус- кается» N a частиц. При реализации этой части программы предполагается, что для всех частиц заведены два вектора -вектор координат, достигнутых каждой из частицей к текущему моменту времени, и вектор состоящий из 0 и 1. Каждая 1 соответствует еще не выбывшей частице, а 0 -уже выбывшей из потока. %%%%%%%%%%%%%%%%%%%%%%%%%%% % Программа расчета потерь пучка % %%%%%%%%%%%%%%%%%%%%%%%%%%% clear; Na=250; k=(1:Na)’; no=100; Section=0.05; dx=0.01; x=zeros(size(k)); % Начальный вектор координат sc=ones(size(k)); % Начальный вектор-счетчик 66 dW=no*Section*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 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Программа Beam1 содержит моделирование выбывания частиц из пучка. Заметим, что отклонения от ( 2 ) делают найденную зависимость более похо- жей на наблюдаемую экспериментально, поскольку здесь моделируется не только среднее, но и наличие флуктуаций относительно среднего. Величина же флуктуа- ций зависит от числа проходящих частиц и тоже может быть определена с помощью моделирования. Число частиц N f , прошедших через достаточно толстый слой, оказывается ма- лым и сильно флуктуирует (изменяется от запуска к запуску). Распределение веро- ятностей для доли прошедших частиц p = N f /N a – это распределение Бернулли. Задание 1. Постройте согласно формуле ( 2 ) зависимость N a (x) в масштабе, удобном для сравнения с модельной зависимостью (гистограммой). Постройте распределение для доли частиц, прошедших толстый слой. На- блюдайте изменение этого распределения при изменениях числа частиц в за- пуске. Подобным же образом моделируются более сложные процессы. Пусть, напри- мер, при столкновении с частицей газа частица пучка A может не только погло- титься, но и превратиться в частицу B, летящую в том же направлении, причем 67 эффективное сечение такого превращения равно σ AB , а эффективное сечение по- глощения частиц A равно σ A Приведем участок программы, моделирующей «судьбу» частиц в этом случае. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Моделирующая часть программы с поглощением 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; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 68 Для подобной модели несложно также составить и решить уравнения, описываю- щие изменения среднего числа частиц A и B, подобно уравнениям ( 1 ), ( 2 ). Стоит также отметить близость рассматривавшейся задачи к задаче о радиоак- тивном распаде атомных ядер (для перехода к ней нужно понимать под x время, прошедшее с начала наблюдения). Задание 2. Добавьте в модель поглощение частиц B и обратные переходы B в A. 9.3. Потери энергии Теперь перейдем к моделированию потерь энергии тяжелыми частицами A при прохождении слоя газа, состоящего из легких. Обозначим массу легкой частицы (O) m, массу тяжелой (A) M , ее скорость V , а энергию E (m M). Оценка, приведенная ниже, показывает, что отклонение траектории частицы A от прямой будет невелико даже при значительной потере энергии. Поэтому будем контролировать прохождения частицей слоя газа тем же способом, что и раньше. Покажем, как можно моделировать потери энергии. Обозначим полное сечение соударений σ, а дифференциальное сечение потери энергии dσ/dε = F (E, ε), где E – энергия частицы перед столкновением, а ε – потеря энергии при столкнове- нии. Сначала следует определить, как и в предыдущих случаях, произошло ли на данном участке dX столкновение. Если же оно произошло, то следует разыграть, какова именно величина потери энергии. Эта операция производится методом бра- ковки, описанным в Приложении D ra=rand(size(k)); ka=find(ra-Na*Sect*dX < 0);% Номера рассеявшихся частиц for ii=ka % Цикл по рассеявшимся на этом шаге частицам % Генерация случайной величины потерянной энергии, % распределенной с плотностью вероятности F(E,eps) while 1 % Бесконечный цикл, выход с помощью break eps(ii)=epsmax*rand; % Выбор величины потери энергии if Fmax*rand < F(E(ii),eps(ii)) break; end; % Конец if 69 end; % Конец while E(ii)=E(ii)-eps(ii); % Вычисление потери энергии end; % Конец for Здесь epsmax=ε max – максимально возможное значение потери энергии, а Fmax =F max – число, не меньшее, чем максимум функции F (E, ε). При столкновениях упругих шариков реализуется особый случай: функция F (E, ε) во всем интервале 0 < ε < ε max фактически от ε не зависит (ε max = 4mME/(m+ M ) 2 ≈ 4(m/M)E) (см. [ 7 , §18, задача 2]). В этом случае вместо приведенно- го выше цикла while - end следует сохранить только отмеченную комментарием строчку вычисления случайной потери энергии. При столкновениях заряженных частиц F (E, ε) неограниченно возрастает при ε → 0. Это связано с очень слабым рассеянием при больших прицельных пара- метрах. В таком случае можно ввести добавочное ограничение, выбрав какое-то значение ε min и отказавшись от учета меньших ε. Тогда F max = F (E, ε min ). С точки зрения физического смысла задачи величины ρ и ε должны быть огра- ничены по ряду причин: при больших прицельных параметрах во взаимодействиях заряженных частиц существенно влияние других частиц, небольшие потери энер- гии не заметны на фоне неизбежного разброса энергий начального пучка и т.п. С точки зрения построения моделирующей программы в предлагаемом способе мо- делирования также не обойтись без такого ограничения. Но если окажется, что полная величина потери, которая при ε < ε min могла бы быть мала в сравне- нии с характерной точностью, принятой при анализе распределения по энергиям ε min X max /dX E/l, то малые значения ε заведомо можно не учитывать. В то же время ставить границу ε min слишком низко невыгодно, так как это приведет к замедлению выбора ε: слишком малую долю будет составлять площадь под кривой F (E, ε) от площади прямоугольника ε min < ε < ε max , 0 < F < F max . Задание 3. Определите распределение по энергиям частиц, прошедших слой Xmax. Предлагается два варианта: • частицы -абсолютно упругие шарики; • частицы -заряженные точки, взаимодействующие по закону Кулона U = α/r. Для этого понадобится выражение для дифференциального эффек- тивного сечения потери энергии при столкновении частиц, которое можно взять в [ 7 , §19]: dσ dε = 2πα 2 M mEε 2 при ε < ε max = 4 m M E. 70 Заметим, что в этой задаче, в отличие от предыдущих, моделирование заведомо является самым простым способом исследования. 9.4. Распределение по углам и энергиям Сначала сделаем грубую оценку, показывающую, что даже при значительной по- тере энергии можно пренебречь отклонением направления движения частиц A от первоначального. Оценка основана на том, что при каждом соударении угол от- клонения изменяется мало, причем направление движения может как удаляться от первоначального, так и приближаться к нему. При столкновении легкая частица получает скорость порядка V , т.е. импульс p ≈ mV и энергию порядка ε ≈ mV 2 2 . Частица потеряет энергию порядка пер- воначальной за N ≈ M V 2 2 /ε ≈ M/m соударений. Угол отклонения тяжелой ча- стицы при одном столкновении θ 1 ≈ p/MV ≈ m/M. Отклонения при разных столкновениях происходят в разные стороны по случайным направлениям, поэтому складываются не углы, а их квадраты – в плоскости (V y , V z ) происходит диффу- зия. Тогда угол отклонения за N соударений θ ≈ θ 1 √ N ≈ m/M 1. Задавать направление движения частицы можно полярными углами вектора ее скорости θ и ϕ. Однако удобнее будет ввести углы θ y = θ cos ϕ ≈ V y /V и θ z = θ sin ϕ ≈ V z /V . Чтобы моделировать отклонение направления скорости частицы при многократ- ных столкновениях, следует воспользоваться дифференциальным эффективным се- чением рассеяния на данный угол dσ dθ = f(θ) 33 Поскольку мы принимаем углы θ небольшими, каждый добавочный угол откло- нения θ 1 можно разыгрывать так же, как первое отклонение от первоначального направления вдоль оси X. Кроме того, следует произвести также выбор азимутального угла ϕ нового от- клонения 34 : phi = 2*pi*rand; 33 Другой вариант расчета основан на том, что угол отклонения частицы в лабораторной системе отсчета θ, угол в системе центра масс χ и потеря энергии ε связаны друг с другом простыми соотношениями (см. [ 7 , |