Практика: решить и построить дифференциальное уравнение в R с помощью ggplot

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

Задача. Некая материальная точка во время движения меняет свою скорость по закону: y'=-2*x так же известно, что в момент времени x0=-4 материальная точка прошла путь в y0=-7. Постройте график зависимости пути от времени для данной материальной точки, а так же определите какой путь пройдет эта точка за время от 0 до 3?

Решение. В среде программирования R такую задачу решить очень просто. Воспользуемся библиотекой для решения дифференциального уравнения library(deSolve) и так же не забудем и о библиотеке library(ggplot2). deSolve — пакет для решения дифференциальных уравнений в R. Он включает функции для различных типов дифференциальных уравнений, включая обыкновенные дифференциальные уравнения, дифференциальные алгебраические уравнения и частичные дифференциальные уравнения. 

library(deSolve)
library(ggplot2)
options(scipen = 999, digits = 10)
#уравнение y'=-2*x и н.у. x0=-4; y0=-7;
f <- function(x, y, params) list(-2*x)

x0=-4; y0=-7; #начальные условия
x <- seq(x0, 4, length.out = 100);

#решение уравнения
df <- as.data.frame(ode(y0, x, f, parms = NULL));

#построение интегральной кривой и начальных условий
solode = data.frame(x = df$time, y = df$`1`)

ggplot(data = solode, aes(x = x, y = y))+
  geom_line(colour = 'red', linewidth = 2)+
  ggtitle(label = "граФик ДУ",subtitle = "y'= -2*x") +
  theme_linedraw()

Теперь перейдем ко второму вопросу задачи: определите какой путь пройдет эта точка за время от 0 до 3? Собственно данный путь можно вычислить с помощью интеграла, но вот как раз какой именно функции мы не знаем. Можно найти эту функцию вручную проинтегрировав y'=-2*x и зная начальные условия найти определенный интеграл. А можно иначе: провести регрессионный анализ данных полученных в ходе решения дифференциального уравнения.

# определить какой путь пройдет матточка от 0 до 3
md = lm(y ~ poly(x, degree = 2, raw = TRUE),data = solode)
regres = summary(md)
regres$coefficients

> regres$coefficients
 Estimate Std. Error
(Intercept) 8.99999809079585233 0.000000027403016952
poly(x, degree = 2, raw = TRUE)1 -0.00000001383809134 0.000000007831202679
poly(x, degree = 2, raw = TRUE)2 -0.99999999160557063 0.000000003754102943

Округлив получаем (Intercept) 8.99999809079585233 = 9, poly(x, degree = 2, raw = TRUE)2 -0.99999999160557063: X^2 = -1. Отсюда наша искомая функция y(x) = 9 - x^2. Вот теперь можно вычислить интеграл и ответить на второй вопрос задачи:

# Create function
f<-function(x) 9 - x^2
# Integrate from 0 to 1
integrate(f, 0, 3)

> integrate(f, 0, 3)
18 with absolute error < 0.0000000000002

Ответ за время от 0 до 3 материальная точка пройдет путь равный 18.

А для еще большей наглядности изобразим график выделив на нем область интегрирования.

#предел интегрирования
lower = 0
upper = 3
xi = seq(lower, upper, 0.01)
yi = f(xi)
dfi = data.frame(xi, yi)

# # # # # # # # # # # # # # # # # # # # # # # # # # # 
I = integrate(f, lower, upper) 
S = I$value
titl = paste('интеграл =', S)
ggplot(data = solode, aes(x = x, y = y))+
  geom_line(colour = 'red', linewidth = 2)+
  ggtitle(label = "граФик ДУ",subtitle = titl) +
  geom_area(data = dfi, fill = 4, aes(x = xi, y = yi),
            alpha = 0.4,
            color = 1,    # Line color
            lwd = 1,    # Line width
            linetype = 1) + # Line type 
  
  geom_hline(yintercept = 0, lwd = 1.5)+
  geom_vline(xintercept = lower, col='red', lwd = 1)+
  geom_vline(xintercept = upper, col='red', lwd = 1)+
  theme_linedraw()