Домашнее задание#
Задача Коши для ОДУ первого порядка.#
I. Метод Эйлера для ОДУ первого порядка.#
Рассмотрим уравнение первого порядка
С начальным условием \(u(t=0) = u_0\).
import numpy as np
import matplotlib.pyplot as plt
%matplotlib notebook
def euler_solve(lam, u0, T, dt):
"""Решает $du/dt = \lambda u$ на $0 < t < T$ с $u(t=0) = u0$ при помощи явного метода Эйлера."""
num_steps = int(T/dt)
tt = np.arange(num_steps+1)*dt
y = np.empty(num_steps+1)
y[0] = u0
for k in range(num_steps):
y[k+1] = y[k] + dt*lam*y[k]
return tt, y
lam = -0.5
tt, y = euler_solve(lam, u0=1.0, T=5, dt=0.1/abs(lam))
plt.plot(tt, y, 'o--', label='numeric solution')
plt.plot(tt, np.exp(lam*tt), '-', lw=2, label='ground truth')
plt.legend(loc='best')
plt.grid(True)
Теперь попробуем задать значение шага \(\tau\) (в коде это dt) такое, что \(|\lambda| \tau > 1\).
lam = -0.5
tt, y = euler_solve(lam, u0=1.0, T=12/abs(lam), dt=2.1/abs(lam))
plt.plot(tt, y, 'o--', label='numeric solution')
plt.plot(tt, np.exp(lam*tt), '-', lw=2, label='ground truth')
plt.legend(loc='best')
plt.grid(True)
Задание 1. Неявный метод Эйлера.#
Решите то же самое уравнение
используя невную схему Эйлера. Сравните поведение неявной и явной схем Эйлера. Постройте решение для нескольких значений шага интегрирования, опишите поведение решения при \(\lambda\tau > 2\).
def implicit_euler_solve(lam, u0, T, dt):
"""Решает du/dt = lambda u на 0 < t < T с u(t=0) = u0 при помощи неявного метода Эйлера."""
# YOUR CODE HERE
raise NotImplementedError()
# Нарисовать графики с разными значениями шага, и отдельно при lambda*tau > 2.
lam = -0.5
tt, y = implicit_euler_solve(lam, u0=1.0, T=8/abs(lam), dt=2.1/abs(lam))
plt.plot(tt, y, 'o--', label='numeric solution')
plt.plot(tt, np.exp(lam*tt), '-', lw=2, label='ground truth')
plt.legend(loc='best')
plt.grid(True)
# Для тестирования
lam = -0.5
tt, y = implicit_euler_solve(lam, u0=1.0, T=8/abs(lam), dt=2.1/abs(lam))
assert (y > 0).all()
II. Системы линейных уравнений#
Рассмотрим систему двух уравнений первого порядка.
где \(\mathbf{u}\) есть вектор длины 2, \(A = \mathrm{const}\) - заданная матрица 2\(\times\)2.
Задание 2. Сравнение явной схемы Эйлера и матричной экспоненты.#
Выполните обобщение алгоритма euler_solve для решения систем линейных уравнений первого порядка с матрицей \(A\), не зависящей от времени, используя явную схему Эйлера.
def euler_solve2(a, u0, T, dt):
"""Solve the system du/dt = Au via an explicit Euler scheme.
Parameters
----------
a : ndarray, shape(ndim, ndim)
The matrix of the l.h.s.
u0 : ndarray, shape(ndim,)
Initial condition
T : float
construct the solution for $t\in [0, T)$
dt : float
Integration step size $\tau$
Returns
-------
t : ndarray, shape (n,)
Integration times
y : ndarray, shape (n, ndim)
Solution of the FD system.
y[k, :] is the solution at t[k].
"""
a = np.asarray(a, dtype=float)
u0 = np.asarray(u0, dtype=float)
num_steps = int(T/dt)
tt = np.arange(num_steps+1)*dt
ndim = a.shape[0]
# YOUR CODE HERE
raise NotImplementedError()
Напишите функцию, возвращающую решение задачи Коши для системы уравнений \(du/dt = A u\) с постоянной матрицей \(A\) через матричную экспоненту. (Используйте библиотечную функцию scipy.linalg.expm)
from scipy.linalg import expm
def mat_exp_solve(a, u0, tt):
"""Construct the solution of $du/dt = A u$ with $u(t=0) = u_0$ at times `tt`.
Parameters
----------
a : ndarray, shape (ndim, ndim)
u0 : ndarray, shape (ndim,)
tt : ndarray, shape (n,)
The values of $t$
Return
------
u : ndarray, shape (n, ndim)
u[:, k] is $\exp(t[k] A)$
"""
# YOUR CODE HERE
raise NotImplementedError()
# Solve via Euler's method, compare to the matrix exponential
from scipy.linalg import expm
a = np.array([[-1, 1],
[0, -1]], dtype=float)
t, y = euler_solve2(a, u0=[1, 1], T=3, dt=0.1)
ym = mat_exp_solve(a, [1, 1], t)
plt.plot(t, y[:, 0], 'o', label='0, Euler')
plt.plot(t, ym[:, 0], '-', label='0, expm')
plt.plot(t, y[:, 1], 's', label='1, Euler')
plt.plot(t, ym[:, 1], '-', label='1, expm')
plt.legend(loc='best')
plt.grid(True)
# Сравните здесь метод Эйлера с методом матричной экспоненты
a = np.array([[-1, 1],
[0, -1]], dtype=float)
III. Жесткие системы#
Рассмотрим линейную систему, \(du/dt = Au\), с матрицей правой части
с начальным условием \(u = (1, 0)^T\).
Система называется жесткой, если для всех \((x, \vec{y}(x))\) (на решениях системы) собственные значения матрицы Якоби системы \(\mathbf{J}(x)\) (в данном случае она равна \(A\)) удовлетворяют условиям
Число \(s=\max _i\left|\operatorname{Re} \lambda_i\right| / \min _k\left|\operatorname{Re} \lambda_k\right|\) называется числом жесткости системы.
Пример.
Для решения задачи Коши системы ОДУ используется численный метод Рунге-Кутты, заданный таблицей Бутчера:
Получим для него функцию и условие устойчивости. Вычислим число жёсткости.
Матрица \(\mathbf{J}(x)\) (см. формулу (7.1.3)) постоянна и ее собственные значения \(\lambda_1=-2, \lambda_2=-4, \lambda_3=-800\) лежат на действительной оси. Так как \(z_i=\lambda_i h\), то и функцию устойчивости достаточно исследовать на действительной оси.
Область устойчивости определяется из условия \(|R(z)| \leq 1\). Получаем
Шаг интегрирования, который удовлетворяет всем условиям \(h \in(0,10 / 800]\).
Число жесткости \(s=800 / 2=400\).
Задание 3. Проверка жёсткости системы и неявные методы.#
Найдите собственные значения матрицы \(A\) (используя np.linalg.eigvals) и прокомментируйте, является ли система жесткой.
Решите систему, используя фиксированный шаг
Стабилен ли метод на шаге такого размера?
# YOUR CODE HERE
raise NotImplementedError()
Постройте графики решения системы на интервале \(0 < t < 1\) с начальным условием \(u = (1, 0)^T\) используя функции euler_solve2 и mat_exp_solve. Используйте несколько значений шага, например \(\tau = 4\cdot 10^{-3}\) и \(\tau = 4.5\cdot 10^{-3}\). Прокомментируйте поведение решений.
# YOUR CODE HERE
raise NotImplementedError()
Реализуйте неявную схему Эйлера для системы линейных уравнений первого порядка с постоянными коэффициентами. Заметьте, что на каждом шаге вам необходимо решать систему линейных алгебраических уравнений (используйте np.linalg.solve ).
Сравните решения, полученные явной и неявной схемами Эйлера.
# YOUR CODE HERE
raise NotImplementedError()
Задание 4*. Исследование на A- и L-устойчивость.#
Определить функцию устойчивости и область устойчивости метода Рунге-Кутты заданного таблицей Бутчера:
Исследовать его на A- и L-устойчивость.
Задание 5*. Исследование на жесткость и область устойчивости.#
Для решения задачи Коши системы ОДУ используется численный метод Рунге-Кутты, заданный таблицей Бутчера:
Получите для него функцию и условие устойчивости. Вычислите число жёсткости.
Задача Коши для ОДУ второго порядка.#
Рассмотрим ОДУ второго порядка, описывающее осцилляции маятника
Задание 6. Законы сохранения и решение ОДУ.#
Преобразуйте данное уравнение второго порядка в систему ОДУ первого порядка.
Решите данную систему уравнений, используя явную схему Эйлера на интервале времени не менее десяти периодов осцилляций.
Мы знаем, что в отсутствии трения выполняется закон сохранения энергии:
Постройте зависимость \(E\) от времени для вашего численного решения. Используйте несколько значений шага. Выполняется ли закон сохранения энергии?
# YOUR CODE HERE
raise NotImplementedError()
Реализуйте схему Рунге-Кутта второго порядка. Используте ее для решения того же уравнения с теми же значениями шага \(\tau\). Сравните решения, полученные методом Рунге-Кутта и методом Эйлера на одинаковых промежутках времени. Проверьте закон сохранения энергии. Обсудите.
# YOUR CODE HERE
raise NotImplementedError()
Проверьте соблюдение закона сохранения энергии после выполнения большого количества шагов. Нарисуйте графики решений вместе с графиком точного решения.
# YOUR CODE HERE
raise NotImplementedError()
Решите задачу используя предиктор и корректор. Предиктор - дифференциальное уравнение первого порядка, корректор - закон сохранения энергии. Убедитесь, что при большом количестве шагов закон сохранения энергии теперь сохраняется. Отличается ли полученное решение от точного решения после большого количества шагов?
Задание 7. Методы предиктора и корректора для ОДУ.#
Используте для решения того же уравнения библиотечную функцию scipy.intergrate_solve_ivp.
Сравните результаты с решениями, полученными методами Рунге-Кутта и Эйлера. Проверьте закон сохранения энергии. Обсудите.
# ... ENTER YOUR CODE HERE ...
Используя библиотечные функции, реализуйте метод прогноза-коррекции Адамса и метод Милна. Проверьте закон сохранения энергии для каждого из них.
# ... ENTER YOUR CODE HERE ...
Задание 8. Нелинейное уравнение Пуассона.#
Напишите программу, которая решает нелинейное уравнение Пуассона:
в области \(0<=x<=10\) с начальными условиями \(\phi(0)=0 \quad , \phi^{\prime }(0) = 0.\)
Постройте аналитическое решение с помощью ряда Тейлора и сравните с численным.
Задание 9. Построение интерфейса и фазовые портреты.#
Откройте файл Физическое моделирование.ipynb . В самом низу есть пример и ссылка на много других подобных примеров.
Замените в предыдущей задачи числа и начальные условия на параметры. Постройте интерфейс для её решений и анимируйте фазовый портрет.
Попробуйте, меняя значение производной в нуле, получить решения, при которых \(\phi(10) = 0\)
Задание 10. Построение интерфейса и фазовые портреты-2.#
Постройте пользовательский интерфейс для решения дифференциального уравнения \(\frac{d^2 y}{d t^2}=-\omega^2 y(t)-\beta y^3(t)\) с начальными условиями: \(y(0)=q, \frac{d y(0)}{d t}=r\).
Вынесите в пользовательский интерфейс \(\omega, \beta, q, r\), число шагов интегрирования, максимальное время интегрирования. Предусмотрите, чтобы можно было бы вводить как положительные, так и отрицательные значения \(\beta\). Пользовательский интерфейс должен быть сверстан в два столбца.
Приложение должно отображать \(y(t), d y(\mathrm{t}) / d t\) и зависимость \(d y / d t\) от у (фазовый портрет). Подумайте, как сохранять результаты вычислительных экспериментов для дальнейшего анализа. Анимируйте фазовый портрет.
Дополнительные задачи.#
Задание 11. Движение планет вокруг Солнца.#
Постановка задачи.#
Задача о движении тела в центральном поле тяготения является хорошим примером, демонстрирующим возможности использования ПК для изучения поведения объекта, подчиняющегося некоторым общим физическим законам. Процесс выведения спутника на орбиту обычно разбивается на два этапа. На первом этапе спутник поднимается над атмосферой практически вертикально на некоторую высоту. Затем обычно последняя ступень ракетоносителя придает спутнику необходимую горизонтальную скорость, и далее он двигается по инерции.
Рассмотрим инерционный полет небольшого тела (спутника) около притягивающего центра с большой массой (Земли). Будем интересоваться тем, какие траектории спутника возможны, какой должна быть его минимальная скорость вблизи поверхности Земли, чтобы он, двигаясь по круговой траектории, не упал на Землю (первая космическая скорость), какой должна быть минимальная начальная скорость спутника, чтобы получилась незамкнутая траектория и спутник ушел от Земли (вторая космическая скорость). В численном эксперименте можно также проверить законы Кеплера:
притягивающий центр расположен в одном из фокусов орбиты;
радиус-вектор от Солнца до планеты “заметает” равные площади за равные промежутки времени;
квадраты периодов обращения спутников вокруг притягивающего центра относятся как кубы больших полуосей эллипсов их орбит.
В основу модели мы положим закон всемирного тяготения. Будем считать, что на тело, движение которого рассматривается, действует только сила тяготения и уравнение Ньютона имеет вид:
Здесь m и M — масса спутника и масса притягивающего центра, G — гравитационная постоянная, r — радиус-вектор, задающий положение спутника относительно притягивающего центра, a — ускорение спутника. Хотя закон всемирного тяготения записан для материальных точек, он имеет такую же форму для сферически-симметричных тел. Движение тела под влиянием центральной силы происходит в одной плоскости, положение которой определяется векторами r0 и v0, задающими начальное положение тела и его начальную скорость. Декартову систему координат с началом в центре тяготения и начало отсчета времени выберем так, чтобы движение происходило в плоскости Oxy и в начальный момент скорость тела была перпендикулярна оси x.
Тогда начальные условия можно записать в виде:
Численная модель.#
Численный анализ задачи удобно проводить, используя в качестве единиц измерения характерные масштабы задачи. В качестве единицы длины удобно взять \(x_{0} .\) Если разговор идет о спутнике Земли, то эта величина имеет порядок радиуса Земли \(R\) и равняется \(R+h\), где \(h\) - высота спутника над поверхностью Земли. Всякое расстояние теперь будет задаваться числом, которое показывает, сколько раз в нем укладывается \(x_{0}\). Безразмерное \(x\) будет равняться \(x\), измеренному в метрах, деленному на \(x_{0}\), также измеренному в метрах. Единицу времени удобно построить, используя гравитационную постоянную и характеристики притягивающего центра. Из уравнения легко видеть, что множитель \(G M / r^{2}\) имеет размерность ускорения \(\left(M / c^{2}\right) .\) Вместо расстояния \(r\) возьмем \(x_{0}\) и сформируем выражение с размерностью времени \((c):\) \(\left(G M / x_{0}^{3}\right)^{-1 / 2}\)
. Его и выберем в качестве единицы времени. В качестве единицы скорости тогда естественно взять \(\left(GM / x_{0}\right)^{1 / 2} .\) Измеренные в этих единицах проекции ускорения определяются следующими уравнениями (здесь и далее для безразмерных физических величин использованы те же обозначения, какие использовались для соответствующих размерных величин):
Все физические величины измеряются теперь в относительных единицах и будут одинаковыми для всех систем “спутник — притягивающий центр”. Уменьшилось также число параметров задачи.
Единственный безразмерный параметр Vo, который остался в задаче, показывает, как соотносятся между собой кинетическая и потенциальная энергии спутника в начальный момент. Действительно, кинетическая энергия спутника в начальный момент равна
\(\mathrm{K}=m v_{0}^{2} / 2\), потенциальная энергия для закона тяготения равна П \(=G M m / x_{0}\), и \(\mathrm{V}_{0}^{2}=2 \mathrm{~K} /\) П.
Для нахождения в различные моменты времени проекций скорости спутника и его координат на временной оси выберем дискретные точки \(t_{n,}\) отстоящие друг от друга на малые интервалы \(\Delta t\).
Тогда проекции скорости \(\mathrm{v}_{x}^{(n+1)}\) и \(v_{x}{ }^{(n+1)}\) в момент времени \(t_{n}+1\) будут приближенно (считаем, что ускорение на этом интервале времени не изменилось) представляться выражениями
а координаты в этот момент будем вычислять, как при равномерном движении (опять считая, что интервал времени Δt мал, и скорость в течение него такая, как в конце интервала):
В начальный момент времени проекции скорости и координаты спутника известны:
Система позволяет шаг за шагом, при малом \(\Delta t\), достаточно точно вычислить траекторию спутника и все ее характеристики.
Для классификации траекторий удобно вычислять их эксцентриситет е. В выбранной системе координат для вычисления эксцентриситета необходимо определить координату \(y=p\) точки на траектории, для которой горизонтальная координата совпадает с фокусом эллипса. Тогда е \(=|p-1|\). Если \(e \leq 1\), то траекторией является замкнутая кривая (окружность при \(e=0\) и эллипс при е 0), при \(e=1\) (в численном эксперименте это значение точно получено быть не может) траекторией является парабола, а при e > 1 траектория есть гипербола. Площадь, “заметаемую” радиус-вектором за некоторое время \(T_{k}=t^{(n+k)}-t^{(i)}=k \cdot \Delta t\), можно вычислить следующим образом:
Интересно сравнить эти площади, вычисленные на разных участках траектории. По второму закону Кеплера они должны быть равны.
По третьему закону Кеплера квадраты периодов обращения относятся как кубы больших полуосей орбит, т.е. для каждой траектории величина \(\mathrm{K} 3=T^{2} / a^{3}\) должна оставаться постоянной. Для нахождения длины большой полуоси эллипса необходимо найти отрицательную координату \(x_{\text {мин }}\) при \(y=0 .\) Момент времени \(t_{M}\), когда это произойдет первый раз, соответствует половине периода обращения спутника. Таким образом, величина большой полуоси равна: \(a=(1+\) \(\left.\left|x_{\text {диин }}\right|\right) / 2\), а период \(T=2 t_{M} .\)
Задание 12. Молекулярная динамика.#
Промоделируйте с помощью формул, описанных ниже, движение 100 частиц, кинетическую энергию их задать по распределению Максвелла, а направлений начальных импульсов - равновероятным образом. Для этого можно сгенерировать 300 чисел с одинаковым нормальным распределением - это будут квадраты проекций скоростей. Затем направление вдоль или против осей - определять с вероятностью 1/2.
Опишите, как изменяется характер движения с изменением температуры.
Взаимодействия между частицами (например, атомами или молекулами) будут считаться с помощью классической физики:
Вот эта формула - на деле полный аналог обычного \(\mathrm{F}=\mathrm{ma}\). (реальная «запутанная» формула длиннее, и выглядит менее эстетично, диссипативную составляющая здесь намеренно опущена).
Математически, моделирование методом частиц, представляет из себя решение задачи Коши для уравнений выше. Вычисление, всего этого дела, лучше всего поручить алгоритму Верле.
Алгоритм Верле используется для вычисления следующего местоположения точки по текущему и прошлому, без использования скорости. Формула получается следующим образом. Записывается разложение в ряд Тейлора для следующего и текущего положения:
Сложив эти 2 уравнения, получим:
Алгоритм оптимален по точности и скорости, но ему надо знать два предыдущих положения частицы. Первый шаг оставим методу Эйлера.
Для моделирования идеального газа самым популярным является потенциал Леннарда-Джонса:
\(r_{\min }\) межатомное расстояние (минимум потенциальной энергии), f - сила взаимодействия, r - расстояние между частицами. После обезразмерирования потенциал выглядит так:
Для ускорения расчётов потенциал Леннард-Джонса часто обрывают:
Ещё одним из способов ускорения вычислений является использование сплайнов (плавно соединяемых между собой многочленов). При этом потенциал взаимодействия разбивается на несколько участков, на каждом из которых он приближается простой функцией. Часто используется следующее приближение:
Задание 13. Закон естественного радиоактивного распада.#
Промоделируйте радиоактивный распад группы атомов. Постройте графики количества оставшихся атомов с течением для случаев \(N_0\) = 10, 100, 1000, 1000000. Постройте их также в логарифмическом масштабе по N.
Эксперименты с тяжелыми атомами говорят, что их ядра распадаются, выбрасывая частицы и превращаясь в более легкие. Распад ядра происходит спонтанно и практически не зависит от внешних условий, в которых оно находится. Число распадов, регистрируемых в единицу времени, не зависит от внешних полей, концентрации радиоактивных атомов, химической структуры молекул, в которые входят радиоактивные атомы, а зависит только от числа атомов, ядра которых способны распадаться. Вероятность естественного распада ядра определяется его устройством и не зависит от того, что происходит вне ядра. Эта вероятность может быть описана некоторым числом, характерным для данных ядер. Таково главное свойство естественного радиоактивного распада.
Вопросы для численного эксперимента на математической модели можно поставить следующим образом:
Как меняется со временем число нераспавшихся ядер, если в начальный момент их было \(N_{0}\) ?
Каково среднее время жизни \(T_{1}\) ядра?
Через какое время \(T_{2}\) распадется половина ядер? Ответы на первый и последний вопросы достаточно легко проверить в натурном эксперименте с быстро распадающимися ядрами.
Приступим к построению математической модели. В принятой концепции явления можно сказать, что на малом интервале времени \(\Delta t\) число распавшихся ядер \(-\Delta N=N(t)-N(t+\Delta t)\) пропорционально \(N^{*} \Delta t .\) Кроме того, это число пропорционально числу нераспавшихся ядер в данный момент времени, поскольку каждое ядро независимо от других может испытать распад на рассматриваемом интервале времени Dt.
Учитывая сказанное, для достаточно малого интервала времени можно написать следующее уравнение радиоактивного распада
Коэффициент пропорциональности записан в таком виде, чтобы константа вещества t имела размерность времени. Она определяет, насколько быстро распадаются ядра данного вещества. В теории радиоактивного распада обычно используют другую постоянную \(\lambda=1 / \tau\), называемую радиоактивной постоянной.
Нам кажется использование \(\tau\) более удобным, тем более что оно, как покажет эксперимент, имеет простой физический смысл. Если \(\tau-\) велико, то за время \(\Delta t\) распадется мало ядер, и наоборот, если \(\tau\) - мало, то ядра распадаются быстро. Во втором случае требуется выбирать достаточно маленькое \(\Delta t\), чтобы \(\Delta N\) было малым и уравнение справедливым. Величина \(\tau\) должна быть измерена в натурном эксперименте по подсчету числа распадов в заданном образце на небольшом интервале времени и вычислена по формуле.
Опыт дает очень разные значения этой константы для разных веществ. Так, для радона \(\tau=5,56\) дня, для радия - 2640 лет, для урана - 7,5\bullet10 \(^{9}\) лет. Для полного определения функции \(N(t)\) нужно задать число ядер вещества в начальный момент:
Уравнения определяют математическую модель радиоактивного распада. Теперь в модели нужно перейти к относительным единицам, чтобы уменьшить число параметров задачи и представить результаты численного эксперимента в относительных единицах. В качестве характерного времени процесса можно взять \(\tau\) - оно будет служить новой единицей времени: размерное время \(t\), измеряемое в секундах, представим в виде \(\mathrm{t}=\tau \cdot \tilde{t}, \tilde{t}-6\) езразмерное время.
В качестве характерного числа атомов возьмем их количество в момент \(t=0\), т.е. представим число атомов в произвольный момент времени в виде \(N=N_{0}\)* \(\widetilde{N} .\) Число \(\widetilde{N}\) показывает долю оставшихся ядер от их числа в начальный момент. Подставим размерные величины:
Все входящие в уравнение величины безразмерные, поэтому в дальнейшем знак безразмерной величины “ « писать не будем.
Перейдем к описанию алгоритма вычисления числа ядер, не распавшихся к моменту времени \(t\). Вдоль оси времени будем продвигаться малыми, но конечными шагами величиной \(\Delta t\), обозначая соответствующие моменты времени \(t_{n}=n \Delta t\) и число ядер \(N_{n} .\) Связь между этими числами задается уравнением:
Запись в правой части числа ядер на левом конце временного интервала \(\Delta t\) может привести к заметной ошибке, если на этом интервале происходит заметное изменение числа ядер.
Система разностных уравнений определяет число нераспавшихся ядер в последовательные моменты времени. Обратим внимание на то, что она не содержит характеристик вещества: полученные при численном эксперименте свойства распада будут универсальными. Индивидуальность ядер проявится только через масштаб времени t.
Вернемся к вопросам, поставленным перед численным экспериментом. Для определения зависимости \(N(t)\) по вычисленным значениям \(N_{a}\) нужно построить плавную кривую процесса. Чтобы вычислить среднее время жизни ядра \(T_{1}\), нужно просуммировать время жизни каждого из ядер и разделить его на полное число распавшихся ядер. Время жизни \(t_{n}\) имеют те ядра, которые распались на интервале \(t_{n+1}-t_{n} .\) Их число есть \(N_{n}-N_{n+1} .\) Отсюда получим
Для вычисления времени \(T_{2}\), за которое распадается половина начального числа ядер, нужно найти для определяемой в численном эксперименте функции \(N(t)\) решение уравнения \(N(t)=0,5 .\)
Задание 14. Моделирование броуновского движения.#
Постановка задачи.#
Броуновским движением называют наблюдающееся под микроскопом непрерывное хаотическое движение мелких частиц, взвешенных в жидкости или газе. Молекулы среды находятся в непрерывном тепловом движении и при столкновении с частицей передают ей соответствующий импульс. Если размеры частицы достаточно малы (меньше \(10^{-4}\) см), то удары, получаемые частицей с разных сторон, не компенсируют друг друга в силу их небольшого количества. Под влиянием этих флуктуаций давления на поверхности частицы она приходит в беспорядочное движение, направление и скорость которого меняются с большой частотой (\(10^{12} c^{-1}\)).
Чтобы анализировать это явление, экспериментатор фиксирует положение частицы через равные промежутки времени, например, через 30 сек., и соединяет их прямыми линиями.
В каждый момент времени частица может с одинаковой вероятностью переместиться по любому направлению. Величина этого смещения также случайна и зависит от свойств жидкости, ее температуры и размеров частицы. Для фиксированной системы имеется некоторое среднее за один такт наблюдения смещение l, которое возрастает с ростом T.
Численная модель.#
Положение броуновской частицы будем задавать в декартовой системе координат. В начальный момент пусть частица находится в начале координат. Смещение частицы вдоль каждой из осей является случайным со средним по величине смещением / и равновероятным в обоих направлениях. Это смещение будем формировать с помощью датчика случайных чисел с равномерным распределением от 0 до \(1 .\)
Координаты нового положения частицы в последовательные моменты времени будем определять так: \(x_{n+1}=x_{n}+d x, y_{n+1}=y_{n}+d y .\) Соответствующий момент времени будет \(t_{n}=\Delta t \cdot n\).
Таким образом, можно сказать, что время пропорционально числу шагов по фиксации положения частицы. Расстояние \(L\), на которое удалилась частица от начала координат, очевидно, равно \(\sqrt{x_{n}^{2}+y_{n}^{2}} .\) За интервал времени \(\mathrm{T}=\mathrm{t}_{n}-\mathrm{t}_{\mathrm{m}}\) квадрат смещения броуновской частицы равен \(\mathrm{S} .\)
В численном эксперименте нужно исследовать зависимость \(S\) от \(\mathrm{T} .\) Ниже приведён пример моделирования броуновского движения на плоскости.
import numpy as np
import random
import matplotlib.pyplot as plt
%matplotlib inline
steps = 5000
random.seed(42)
x,y = np.cumsum(np.random.randn(steps)), np.cumsum(np.random.randn(steps))
points = 10
ip = lambda x, steps, points: np.interp(np.arange(steps*points),
np.arange(steps)*points,
x)
X, Y = ip(x, steps, points), ip(y, steps, points)
fig, ax = plt.subplots(1, 1, figsize=(10, 10))
ax.set_title('Brownian Motion')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.plot(X, Y, color='blue',
marker='o', markersize=1)
[<matplotlib.lines.Line2D at 0x7f1c584514d0>]
Задание 15*. Исследование аппроксимации и сходимости разностной схемы.#
Для предложенных разностных схем вычислить порядок аппроксимации, найти точное решение и исследовать его сходимость к решению дифференциальной задачи Коши \(y^{\prime}=\alpha x, y(0)=y_{0}, 0 \leqslant t \leqslant T, \alpha=\) const. У казать дополнительные краевые условия в схемах в) и г) и выяснить их влияние на сходимость и требования к точности их задания.
a) \(\frac{y_{n+1}-y_{n}}{\tau}=\alpha y_{n}\),
б) \(\frac{y_{n+1}-y_{n}}{\tau}=\alpha \frac{y_{n}+y_{n+1}}{2}\),
в) \(\frac{y_{n+1}-y_{n-1}}{2 \tau}=\alpha y_{n}\),
г) \(\frac{y_{n+1}-y_{n}}{\tau}=\frac{\alpha}{2}\left(3 y_{n}-y_{n-1}\right)\).
Решение.
Аппроксимация
a)
б)