Заметки по R: Дифференциальные уравнения

library("knitr")
opts_chunk$set(
  # cache=FALSE,
               message=FALSE, warning=FALSE) 

library("ggplot2") # для построения графиков
library("rasterVis")
library("fields")
library("deSolve") # решение дифф. уравнений с начальными условиями
library("bvpSolve") # решение дифф. уравнений с краевыми условиями
library("dplyr") # манипуляции с данными

Пакет 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 <- mutate(df, y1dot = y2, y2dot = y1 + cos(y2))

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

df <- mutate(df, len = sqrt(y1dot^2 + y2dot^2), 
             angle = atan2(y1dot, 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 <- mutate(df, y1dot = y2, y2dot = y1 + cos(y2))
plot(df$y1, df$y2, pch = ".", xlab = expression(paste(y[1])),
  ylab = expression(paste(y[2])), main = "График векторного поля")
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.000000 4.000000
## 2 0.01 1.040018 4.003678
## 3 0.02 1.080076 4.007785
## 4 0.03 1.120176 4.012326
## 5 0.04 1.160324 4.017305
## 6 0.05 1.200524 4.022725

str(sol)

## 'data.frame':    1001 obs. of  3 variables:
##  $ time: num  0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 ...
##  $ y1  : num  1 1.04 1.08 1.12 1.16 ...
##  $ y2  : num  4 4 4.01 4.01 4.02 ...

ggplot(sol) + geom_line(aes(time, y1), size = 2) + labs(x = "t",
  y = expression(paste(y[1])), title = "Решение ОДУ с начальными условиями")

Функция 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.0000000 -1.553150
## 2 0.01 0.9845193 -1.543001
## 3 0.02 0.9691398 -1.532904
## 4 0.03 0.9538610 -1.522860
## 5 0.04 0.9386824 -1.512868
## 6 0.05 0.9236035 -1.502928

ggplot(sol) + geom_line(aes(x, y1), size = 2) +  labs(x = "x", 
  y = expression(paste(y[1])), title = "Решение ОДУ с краевыми условиями")

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

Есть несколько способов представить себе функцию от двух переменных, \(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 <- mutate(df, z = (3 * x^2 + y) * exp(- x^2 - y^2))

r <- rasterFromXYZ(df, crs = proj)

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

contour(r)

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

streamplot(r)

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

vectorplot(r)