Практикум по матлабу. практикум по матлабу. Физических процессов с использованием
 Скачать 1.13 Mb.
|
|
Для случайных блужданий в направлении оси x вместо ( 12 ) имеем x 2 (t) = 2D t. 44 6.3. Программа, изображающая случайные блуждания В приводимой ниже программе смещения задаются датчиком случайных чисел. В качестве смещения задается величина ∆x = dh ∗ (2.0 ∗ rand − 1.0), где rand – квазислучайное число, вырабатываемое компьютером, dh – параметр случайно- го блуждания, масштаб «шага». При каждом обращении к процедуре rand(m,n) компьютер выдает матрицу m x n случайных чисел, лежащих в интервале от 0 до 1, причем для нашей задачи (и для множества других) выбор этих чисел неотли- чим от случайного, а распределение этих чисел в указанном интервале равномерное. Легко видеть, что распределение величин ∆x окажется равномерным в интервале −dh < ∆x < dh: dw d ∆x = 0 при |∆x| > dh, 1 2dh при |∆x| < dh. (14) Средний квадрат смещения при одном шаге (∆x) 2 = ∞ −∞ (∆x) 2 dw d ∆x d ∆x = dh 2 3 . (15) Получаемые нами смещения на первых одном – двух шагах совсем не похожи на броуновское движение, так как функция распределения еще далека от гауссовой, однако спустя три – пять шагов функция распределения начинает отлично имити- ровать гауссову и смещения становятся практически такими же, как блуждания броуновских частиц. Так и должно получиться согласно теории вероятностей. Предлагаемые ниже задания имеют главной целью иллюстрацию описанных за- кономерностей. Задание 1. Получите на экране картину движения точек, моделирующих слу- чайные блуждания в соответствие с описанной выше функцией распределе- ния частиц по координате x. Для этого воспользуйтесь начальным вариантом программы diffus.m, имеющимся в пакете MPP, который приведен далее. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Стартовая программа для демонстрации случайных % % блужданий % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% n =500; % Число частиц dh =.02; % Параметр случайного распределения % Задание вектор-столбцов координат точек 45 y =1:n; y=y‘; x =zeros(size(y)); h=plot(x,y,’k.’); % Вывод начального положения точек axis([-2 2 0 n+1 ]); % Задание осей % Определение режима перерисовки и размера точек set(h,’EraseMode’,’background’,’MarkerSize’,3); pause % Пауза для вывода графика на экран i=0; % Начальное значение while 1 % Бесконечный цикл i=i+1; x=x+dh*(2*rand(n,1)-1); % Случайные смещения x-координаты % каждой точки % Смена координат точек на рисунке set(h,’XData’,x,’YData’,y,’Color’,’k’) end; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Задание 2. Выведите на экран гистограмму распределения частиц по коорди- нате x. Это можно сделать на том же рисунке, что и вывод самих частиц, или открыть для отрисовки гистограммы отдельное подокно (см. 2.1 ) или, по- дробно, в Дополнении (п. 8 ) описание функции subplot. При этом следует обратить внимание на то, что функция для отрисовки гистограммы (например hist) при использовании в стандартном виде отрисовывают гистограмму сра- зу после вызова, но как все функции верхнего уровня постоянно перерисовы- вают оси и, следовательно, вместо анимации получаются мигающая картина. Для построения нормально работающей динамической картины необходимо использовать функцию hist в виде [n,x]=hist(y,m), что позволяет сначала насчитать параметры гистограммы, а потом, используя функции stairs, line, построить динамическую гистограмму с помощью оператора set, как это де- лалось ранее в п. 2.4.2 . Подробнее функции, используемые при отрисов- ке гистограмм, описаны в Дополнении, п. 8.2 . Если при попытке нарисовать динамическую гистограмму возникнут проблемы, то можно воспользоваться помещенной в директорию MPP функцией Hist_my в качестве образца. Выведите на экран кроме гистограммы «теоретическую» функцию вида ( 7 ). Для отрисовки функции можно либо насчитывать точки функции в виде век- тора и использовать line, либо использовать функцию fplot (Дополнение, 46 п. 8 ). При этом следует выбирать масштабы так, чтобы на экране площадь под этой кривой была равна площади под гистограммой. Площадь под гистограммой равна S гист = (xmax − xmin)N/L, где xmin, xmax -область отрисовки гистограммы по оси x, N -число ча- стиц, а L -число бинов, на которые разбивается ось x. Чтобы согласовать указанным образом масштаб функции ( 7 ), следует изображать функцию w = S гист f (x) в том же масштабе, что и гистограмму: (xmin ≤ x ≤ xmax, 0 ≤ w ≤ wmax). Сохранится ли описанное согласование масштабов, если значительная доля частиц «расползется» за пределы интервала [ xmin, xmax]? Задание 3. Получите на экране графики x н и x 2 н в зависимости от числа шагов. Напишите функцию, реализующую метод наименьших квадратов (см. При- ложение C ), и с ее помощью рассчитайте коэффициенты соответствующих аппроксимирующих прямых. Постройте прямые, наилучшим образом аппрок- симирующие зависимости x н и x 2 н от числа шагов. Обратите внимание, что при небольшом числе частиц начинаются существен- ные отклонения от закона ( 5 ). Попытайтесь сделать соответствующие оценки и проверить их в ходе машинного эксперимента. Задание 4. Получите на экране двумерную картину случайных блужданий ча- стиц, вышедших из одной точки. Если мы разобьем плоскость x, y на кольца одинаковой ширины и подсчитаем количества частиц в каждом из колец, то получим функцию распределения по расстоянию от начала координат R = √ x 2 + y 2 . Постройте функцию распределения по R. Заметим, что предложенный закон блужданий обладает анизотропией. На- пример, после первого шага частицы заполнят квадрат (а не круг, как было при изотропных блужданиях). Однако спустя несколько шагов облако частиц становится изотропным. Это легко проверить, анализируя формулу ( 7 ). Используя ( 7 ), получите и выведите на экран теоретическую функцию рас- пределения по R. 47 Получите на экране зависимость концентрации частиц от R (в виде гисто- граммы). Для этого количество частиц в каждом из колец, найденное с по- мощью процедуры hist, следует разделить на площадь этого кольца и лишь затем воспользоваться процедурой для отрисовки полученной зависимости. 6.4. Броуновские частицы в поле тяжести Если частицы находятся в поле тяжести, то кроме случайных блужданий они со- вершают еще и регулярное движение, направленное вниз. Если частицы к тому же отскакивают от дна сосуда, то в результате конкуренции случайных блужданий и регулярного смещения вниз устанавливается распределение концентрации частиц, быстро убывающее с высотой, называемое распределением Больцмана: n (z) ∝ exp(−mgz/k Б T ). (16) Здесь n (z) -концентрация частиц на высоте z над уровнем дна. Подобное распределение концентрации с высотой можно получить и в нашей компьютерной модели. Выберем теперь интервал времени, отвечающий одному «шагу», заметно б`ольшим: ∆t τ. Можно считать, что за такое время средняя скорость движения вниз v дрейф успеет установиться и будет определяться условием 23 mg ∼ αv дрейф , (17) означающим, что в среднем сила трения компенсирует силу тяжести. Регулярное смещение b частицы за время ∆t оценивается как b ∼ v дрейф ∆t ∼ gτ∆t. (18) Однако величина ∆t должна быть выбрана и не слишком большой: нужно, что- бы регулярное смещение было мало в сравнении с интервалом высоты, на котором заметно изменяется концентрация частиц. Для этого согласно ( 16 ) должно быть выполнено условие mgb k Б T , т.е. gb v 2 T . Умножив обе стороны этого нера- венства на τ ∆t и учитывая ( 18 ), ( 11 ) (с заменой t → ∆t), получаем b 2 a 2 , т.е. регулярное смещение должно быть мало в сравнении со случайным. 23 Для упрощения формул мы не выписываем архимедову выталкивающую силу. 48 Соответствующие численные оценки, например для условий работы, выполня- емой в лабораторном практикуме, предлагаем проделать самостоятельно. Для моделирования такого движения надо ввести на каждом шаге кроме случай- ного еще и постоянное смещение b, направленное вниз, а также обеспечить «упру- гое отражение» частиц от дна. Что касается отражения от «дна», то закон такого отражения не предопреде- лен моделью и его можно выбирать по-разному, нужно только, чтобы частицы не терялись. При этом в слое, прилегающем к поверхности, может оказаться «неесте- ственно» много или мало частиц, однако отступая от «дна», экспоненциальное убы- вание концентрации с высотой воспроизводится хорошо. (Разумеется, такое рас- пределение устанавливается не сразу.) В итоге устанавливается распределение n (x) = C · exp (−x/h). (19) Зависимость h от b и a легко установить, сопоставляя формулы ( 16 ) и ( 19 ) и поль- зуясь оценками для броуновского движения h ∼ k Б T mg ∼ v 2 T τ ∆t gτ ∆t ∼ a 2 b . (20) Задание 5. Получите на экране распределение частиц «в поле тяжести» и соот- ветствующую наблюдаемую функцию распределения. Удобно направить «по- ле тяжести» вдоль оси x, чтобы затем изображения частиц согласовались с изображением гистограммы. Найдите с помощью компьютерного эксперимента коэффициент пропорцио- нальности в формуле ( 20 ). Для этого заметим, что величина h равна «средней высоте столбца частиц», которая определяется соотношением x = ∞ 0 xn (x)dx ∞ 0 n (x)dx = ∞ 0 e −x/h xdx ∞ 0 e −x/h dx = h. Постройте теоретическую кривую вида ( 19 ) и гистограмму, согласовывая их масштабы и используя режим накопления данных в процедуре hist. (Удоб- но вывести на экран также точку, показывающую положение центра тяжести столба частиц, чтобы заметить момент, начиная с которого высота центра тя- жести не будет изменяться регулярно, а будет лишь флуктуировать; тогда и можно будет начать накопление данных.) При этом в процедуре hist по мере 49 накопления данных удобно будет изменять масштаб, согласуя его с полным числом учтенных точек. 7. Броуновское движение В гл. 6 дано качественное описание движения броуновской частицы. В этом описа- нии существенную роль играет время τ , за которое частица «забывает» направле- ние своего движения. Рассматривая положение частицы через интервалы времени ∆t τ, получаем реализацию процесса случайных блужданий, т.е. случай, изу- ченный в гл. 6 7.1. Случайные силы В данной работе предлагается моделировать движение броуновской частицы с мас- штабом времени ∆t τ. В то же время величина ∆t не слишком мала. В наибо- лее простом для понимания случае, движении броуновской частицы в разреженном газе, будем считать число ударов молекул о броуновскую частицу L за время ∆t больш`им, L 1. Изменение импульса частицы за время ∆t за счет взаимодей- ствия с молекулами среды ∆p = m∆v можно представить в виде ∆p + ∆p случ , где среднее значение определяет силу трения: ∆p = f тр ∆t, f тр = −αv, (1) а добавка ∆p случ ответственна за безостановочное броуновское движение. Такое безостановочное движение есть тепловое движение, и скорость его определяется температурой среды: 1 2 m v 2 = 3 2 k Б T. (2) Эта скорость устанавливается в результате «компромисса» между случайными толч- ками, в среднем разгоняющими частицу, и воздействием силы трения, тормозящим ее. Поэтому величина (∆p случ ) 2 должна быть тем больше, чем выше температура и чем больше коэффициент α, определяющий силу трения. Проведем соответствующую этим рассуждениям количественную оценку. Из- менение скорости частицы за время ∆t v → v − α m v∆t + 1 m ∆p случ (3) 50 не должно приводить к изменению v 2 : v 2 → v 2 (1 − ∆t τ ) 2 + 1 m (∆p случ ) 2 , (4) откуда (∆p случ ) 2 = 2mαv 2 ∆t. (5) (Заведомо малое слагаемое с ( ∆t τ ) 2 1 отброшено.) Таким образом, (∆p случ ) 2 = 6αk Б T ∆t. (6) Можно ввести «случайную силу» f случ = ∆p случ ∆t , нужно только иметь в виду, что ее «амплитуда» зависит от ∆t: f 2 случ = 6αk Б T ∆t . (7) Характер зависимости ( 6 ) от ∆t очевиден для случая, когда в качестве среды рас- сматривается разреженный газ. Тогда среднее значение числа ударов молекул за время ∆t: L ∝ ∆t, а флуктуации этого числа, определяющие ∆p случ , пропорцио- нальны √ L ∝ √ ∆t. Соотношение ( 3 ) можно интерпретировать следующим образом: точка, изоб- ражающая состояние броуновской частицы в пространстве скоростей, совершает случайные блуждания в «потенциальном поле» 1 2 αv 2 Описанный подход к изучению движения броуновской частицы, позволяющий продвинуться в области масштабов ∆t τ, принадлежит П.Ланжевену. Смещение частицы за время ∆t τ ∼ m α находится, естественно, как ∆r = v∆t. (8) Подчеркнем различие в характере движения броуновской частицы, рассматрива- емого в разных масштабах времени. При ∆t τ последовательные положения частицы образуют ясно выраженную траекторию, и при уменьшении ∆t движе- ние приближается к равномерному. Если же уменьшить ∆t, не выходя из области ∆t τ, то можно видеть, что каждый «шаг» является результатом нескольких более коротких, но столь же хаотических шагов. Во избежание недоразумений отметим, что раздел теории вероятностей, на- зываемый теорией броуновского движения, рассматривает случайные блуждания, воспроизводящиеся для сколь угодно малых ∆t (что соответствует пределу τ → 0). 51 Таким же образом можно изучать тепловые флуктуации гармонического осцил- лятора. Вводя в выражение ( 3 ) вклад возвращающей силы mω 2 x, получим v → v − α m v ∆t − ω 2 x ∆t + 1 m ∆p случ . (9) Задание 1. Получите на экране траекторию броуновской частицы и «траекто- рию» в пространстве скоростей. Отметьте на траектории другим цветом поло- жения частицы с интервалом времени τ . Как изменяется характер траекторий с изменением τ , L, ∆t? Задание 2. Получите зависимость x (t), v(t) для гармонического осциллятора. Выведите для сравнения одновременно графики для двух одинаковых осцил- ляторов. Задание 3. Рассматривая движение N ∼ 200 броуновских частиц, получите зависимости v 2 (t) и r 2 (t). Определите коэффициент диффузии. 7.2. Корреляционные функции Значения компоненты скорости частицы в моменты t 1 и t 2 , разделенные интерва- лом ξ = t 2 − t 1 , при | ξ | τ статистически независимы: v x (t 1 )v x (t 2 ) = v x (t 1 )v x (t 2 ) = 0. (10) (Перемножаются компоненты скорости одной и той же частицы, усреднение же подразумевается по очень большому числу частиц.) При значениях | ξ |≤ τ ком- поненты скорости не успевают сильно измениться за время ξ; мерой их взаимной зависимости служит корреляционная функция 24 ϕ (ξ) = v x (t)v x (t + ξ). (11) Поскольку мы рассматриваем движение броуновских частиц, статистические свой- ства которого не изменяются со временем (например, температура постоянна), функ- ция ϕ фактически зависит лишь от ξ, а не от моментов t и t + ξ по отдельности. Отсюда следует, в частности, что ϕ (ξ) -четная функция; чтобы проверить это, достаточно заменить в ( 11 ) t на t − ξ. 24 |