Как решить эти уравнения ode в R?

Question

Мне нравится решать следующие обыкновенные дифференциальные уравнения в R методом Рунге-кутты.

dy1 <- (-51.33) * ((1-y[2])/y[1]) 

dy2 <- 1.54 * y[1] * (1-y[2]) - 2.14 * y[2] 

Когда y1 становится равным нулю, то dy1 станет бесконечность. Чтобы этого избежать, мне нужно написать R-код, указывающий, когда y[1] станет меньше 0,001, остановите производную y[1] и сохраните нуль. Я вставил код ниже: R

yini <- c(1,0) 

intabs <- function (t, y, parms) { 
    ifelse (y[1] <= 0.01, dy1 <- 0, no) 
    dy1 <- -51.33 * ((1-y[2])/y[1]) 
    dy2 <- 1.54 * y[1] * (1-y[2]) - 2.14 * y[2] 
    list(c(dy1, dy2)) 
} 

times <- seq(from = 0, to = 1, by = .002) 

out <- ode (times = times, y=yini, func = intabs, parms = NULL, method = "rk4") 

head (out, n=50) 

Я использовал ifelse заявление для обозначения y[1] меньше или равно 0,001, а затем сохранить dy1 как ноль. Я получаю результаты. Но, похоже, я сделал некоторую ошибку, которая приводит к отрицательным значениям dy1. Я очень новичок в написании программ. Пожалуйста, помогите, если вы заметили ошибку.

Answer 1

Такой скачок в функции ODE (если правильно закодирован) является ядом для любого контроллера размера шага. Функция ODE должна быть как минимум такой же гладкой, как и порядок ошибок метода интеграции для управления размером шага для работы по назначению. Он все еще может выйти из строя, если ODE является жестким.

Менее инвазивный метод де-выражение образовывать форму единственного числа 1/y является использование

y/(1e-6+y*y) 

который для y>1e-2 даст хорошее приближение 1/y и гладкая везде. Измените константу, чтобы адаптировать регион, в котором он отклоняется от 1/y.