Решаем дифференциальные уравнения
Прежде, чем приступить к практической части, следует сказать, что для решения ОДУ используется функция ode из R-пакета deSolve, предназначенная для решения обыкновенных дифференциальных уравнений и включающая в себя различные методы решения (например, метод Рунге-Кутта).
Описание функции ode: ode(y, times, func, parms, method = c(lsoda,...,ode45,...))
y - начальное значение искомой функции
times - интервал значений независимой переменной
func - функция, описывающая дифференциальное уравнение (определяется как function(times, y, parms))
parms - параметры (вектор параметров) или значение NULL
method - численный метод (по умолчанию lsoda)
Для того, чтобы использовать эту функцию, дифференциальное уравнение (система дифференциальных уравнений) определяется как R-функция (func), которая вычисляет производные в системе ОДУ, в соответствии с независимой переменной t. Возвращаемым значением func должен быть список, где первый элемент представляет собой вектор, содержащий производные от y по t, а следующий - дополнительные значения, записанные для каждого момента времени. Производные должны быть указаны в том же порядке, как и компоненты вектор-функции y = {y[1], y[2], ..., y[n]}.
Пример 1. Требуется найти решение обыкновенного дифференциального уравнения y ′ = y 2 + t (4) с начальным условием y(0) = 0.1 на интервале t ∈ [0, 1].
Решение. Для того, чтобы получить численное решение данного уравнения в R, необходимо загрузить пакет deSolve (строка 1, Listing1). Простое дифференциальное уравнение (4) реализуется в R с помощью функции f1 (строка 2), аргументами которой являются независимая переменная t, зависимая 8 переменная y и возможные параметры parms, которые в данном случае отсутствуют. Функция f1 возвращает производные как список (строка 4). Отметим, что переменная dY.dt есть обозначение производной y ′ (t) в коде программы. Далее, определяем начальные условия y(0) = 0.1 и интервал независимой переменной t ∈ [0, 1], как последовательность точек от 0 до 1 с шагом 0.1 (строка 5). Решение задачи в R дает функция ode (строка 6). Ответ получаем в виде таблицы ODE, которая состоит из двух столбцов (t; y):
f1 = function(t ,y , parms ){
dY.dt = y^2 + t
return (list(dY.dt))};
y0 = 0.1; t0 = seq(0, 1, 0.1)
sol = ode(y = y0, t = t0, func = f1, parms = NULL )
plot(sol[,1], sol[,2], type ="l", xlab ="t", ylab ="y(t)" ,
col ="red", lwd =4)
ODE = data.frame(t = sol[,1], y = sol[,2])
ODE
> ODE t y 1 0.0 0.1000000 2 0.1 0.1060456 3 0.2 0.1223322 4 0.3 0.1491583 5 0.4 0.1869647 6 0.5 0.2364260 7 0.6 0.2985538 8 0.7 0.3748554 9 0.8 0.4675445 10 0.9 0.5798907 11 1.0 0.7167838 |
Зависимость y от t (рис.1) можно легко изобразить при помощи команды plot (строка 7-8).
Пример 2. Решить уравнение уравнение y'=2*x, начальные условия x0=0; y0=0. Вычислить y(5)
library(deSolve)
#уравнение y'=2*x
f <- function(x, y, params) list(2*x)
x0=0; y0=0; #начальные условия
x <- seq(x0, 5, length.out = 100);
#решение уравнения
df <- as.data.frame(ode(y0, x, f, parms = NULL));
#построение интегральной кривой и начальных условий
plot(NULL,xlim=c(-5,6),ylim=c(0,30),ylab="", xlab="")
lines(df,col="red")
points(x0,y0)
grid(col='black')
95 4.74747475 22.538516483 96 4.79797980 23.020610147 97 4.84848485 23.507805331 98 4.89898990 24.000102036 99 4.94949495 24.497500260 100 5.00000000 25.000000005
Ответ y(5) = 25.000000005
Пример 3. Решим следующее уравнение второго порядка для t ∈ [0, 20]: y ′′ = −0.1y (1) с начальным условием y(0) = 1, y ′ (0) = 0.
Решение. Чтобы решить данную задачу, необходимо, сначала, преобразовать уравнение (1) к системе дифференциальных уравнений первого порядка, с помощью замены y ′ = y1: y ′ = y1 y ′ 1 = −0.1y с начальными условиями y(0) = 1 и y1(0) = 0. Реализация данной задачи в R отображена в скрипте:
f3 = function (t, y, parms ){
with(as.list(y),{
dY.dt = Y1
dY1.dt = -0.1*Y
list(c(dY.dt, dY1.dt))})};
y0 = c(Y = 1, Y1 = 0); t0 = seq(0, 20, 0.1)
out = ode(y = y0, t =t0, f3, parms = NULL )
plot(out[,1], out[,2], type ="l", xlab =" t " , ylab =" y( t )" ,
col ="red", lwd =4)
