Прикладные дифференциальные уравнения

Занятие 1

Илья Щуров

Факультет компьютерных наук, Прикладная математики и информатика, 2021-22 учебный год

Страница курса

Задача 1

Реализовать метод Эйлера.

In [1]:
import numpy as np
In [2]:
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)
In [3]:
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()
Out[3]:
<matplotlib.legend.Legend at 0x7f8e29e66780>

Задача 2

Решить уравнение $$\dot x = x ^ 2$$ с начальным условием $x(0)=1$ на отрезке $[0, 3]$, построить график.

In [4]:
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')

Вопрос: что случилось?

Задача 3: логистический рост популяции

$$\dot x = (1 - x)x$$
In [5]:
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")

Задача 4

Адаптировать функцию euler, чтобы она могла работать с многомерными фазовыми пространствами

In [6]:
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)

Задача 5

Построить и визуализировать решение модели Лотки — Волтерры:

$$\dot x = x - xy,\quad \dot y = -y + xy,$$

где $x$ — число кроликов, $y$ — число лис.

In [7]:
def lotka_volterra(t, x):
    rabbits, foxes = x
    return np.array([rabbits - foxes * rabbits, -foxes + foxes * rabbits])
In [8]:
t, x = euler_multidim(lotka_volterra, 0, np.array([0.2, 0.2]), 15, 5000)
In [9]:
import plotly.graph_objects as go
In [10]:
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"),
        ),
    ),
)