Задания для самостоятельной работы

Федеральное государственное образовательное бюджетное учреждение

Высшего образования

«ФИНАН

«ФИНАНСОВЫЙ УНИВЕРСИТЕТ ПРИ ПРАВИТЕЛЬСТВЕ

РОССИЙСКОЙ ФЕДЕРАЦИИ»

Департамент анализа данных, принятия решений и финансовых технологий

 

С.А.Зададаев

 

Численное решение дифференциальных

Уравнений в R (RStudio)

Учебно-методические рекомендации для проведения

семинара №24 по компьютерному практикуму

 

Для бакалавров направления 38.03.01 «Экономика»

 

 

Электронное издание

 

Москва 2017

Численное решение дифференциальных

Уравнений в R (RStudio)

Обратимся к еще одному применению вычислительных компьютерных средств в решении математических задач. В данном разделе мы познакомимся с общей идеей численных методов решения дифференциальных уравнений на примере схемы Эйлера и ее реализации на языке R.

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

Мы же остановимся лишь на вычислительном аспекте применения схемы Эйлера и прокомментируем основные теоретические проблемы по возможности на простых примерах.

Задача Коши

Задание 1. Изобразить на плоскости семейство интегральных кривых для уравнения  и выделить среди них частное решение, удовлетворяющее начальному условию .

Решение.  Рассмотрим заданное обыкновенное дифференциальное уравнение первого порядка:

.

Общее решение (совокупность всех частных решений) этого уравнения несложно получить путем непосредственного интегрирования его правой и левой частей:

.

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

Код процедуры построения данного рисунка на R здесь несложный:

x <- seq(-5, 5, by = 0.01)                                       # Задаем последовательность x

plot(x,atan(x), type = "l", lwd = 3, col = "red")     # Рисуем arctg(x) жирной линией (lwd = 3)

abline(h = 0, v = 0, col = "gray40")                      # Рисуем оси координат

for (C in seq(-2, 2, by=0.4)) {                               # Выбираем последовательность для C

lines(x,atan(x)+C, type = "l", lwd = 1, col = "red") # Рисуем другие интегральные кривые

}

Заметим, что по сути мы решили задачу Коши:

 

определив её точное единственное решение .

    В следующем задании получим это решение с помощью приближенного метода.

Схема Эйлера

    Идея применения схемы Эйлера численного решения задачи Коши для дифференциального уравнения первого порядка:

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

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

Отметим, что расстояние между соседними точками при этом одинаково и равно:

.

    Вторым шагом записывают разложение точного решения  в ряд Тейлора в точке :

и отбрасывают слагаемые старше первой степени . Здесь-то и происходит возникновение приближенного равенства:

.

Ясно, что чем меньше выбор шага , тем точнее приближение.  Далее из полученного приближения нетрудно выразить производную в i– ом узле сети :

,

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

,

то получится некоторое рекуррентное соотношение между значениями неизвестной функции  в соседних узлах сети:

.

Мы не будем переобозначать неизвестную функцию  на ее приближенный аналог, но со следующей строчки мы вернем точный знак равенства, и это будет верно уже исключительно для приближенного решения, а не точного:

.

Данный тип уравнений называют разностным уравнением. Решается такое уравнение достаточно просто, т.к. для нахождения  необходимо всего лишь знать предыдущее значение . Поэтому сначала задается начальное условие  и с его помощью находят последующее значение :

Далее подставляют найденное значение , как уже известное, в ту же правую часть рекуррентного соотношения и определяют следующее значение и т.д. пока не доберутся до конца отрезка .

Запишем в итоге получившуюся схему Эйлера в терминах – приближенного сеточного аналога точного решения:

Итог по схеме Эйлера

 Если требуется численно решить задачу Коши

                                                                                (1)

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

                                                (2)

Полученная числовая совокупность значений функции  в узлах сети и будет дискретным (разностным) аналогом решения задачи Коши, а если ее соседние точки на графике соединить отрезками прямых, то такое приближенное решение будет уже непрерывным аналогом решения, называемым ломаной Эйлера.

Прокомментируем конкретное применение схемы Эйлера, ее фактическую суть и «подводные камни» на примерах.

 

Задание 2. Решить приближенно с помощью схемы Эйлера задачу Коши

Решение.   Для удобства запишем систему (2) в нашем конкретном случае:

                                   (3)

Теперь нам остается аккуратно составить алгоритм вычислений .

1. Будем искать решение справа от начальной точки , скажем, на отрезке . Введем переменные концов отрезка  и ,  длину шага , само разбиение отрезка с помощью оператора seq и n – общее количество точек разбиения:

xo <- -5                                  # Начальное значение x

xn <- 5                                   # Конечное значение x

dx <- 1.2                                # Величина шага dx

x <- seq(xo, xn, by = dx)          # Разбиение отрезка на xi

n <- length(x)                         # Количество получившихся xi

2. Образуем наш численный аналог решения – последовательность значений функции в узлах сети, первоначально заполнив этот массив (вектор) нулями:

y <- rep(0,n)                            # Заполняем массив y[i] нулями для всех xi

y[1] <- yo <- atan(-5)              # Задаем начальное значение как первое значение y[1]

3. Вычисляем в цикле последующие y[i+1] согласно нашей формуле (3). Обратите внимание, что переменная цикла изменяется не до n, а до (n-1):

for (i in 1:(n-1)) {

y[i+1] <- y[i] + dx/(1 + x[i]^2)

}

4. Последним шагом отображаем полученный результат на графике:

plot(x, y, type = " l ", lwd = 5, col = "blue")     # Рисуем y(x) - численное решение

abline(h = 0, v = 0, col = "gray40")           # Рисуем оси координат

Если при этом мы хотим сравнить полученное приближенное решение с точным , то добавим еще пару строк с построением точного решения поверх предыдущего рисунка:

t <- seq(xo, xn, by = 0.01)                               # Разбиваем вспомогательную ось ox (~t)

lines(t, atan(t), type = " l ", lwd = 4, col = "red") # Рисуем точное решение

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

for (j in c(0.15,1/2,1/3,-1/2,-1/4,-1,1)) {

lines(t, atan(t)-j, type = " l ", lwd = 1, col = "red") # Рисуем интегральные кривые

}

Теперь смотрим:

    Мы специально сделали шаг  – просто огромным, чтобы показать: что происходит в схеме Эйлера. Наше приближенное решение (синий цвет) в каждой итерации, т.е. в каждой точке на каждом изломе принимает вид касательной к графику той интегральной кривой, с которой пересекается в текущем узле сети. Таким образом, первоначально стартует приближенное решение с касательной к настоящему решению, но неизбежно съезжает с него на соседнюю интегральную кривую и в следующей точке уже ориентируется на касательную к другой (пусть даже и очень близкой к исходной) интегральной кривой.

    Это один из самых главных недостатков схемы Эйлера. Нет никакой гарантии, что, съехав на соседнюю интегральную кривую, решение не унесет нас совсем в другую сторону. Конечно, можно сильно уменьшить шаг , но все-равно есть шанс рано или поздно отклониться сколь угодно далеко от истинного решения. Здесь важную роль играет устойчивость самого решения и взаимная конфигурация поля интегральных кривых (см. задание 4).

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

# Пример с y'=1/(1+x^2); y(-5)=arctg(-5) Точное решение y(x) = arctg(x)

 

xo <- -5                                             # Начальное значение x

xn <- 5                                              # Конечное значение x

dx <- 0.01                                         # Величина шага dx

x <- seq(xo, xn, by = dx)                     # Разбиение отрезка на xi

n <- length(x)                                    # Количество получившихся точек разбиения xi

y <- rep(0,n)                                      # Заполняем массив y[i] нулями для всех xi

y[1] <- yo <- atan(-5)                        # Задаем начальное значение как первое значение y[1]

for (i in 1:(n-1)) {

y[i+1] <- y[i] + dx/(1 + x[i]^2)       # Программируем правую часть уравнения

}

plot(x,y, type = " l ", lwd = 5, col = "blue")                  # Рисуем y(x) - численное решение

abline(h = 0, v = 0, col = "gray40")                        # Рисуем оси координат

t <- seq(xo, xn, by = 0.01)

lines(t,atan(t), type = "l", lwd = 1, col = "yellow")      # Рисуем точное решение

Если увеличить масштаб отображения окна (Ctrl + колесо мыши), то можно увидеть «точное» попадание желтой кривой в синий коридор.

Задание 3. Решить приближенно с помощью схемы Эйлера задачу Коши

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

Решение.  Заметим, что само дифференциальное уравнение здесь является линейным неоднородным первого порядка с постоянными коэффициентами и может быть легко решено одним из трех методов: методом Бернулли, методом вариации постоянной или методом неопределенных коэффициентов (проделайте это самостоятельно). Общее решение принимает вид:

,

а подстановка начального условия  приводит к уравнению на C:

,

из которого  и точное решение задачи Коши получается следующим:

.

Теперь решим задачу численно. Будем искать приближенное решение на отрезке, например, . Для удобства запишем систему (2) в нашем конкретном случае:

                                   

и внесем необходимые изменения в итоговый код программы для задания 2 (выделено синим жирным):

# Пример с y'+ y = 1 + x; y(-2)= -2 + exp(2) Точное решение y(x) = x + exp(-x)

 

xo <- -2                                              # Начальное значение x

xn <- 7                                               # Конечное значение x

dx <- 0.01                                         # Величина шага dx

x <- seq(xo, xn, by = dx)                     # Разбиение отрезка на xi

n <- length(x)                                    # Количество получившихся точек разбиения xi

y <- rep(0,n)                                      # Заполняем массив y[i] нулями для всех xi

y[1] <- yo <- -2 + exp(2)                    # Задаем начальное значение как первое значение y[1]

for (i in 1:(n-1)) {

y[i+1] <- y[i] + dx*(1 + x[i] - y[i])  # Программируем правую часть уравнения

}

plot(x,y, type = " l ", lwd = 5, col = "blue")                  # Рисуем y(x) - численное решение

abline(h = 0, v = 0, col = "gray40")                        # Рисуем оси координат

t <- seq(xo, xn, by = 0.01)

lines(t, t+exp(-t), type = " l ", lwd = 1, col = "yellow") # Рисуем точное решение

    В итоге получим прекрасное соответствие точного и приближенного решений для не очень малого шага :

Переходя ко второй части задания (определение минимума полученного решения), заметим, что точка истинного минимума для точного решения задачи Коши

определяется из уравнения

и оказывается равна:

,

что дает значение самого минимума:

.

Замечание. Здесь , так что речь идет, действительно, о точке минимума.

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

# Нахождение минимума y:

min(y)                               # Минимальное значение среди всех y[i]

Но если нам необходимо определить и соответствующее значение самой точки , то необходимо отследить для какого номера i достигается это минимальное значение:

# Нахождение точки минимума и самого минимума:

i <- which.min(y)               # Номер в массиве разбиений отрезка [xo, xn] для минимума

x[i]                                 # Точка минимума (xmin)

y[i]                                 # Значение минимума (ymin). Можно ввести просто min(y)

    Кстати, получающееся значения  и  более, чем убедительно, совпадают с точными:

> x[i]              # Точка минимума (xmin)

[1] -0.01

> y[i]              # Значение минимума (ymin). Можно ввести просто min(y)

[1] 0.9899832

Здесь точность определения точки минимума, конечно, не превышает точности самого шага .

Аналогично можно получить максимальное значение путем замены min на max:

# Нахождение точки максимума и самого максимума:

i <- which.max(y)              # Номер в массиве разбиений отрезка [xo, xn] для максимума

x[i]                                 # Точка максимума (xmax)

y[i]                                  # Значение максимума (ymax). Можно ввести просто max(y)

Замечание. Функции which.max(y)  и which.min(y) в случае нескольких одинаковых экстремумов возвращают все соответствующие номера в виде массива по возрастанию.

                      

Задание 4. Решить приближенно с помощью схемы Эйлера задачу Коши для уравнения Риккати:

и сравнить графически полученный результат с точным решением задачи:

.

      Решение.  Желающие получить общее решение данного уравнения могут воспользоваться тем, что известно одно из его частных решений , и свести уравнение с помощью замены  к уравнению Бернулли.

Мы же, как в предыдущих случаях выберем отрезок и запишем рекуррентную систему для нахождения сеточного решения. Будем рассматривать решение на отрезке . Подставляя в (2) правую часть нашего приведенного уравнения:

,

получим систему:

    Соответственно это приведет к следующим изменениям кода:

# Пример с y'+ y^2 = 1/2 + cos(x) + (x/2+sinx)^2;       y(-4)=sin(-4) - 2 

# Точное решение y(x) = sin(x) + x/2

 

xo <- -4                                              # Начальное значение x

xn <- 7                                               # Конечное значение x

dx <- 0.01                                         # Величина шага dx

x <- seq(xo, xn, by = dx)                     # Разбиение отрезка на xi

n <- length(x)                                    # Количество получившихся точек разбиения xi

y <- rep(0,n)                                      # Заполняем массив y[i] нулями для всех xi

y[1] <- yo <- -2 +sin(-4)                    # Задаем начальное значение как первое значение y[1]

for (i in 1:(n-1)) {

 y[i+1] <- y[i] + dx*(1/2 + cos(x[i]) + (x[i]/2 + sin(x[i]))^2 - y[i]^2)   

}

plot(x,y, type = " l ", lwd = 5, col = "blue")             # Рисуем y(x) - численное решение

abline(h = 0, v = 0, col = "gray40")                       # Рисуем оси координат

t <- seq(xo, xn, by = 0.01)

lines(t, sin(t) + t/2, type = " l ", lwd = 2, col = " red ") # Рисуем точное решение

    Однако результат нас не может порадовать из-за существенного отклонения приближенного решения от истинного:

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

    В последнем, уже удовлетворительном случае, нам пришлось разбить отрезок более, чем на миллион точек. И это не предел для R в резервировании массивов. Однако десять миллионов элементов массива – это уже возможный предел. В этих случаях необходимо оптимизировать вычисления с динамической записью результатов в файл. 

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

    В настоящее время схемой Эйлера практически не пользуются, однако, быстро и просто «прикинуть» регулярное решение ей всегда можно.

Задание 5.* Исследовать приближенное решение, полученное в предыдущем задании, на локальный максимум в окрестности точки  и локальный минимум в окрестности точки .

    Решение.  Обратимся еще раз к графику, полученному в задании 4, итогового приближенного решения задачи Коши:

    Если мы введем код, подобный тому, который мы использовали в задании 3:

# Нахождение точки минимума и самого минимума:

i <- which.min(y)               # Номер в массиве разбиений отрезка [xo, xn] для минимума

x[i]                                 # Точка минимума (xmin)

y[i]                                 # Значение минимума (ymin). Можно ввести просто min(y)

 

# Нахождение точки максимума и самого максимума:

i <- which.max(y)              # Номер в массиве разбиений отрезка [xo, xn] для максимума

x[i]                                 # Точка максимума (xmax)

y[i]                                  # Значение максимума (ymax). Можно ввести просто max(y)

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

> x[i]                  # Точка минимума (xmin)

[1] -2.09615

> y[i]                  # Значение минимума (ymin). Можно ввести просто min(y)

[1] -1.912823

> x[i]                  # Точка максимума (xmax)

[1] 7

> y[i]                  # Значение максимума (ymax). Можно ввести просто max(y)

[1] 4.156987

а нам требуется определить локальные экстремумы в окрестностях точек

и .

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

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

k <- which(x>=4)[1]            # Номер в массиве, до которого ищем максимум

i <- which.max(y[1:k])            # Номер глобального max в массиве разбиений отрезка [xo, 4]

x[i]                                       # Точка максимума (xmax)

y[i]                                       # Значение максимума (ymax). Можно ввести max(y[1:k])

> x[i]                # Точка максимума (xmax)

[1] 2.09279

> y[i]                # Значение максимума (ymax). Можно ввести max(y[1:k])

[1] 1.913587

С другой стороны, из графика нетрудно видеть, что на отрезке  искомый локальный минимум является глобальным. Аналогично ограничиваем   и получаем:

k <- which(x>=3)[1]            # Номер в массиве, начиная с которого ищем минимум

i <- which.min(y[k:n]) + k - 1  # Номер в массиве разбиений отрезка [3, xn]

x[i]                                       # Точка минимума (xmin)

y[i]                                       # Значение минимума (ymin). Можно ввести min(y[k:n])

> x[i]                # Точка минимума (xmin)

[1] 4.18879

> y[i]                # Значение минимума (ymin). Можно ввести min(y[k:n])

[1] 1.228369

    Но почему такие странные условия в первых строчках? Например:

k <- which(x>=4)[1]

Дело в том, что абсолютно точного значения  в нашем разбиении

x <- seq(xo, xn, by = dx)

может не оказаться (это зависит от кратности шага ). Поэтому мы ищем все номера, для которых  переваливает за 4. Функция which(x>=4) возвращает нам массив таких номеров, из которого мы выбираем первое значение which(x>=4)[1], соответствующее, либо самой точке , либо максимально близкой к ней.

 

Задания для самостоятельной работы

1. Решить приближенно с помощью схемы Эйлера задачу Коши:

на отрезке  и сравнить графически полученный результат с точным решением задачи:

.

2. Решить приближенно с помощью схемы Эйлера задачу Коши для уравнения Бернулли:

на отрезке  и сравнить графически полученный результат с точным решением задачи:

.

Указание: выбрать необходимый шаг .

3. Какую задачу Коши и на каком отрезке численно решает следующий код:

xo <- -2

xn <- 2

dx <- 0.1

yo <- sqrt(6)

x <- seq(from = xo, to = xn, by = dx)

n <- length(x); n

y <- rep(0,n)

y[1] <- yo

for (i in 1:(n-1)) {

y[i+1] <- y[i] + dx*(2*x[i]^3 - 3*x[i]+ 3/4)/y[i]

}

plot(x,y, type = " l ", lwd = 5, col = "blue")                          

abline(h = 0, v = 0, col = "gray40")                                     

lines(x, sqrt(x^4-3*x^2+1.5*x+5), type = " l ", lwd = 1, col = "yellow")

Каково точное решение данной задачи Коши? Корректен ли выбор шага в схеме Эйлера? Если некорректен – исправьте на достаточный.

4. * Определить приближенно точки локальных экстремумов и их значения для решения задачи Коши из предыдущего задания. Указание: проверяйте полученные ответы на соответствие графику.