Реализовать метод Эйлера.
import numpy as np
def euler(f, t0, x0, T, steps):
"""
Решает уравнение \dot x = f(t, x) с начальным условием
x(t0) = x0
Возвращает np.array([[t0, x0], [t1, x1], ..., [t_steps, x_steps]]),
где t_steps == T
"""
points = [[t0, x0]]
x = x0
delta_t = (T - t0) / steps
for i in range(steps):
t = t0 + delta_t * (i + 1)
x = x + delta_t * f(t, x)
points.append([t, x])
return np.array(points)
def x_identity(t, x):
return x
assert euler(x_identity, 0, 1, 5, 10).shape == (11, 2)
import matplotlib.pyplot as plt
%matplotlib inline
t0 = 0
x0 = 1
T = 2
for steps in [2, 5, 10, 100]:
solution = euler(x_identity, t0, x0, T, steps)
plt.plot(
solution[:, 0], solution[:, 1], label=f"Euler approximation, {steps} steps"
)
t_true = np.linspace(t0, T, 500)
plt.plot(t_true, x0 * np.exp(t_true), label="True solution")
plt.legend()
<matplotlib.legend.Legend at 0x7f8e29e66780>
Решить уравнение $$\dot x = x ^ 2$$ с начальным условием $x(0)=1$ на отрезке $[0, 3]$, построить график.
t0 = 0
x0 = 1
T = 1.1
steps = 200
solution = euler(lambda t, x: x ** 2, t0, x0, T, steps)
plt.plot(solution[:, 0], solution[:, 1], label=f"Euler approximation, {steps} steps")
--------------------------------------------------------------------------- OverflowError Traceback (most recent call last) <ipython-input-4-18f2dd43aa5b> in <module> 3 T = 1.1 4 steps = 200 ----> 5 solution = euler(lambda t, x: x ** 2, t0, x0, T, steps) 6 plt.plot(solution[:, 0], solution[:, 1], label=f"Euler approximation, {steps} steps") <ipython-input-2-c835d1b0d8e4> in euler(f, t0, x0, T, steps) 12 for i in range(steps): 13 t = t0 + delta_t * (i + 1) ---> 14 x = x + delta_t * f(t, x) 15 points.append([t, x]) 16 return np.array(points) <ipython-input-4-18f2dd43aa5b> in <lambda>(t, x) 3 T = 1.1 4 steps = 200 ----> 5 solution = euler(lambda t, x: x ** 2, t0, x0, T, steps) 6 plt.plot(solution[:, 0], solution[:, 1], label=f"Euler approximation, {steps} steps") OverflowError: (34, 'Result too large')
Вопрос: что случилось?
t0 = 0
x0 = 1
T = 10
steps = 200
for x0 in [0, 0.001, 0.01, 0.1, 0.3, 0.5, 0.7, 1, 1.5]:
solution = euler(lambda t, x: (1 - x) * x, t0, x0, T, steps)
plt.plot(solution[:, 0], solution[:, 1], color="C0")
Адаптировать функцию euler
, чтобы она могла работать с многомерными фазовыми пространствами
def euler_multidim(f, t0, x0, T, steps):
"""
Решает уравнение \dot x = f(t, x) с начальным условием
x(t0) = x0
Возвращает tpoints, xpoints
tpoints == [t0, t1, ..., t_steps],
xpoints == [x0, x1, ..., x_steps]
где t_steps == T
"""
tpoints = [t0]
xpoints = [x0]
x = x0
delta_t = (T - t0) / steps
for i in range(steps):
t = t0 + delta_t * (i + 1)
x = x + delta_t * f(t, x)
tpoints.append(t)
xpoints.append(x)
return np.array(tpoints), np.array(xpoints)
Построить и визуализировать решение модели Лотки — Волтерры:
$$\dot x = x - xy,\quad \dot y = -y + xy,$$где $x$ — число кроликов, $y$ — число лис.
def lotka_volterra(t, x):
rabbits, foxes = x
return np.array([rabbits - foxes * rabbits, -foxes + foxes * rabbits])
t, x = euler_multidim(lotka_volterra, 0, np.array([0.2, 0.2]), 15, 5000)
import plotly.graph_objects as go
go.Figure(
go.Scatter3d(x=t, y=x[:, 0], z=x[:, 1], mode="lines"),
layout=go.Layout(
scene=dict(
xaxis=dict(title="time"),
yaxis=dict(title="rabbits"),
zaxis=dict(title="foxes"),
),
),
)
plt.plot(x[:, 0], label="rabbits")
plt.plot(x[:, 1], label="foxes")
plt.legend()
<matplotlib.legend.Legend at 0x7f8e2d68b978>
plt.plot(x[:, 0], x[:, 1])
plt.xlabel("rabbits")
plt.ylabel("foxes")
Text(0, 0.5, 'foxes')
### FROM: https://stackoverflow.com/a/27666700/3025981
import matplotlib.lines as mlines
import matplotlib.patches as mpatches
def add_arrow_to_line2D(
axes, line, arrow_locs=[0.2, 0.4, 0.6, 0.8],
arrowstyle='-|>', arrowsize=1, transform=None):
"""
Add arrows to a matplotlib.lines.Line2D at selected locations.
Parameters:
-----------
axes:
line: Line2D object as returned by plot command
arrow_locs: list of locations where to insert arrows, % of total length
arrowstyle: style of the arrow
arrowsize: size of the arrow
transform: a matplotlib transform instance, default to data coordinates
Returns:
--------
arrows: list of arrows
"""
if not isinstance(line, mlines.Line2D):
raise ValueError("expected a matplotlib.lines.Line2D object")
x, y = line.get_xdata(), line.get_ydata()
arrow_kw = {
"arrowstyle": arrowstyle,
"mutation_scale": 10 * arrowsize,
}
color = line.get_color()
use_multicolor_lines = isinstance(color, np.ndarray)
if use_multicolor_lines:
raise NotImplementedError("multicolor lines not supported")
else:
arrow_kw['color'] = color
linewidth = line.get_linewidth()
if isinstance(linewidth, np.ndarray):
raise NotImplementedError("multiwidth lines not supported")
else:
arrow_kw['linewidth'] = linewidth
if transform is None:
transform = axes.transData
arrows = []
for loc in arrow_locs:
s = np.cumsum(np.sqrt(np.diff(x) ** 2 + np.diff(y) ** 2))
n = np.searchsorted(s, s[-1] * loc)
arrow_tail = (x[n], y[n])
arrow_head = (np.mean(x[n:n + 2]), np.mean(y[n:n + 2]))
p = mpatches.FancyArrowPatch(
arrow_tail, arrow_head, transform=transform,
**arrow_kw)
axes.add_patch(p)
arrows.append(p)
return arrows
### END FROM
ax = plt.gca()
line, = plt.plot(x[:, 0], x[:, 1])
add_arrow_to_line2D(ax, line, arrowsize=2, arrowstyle='->')
plt.xlabel("rabbits")
plt.ylabel("foxes")
Text(0, 0.5, 'foxes')
Изучите с помощью компьютерной симуляции, как зависит ошибка метода Эйлера от длины шага $\Delta t = (T-t_0) / N$, где $N$ — количество шагов? Возьмите уравнение с известным решением (например, $\dot x = x$) и для какого-то фиксированного начального условия и фиксированного $T$ найдите точное решение $x(T)$ и его приближение $x_{Euler}(t; N)$ для $N$ шагов. Постройте график $|x(T) - x_{Euler}(t; N)|$ как функции от $N$ при больших $N$. Умножьте ошибку на $N$ и на $N^2$ и постройте графики получающихся функций. Возьмите логарифм ошибки и поделите его на логарифм $N$, постройте график этой функции при больших $N$. Сделайте вывод о скорости уменьшения ошибки в терминах $O$-больших или $\Theta$ от $N$ и $\Delta t$.
Модифицируйте функцию euler_multidim
таким образом, чтобы она обрабатывала случай T < t_0
.
Найдите все положения равновесия в модели Лотки — Вольтерры. Покажите, что существует единственное положение равновесия в области $x>0$, $y>0$. Изучите с помощью компьютерных экспериментов, как ведут себя решения вблизи этого положения равновесия и вдалеке от него. Является ли на ваш взгляд это положение равновесия устойчивым? Объясните, почему. В каких случаях решения похожи на решения уравнения осциллятора (близки к гармоническим колебаниям)? В каких случаях не похожи? Зависит ли период от выбора начального условия? Как вы думаете, как объясняются эти результаты?