Построение интегральных кривых обыкновенных дифференциальных уравнений в R

Численно решить дифференциальное уравнение и построить интегральные кривые в R можно с помощью пакета deSolve.

Для начала рассмотрим пример построения интегральной кривой для уравнения $y'=y-x^2$ с начальными условиями $y(x_0)=y_0$:

library(deSolve)

#уравнение y'=y-x^2
f <- function(x, y, params) list(y-x^2)

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

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

#построение интегральной кривой и начальных условий
plot(NULL,xlim=c(-5,5),ylim=c(-3,20),ylab="", xlab="")
lines(df,col="red")
points(x0,y0)

Получается следующий график:

Для построения семейства интегральных кривых можно задать массив начальных условий:

library(deSolve)

#уравнение y'=y-x^2
f <- function(x, y, params) list(y-x^2)

x <- seq(-6, 5, length.out = 100);
y0sec <- seq(-1,29,0.3)

plot(NULL,xlim=c(-5,5),ylim=c(-3,20),ylab="", xlab="")
for(y0 in y0sec)
  {  df <- as.data.frame(ode(y0, x, f, parms = NULL));
     lines(df,col="red")
}

Если задавать массив начальных условий таким образом, то в некоторых областях интегральные кривые ложатся слишком плотно, в других наоборот слишком редко. В таких случаях лучше подобрать массив начальных условий так, чтобы график интегральных кривых был более показательным:

library(deSolve)

#уравнение y'=y-x^2
f <- function(x, y, params) list(y-x^2)

x <- seq(-6, 6, length.out = 100);
y0sec <- c(14,19,22.5,24.5,25.5,25.8,25.95,25.99,
           25.995,25.997,25.999,25.9996,26,26.005,
           26.015,26.035,26.077,26.17,26.35,26.65,27,27.5)

plot(NULL,xlim=c(-5,5),ylim=c(-3,20),ylab="", xlab="")
for(y0 in y0sec)
{  df <- as.data.frame(ode(y0, x, f, parms = NULL));
lines(df,col="red")
}

В результате получится следующее:

И если ранее нанести на график изоклины, получим следующее: