Визуализация и решение обыкновенных дифференциальных уравнений в R

require(knitr)
opts_chunk$set(cache = FALSE, message = FALSE, warning = FALSE)

require(ggplot2)  # для построения графиков
require(rasterVis)
require(fields)
require(deSolve)
require(bvpSolve)

Пакет rasterVis предназначен для изображения данных на реальных географических картах, поэтому там нужно понятие проекции. Мы пока просто введем это шаманское заклинание

proj <- CRS("+proj=longlat +datum=WGS84")

Построим график векторного поля для системы:

\[ \left\{ \begin{array}{l} \dot{y}_1=y_2 \\ \dot{y}_2=y_1+\cos(y_2) \end{array} \right. \]

Задаем решетку, и рассчитываем \( \dot{y}_1 \) и \( \dot{y}_2 \) в точках решетки:

y1 <- seq(-6, 6, 0.05)
y2 <- seq(-6, 6, 0.05)
df <- expand.grid(y1 = y1, y2 = y2)
df$y1dot <- df$y2
df$y2dot <- df$y1 + cos(df$y2)

Рассчитываем длины и углы для стрелочек, помещаем результат в объект Raster.

df$len <- sqrt(df$y1dot^2 + df$y2dot^2)
df$angle <- atan2(df$y1dot, df$y2dot)

df2 <- df[c("y1", "y2", "len", "angle")]

rast <- rasterFromXYZ(df2, crs = proj)

Строим классический график со стрелочками

vectorplot(rast, isField = TRUE)

Строим няку с капельками

streamplot(rast, isField = TRUE)

Простой график можно руками построить без доп. пакетов. При этом нам нужно самостоятельно уменьшить количество стрелочек.

y1 <- seq(-6, 6, 0.5)
y2 <- seq(-6, 6, 0.5)
df <- expand.grid(y1 = y1, y2 = y2)
df$y1dot <- df$y2
df$y2dot <- df$y1 + cos(df$y2)
plot(df$y1, df$y2, pch = ".")
arrow.plot(df$y1, df$y2, df$y1dot, df$y2dot, arrow.ex = 0.03, length = 0.05)

Решим ОДУ с начальным условиями

Решим систему ОДУ с начальными условиями

Описываем саму систему:

eq1 <- function(t, y, parampampam) {
    return(list(c(y[2], y[1] + cos(y[2]))))
}

Начальные условия:

y.start <- c(y1 = 1, y2 = 4)

Точки, в которых компьютер будет считать функцию:

t <- seq(0, 10, by = 0.01)

Решаем

sol <- ode(y = y.start, times = t, func = eq1)
sol <- data.frame(sol)
head(sol)

##   time    y1    y2
## 1 0.00 1.000 4.000
## 2 0.01 1.040 4.004
## 3 0.02 1.080 4.008
## 4 0.03 1.120 4.012
## 5 0.04 1.160 4.017
## 6 0.05 1.201 4.023

qplot(data = sol, time, y1)

Функция ode возвращает матрицу, а для рисования графиков удобнее табличка с данными, data.frame. Строчка sol <- data.frame(sol) переделывает матрицу в таблицу с данными.

Решим систему ОДУ с краевыми условиями

Описываем саму систему:

eq1 <- function(t, y, parampampam) {
    return(list(c(y[2], y[1] + cos(y[2]))))
}

Граничные условия:

y.start <- c(y1 = 1, y2 = NA)
y.final <- c(y1 = 42, y2 = NA)

Точки, в которых компьютер будет считать функцию:

t <- seq(0, 10, by = 0.01)

Решаем

sol <- bvptwp(yini = y.start, yend = y.final, x = t, func = eq1, nmax = 2000)
sol <- data.frame(sol)
head(sol)

##      x     y1     y2
## 1 0.00 1.0000 -1.553
## 2 0.01 0.9845 -1.543
## 3 0.02 0.9691 -1.533
## 4 0.03 0.9539 -1.523
## 5 0.04 0.9387 -1.513
## 6 0.05 0.9236 -1.503

qplot(data = sol, x, y1)

Бесплатное приложение. Изображение функций двух переменных

Есть несколько способов представить себе функцию от двух переменных, \( z(x,y) \):

• 3D график
• Линии уровня
• Векторное поле градиентов функции

Создаем data.frame с декартовым произведением двух векторов

df <- expand.grid(x = seq(-2, 2, 0.01), y = seq(-2, 2, 0.01))

Изобразим функцию \( z(x,y)=(3\cdot x^2+y)\cdot e^{-x^2-y^2} \).

Cоздаем переменную z как функцию от x и y

df$z <- with((3 * x^2 + y) * exp(-x^2 - y^2), data = df)

r <- rasterFromXYZ(df, crs = proj)

Линии уровня функции z

contour(r)

Капельки текущие по градиенту

streamplot(r)

Направление градиентов, заодно вид сбоку для графика функции

vectorplot(r)