Численное решение оду и систем оду: Численное решение системы дифференциальных уравнений (Лекция 14) Курс лекций по Информатике

Содержание

Численное решение системы дифференциальных уравнений (Лекция 14) Курс лекций по Информатике

Системой дифференциальных уравнений называется система вида

где x — независимый аргумент,

y_i — зависимая функция, ,

y_i|_x=x0 =y_i0 — начальные условия.

Функции y_i(x), при подстановке которой система уравнений обращается в тождество, называется решением системой дифференциальных уравнений.

Численные методы решения систем дифференциальных уравнений.

Метод Эйлера.
y_ij+1=y_ij+hf_i(x_i,y_1j y_2j..y_nj)

j — номер шага.

x_j+1=x_j+h
Модифицированный метод Эйлера.
k_i1=h*f_i(x_j,y_1j..y_nj)

k_i1=h*f_i(x_j+h,y_1j+k_i1.
.y_nj+k_i2)

y_ij+1=y_ij+(k_i1+k_i2)/2

x_j+1=x_j+h
Метод Рунге-Кутта четвертого порядка.
k_i1=hf_i(x_j,y_1j..y_nj)

k_i2=hf_i(x_j+h/2,y_2j+k_i1/2,..,y_nj+k_n1/2)

k_i3=hf_i(x_j+h/2,y_2j+k_i2/2,..,y_nj+k_n2/2)

k_i4=hf_i(x_j+h,y_1j+k_i2,..,y_nj+k_n3)

y_ij+1=y_ij+(k_i1+2k_i2+2k_i3+k_i4)/6

x_j+1=x_j+h

Дифференциальным уравнением второго порядка называется уравнение вида

F(x,y,у’,y»)=0

(1)

или

y»=f(x,y,y’).

(2)

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

Численно ищется частное решение уравнения (2), которое удовлетворяет заданным начальным условиям, то есть решается задача Коши.

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

(3)

Функция f₂(x, y₁, y) в систему (3) введена формально для того, чтобы методы, которые будут показаны ниже, могли быть использованы для решения произвольной системы дифференциальных уравнений первого порядка. Рассмотрим несколько численных методов решения системы (3). Расчетные зависимости для i+1 шага интегрирования имеют следующий вид. Для решения системы из n уравнений расчетные формулы приведены выше. Для решения системы из двух уравнений расчетные формулы удобно записать без двойных индексов в следующем виде:

Метод Эйлера.
у_1,i+1=у_1,i+hf₁(x_i, y_1,i, y_i),

у_i+1=у_i+hf₂(x_i, y_1,i, y_i),

x_i+1=x_i+h.
Метод Рунге-Кутта четвертого порядка.
у_1,i+1=у_1,i+(m₁+2m₂+2m₃+m₄)/6,

у_i+1=у_i+(k₁+2k₂+2k₃+k₄)/6,

m₁=hf₁(x_i, y_1,i, y_i),

k₁=hf₂(x_i, y_1,i, y_i),

m₂=hf₁(x_i+h/2, y_1,i+m₁/2, y_i+k₁/2),

k₂=hf₂(x_i+h/2, y_1,i+m₁/2, y_i+k₁/2),

m₃=hf₁(x_i+h/2, y_1,i+m₂/2, y_i+k₂
/2),

k₃=hf₂(x_i+h/2, y_1,i+m₂/2, y_i+k₂/2),

m₄=hf₁(x_i+h, y_1,i+m₃, y_i+k₃),

k₄=hf₂(x_i+h, y_1,i+m₃, y_i+k₃),

x_i+1=x_i+h,

где h — шаг интегрирования. Начальные условия при численном интегрировании учитываются на нулевом шаге: i=0, x=x₀, y₁=y₁₀, y=y₀.

Контрольное задание по зачетной работе.

Колебания с одной степенью свободы

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

Задание.

Численно и аналитически найти:

закон движения материальной точки на пружинке х(t),
закон изменения силы тока I(t) в колебательном контуре (RLC — цепи) для заданных в табл.1,2 режимов. Построить графики искомых функций.

Варианты заданий.

Таблица режимов

№	Режим
1	Свободные незатухающие колебания
2	Затухающее колебательное движение
3	Апериодическое движение
4	Предельное апериодическое движение
5	Вынужденное колебание без сопротивления
6	Вынужденное колебание без сопротивления, явление резонанса
7	Вынужденное колебание с линейным сопротивлением
8	Вынужденное колебание с линейным сопротивлением, явление резонанса

Таблица

Варианты заданий и номера режимов:

движение точки
RLC — цепь

Вар.	Задание	Вар.	Задание
1	а) 1,2,5	16	б)1,2,6
2	а) 1,3,6	17	а) 1,4,7
3	б)1,3,7	18	б)1,2,7
4	а) 1,4,8	19	а) 1,2,5
5	б)1,2,8	20	б)1,4,6
6	а) 1,4,7	21	а) 1,3,5
7	б)1,3,6	22	б)1,3,8
8	а) 1,4,5	23	б)1,4,5
9	б)1,3,8	24	а) 1,3,6
10	а) 1,3,5	25	б)1,4,7
11	б)1,4,6	26	а) 1,2,8
12	а) 1,2,7	27	б)1,4,8
13	б)1,2,5	28	а) 1,3,6
14	а) 1,2,6	29	б)1,3,7
15	б)1,4,7	30	а) 1,2,5

Рассмотрим более подробно порядок составления дифференциальных уравнений и приведения их к машинному виду для описания движения тела на пружинке и RLC-цепи.

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

Рис. Материальная точка на пружинке.

Уравнение движения (второй закон Ньютона) для материальной точки с учетом действия сил линейного сопротивления (-βx’), упругости пружины (-cx) и синусоидальной силы F₀∙sin(pt) может быть записано следующим образом

,

где m=1+int(n/2) — масса материальной точки, β — коэффициент сопротивления, с=2+int(n/3) — жесткость пружины, х — координата (х=0 в положении равновесия точки), t — время, p — частота вынужденных колебаний, F₀=n — амплитуда силы, n — номер варианта, int — целая часть числа. Параметры β, p и начальные условия выбираются самостоятельно с учетом рассматриваемого режима.
Колебательный контур (RLC цепь) показан на рис.

Уравнение падения напряжения в цепи переменного тока имеет вид

,

где L=1+int(n/2) — индуктивность, R — сопротивление, C=2+int(n/3) — емкость конденсатора, q — заряд, U₀=n — амплитуда напряжения, p — частота, — сила тока. Параметры R, p и начальные условия выбираются самостоятельно с учетом рассматриваемого режима.

Содержание отчета:

Название, цель работы и задание.
Математическое описание, алгоритм (структограмма) и текст программы.
Шесть графиков зависимости (три точные и три приближенные) x(t) или I(t), выводы по работе.

Сообщество Экспонента

вопрос
14.02.2023

Другое, Системы управления

Гидроцилиндр

Гидравлика

14. 02.2023

вопрос
12.02.2023

Системы управления, Электропривод и силовая электроника, Другое

Есть модель двигателя https://www.mathworks.com/help/sps/ref/bldc.html Мне необходимо построить такую же модель только из стандартных блоков. Mask -> Look under mask не работает. Как можно заглянут…

3 Ответа

Электропривод
BLDC

12.02.2023

вопрос
11.02.2023

Автоматизация испытаний

Как из MatLab опрашивать датчики, подключённые через переходник USB <-> I2C на основе микросхемы Ch441A ? Какие драйвера или библиотеки в ОС Windows необходимо установить для такой работы ? Датч…

Ch441A
USB
I2C

11.02.2023

вопрос
09.02.2023

Электропривод и силовая электроника

Здравствуйте, а существуют ли модели преобразователей частоты построенных по трехуровневой топологии с цепью заряда конденсаторов (фильтров) в цепи постоянного тока?

Публикация
07.02.2023

Больше ядер — больше возможностей! На предстоящем вебинаре мы расскажем о важной и актуальной теме: использование технологии многоядерных вычислений при моделировании энергосистем в режиме реального времени. При построении цифровых двойников энергосистем…

Приглашаем Вас на вебинар «Использование технологии многоядерных вычислений при моделировании энергосистем в режиме реального времени» 16 марта 2023 года.

Публикация
25.01.2023

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

Приглашаем вас на вебинар «Методы суррогатного моделирования сложных динамических систем», который пройдет 16 февраля в 10:00 по московскому времени.

MATLAB
Simulink
нейронные сети

25.01.2023

вопрос
18.01.2023

Есть входной аудиосигнал. Его надо пропустить через фильтр НЧ (600 Гц) в MATLAB. Как это сделать?

9 Ответов

Публикация
18.01.2023

Вебинар будет состоять из двух частей. В первой части будет обсуждаться роль цифровых двойников в предсказательном обслуживании. Далее будет построен цифровой двойник настоящего трансформатора малой мощности, используя MATLAB/Simulink, усилитель и КПМ РИТМ. Во…

Приглашаем на первый вебинар в этом году по теме: «Цифровой двойник трансформатора: на пути к интеллектуальному мониторингу» 9 февраля в 10:00.

MATLAB
Simulink
Машинное обучение
Predictive Maintenance
РИТМ

18.01.2023

вопрос
16.01.2023

Всем здравствуйте, стоит задача сделать генератор сигналов в Matlab, который формирует сигнал и выводит его через звуковую карту. Есть вот такой код Tm = 5;% Длина сигнала (с)Fd = 22050;% Частота диск…

5 Ответов

MATLAB
Обработка сигналов

16. 01.2023

Отвеченный вопрос
11.01.2023

Здравствуйте! Получил задание на разработку алгоритма и программы, реализующих оценку распределения модуля мгновенных значений фонограммы. 1) Разработать методику, алгоритм и программу оценки распреде…

8 Ответов

Численное решение систем обыкновенных дифференциальных уравнений (ОДУ 1-го порядка) | Сачин Чандрасекара

Решение систем ОДУ с помощью MATLAB

Фото Джонатана Кемпера на Unsplash

Обыкновенные дифференциальные уравнения, или ОДУ, играют важную роль в различных приложениях. И часто для их решения лучше всего подходят численные методы. Давайте посмотрим, как мы делаем это в MATLAB, решая системы обыкновенных дифференциальных уравнений первого порядка. Здесь мы обеспокоены применением распространения ВИЧ-инфекции.

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

Три дифференциальных уравнения о распространении ВИЧ

Модель содержит шесть параметров kr1, kr2, kr3, kr4, kr5 и kr6, , которые обеспечивают скорость гибели клеток, распространение инфекции, репликацию вируса и другие процессы, влияющие на распространение ВИЧ в организме.

Как упоминалось выше, имеется шесть параметров с соответствующими значениями.

Значения параметров

Я надеюсь, вы все знаете, что мы можем получить решение одного дифференциального уравнения с помощью MATLAB. Затем давайте посмотрим, как мы делаем это с помощью функции ode45 в MATLAB, решая систему ОДУ. Как указывалось ранее, функция ode45 является наиболее часто используемым решателем ОДУ в MATLAB.

Для этого сначала нужно преобразовать систему ОДУ в новый формат. Следовательно, первым шагом является объединение всех этих ОДУ в одно уравнение для одной переменной, которую мы назовем X.

Значения в терминах компонентов вектора «X»

Стратегия состоит в том, чтобы создать X , вектор-столбец с тремя компонентами, каждая из которых соответствует одной из переменных H, I, и V . Затем мы берем производную каждого элемента X . Соответствующее уравнение определяет каждое значение исходной системы.

Теперь, когда у нас есть векторизованная система ОДУ, следующим шагом будет создание функции MATLAB. Эта функция вычисляет значение вектора производной для использования с функцией ode45. Наша функция « модель» будет генерировать один выход dXdt для заданных входных значений t и X.

Функция, построенная с соответствующими значениями параметров

Последний шаг — вызвать функцию ode45, передав дескриптор нашей производной функции, вектор, содержащий начальное и конечное значения времени. Здесь начальные условия: H_0 = 1000000, I_0 = 0 и V_0 = 100.

Заданные начальные условия

Решение нашего дифференциального уравнения записывается в выходные переменные0010 и ‘tsol’ , когда мы запускаем этот скрипт. Первым выходом является вектор, содержащий значения времени, используемые для вычисления численного решения. Второй вывод представляет собой матрицу, содержащую рассчитанное значение переменной X для каждого времени в ‘tsol’ . Затем отдельные столбцы могут быть извлечены для дополнительного анализа или визуализации.

Новак М. и Мэй Р. М. Вирусная динамика: математические принципы иммунологии и вирусологии: математические принципы иммунологии и вирусологии . Oxford University Press, 2000.

Обыкновенные дифференциальные уравнения в R

Введение

Часто бывает естественным сформулировать ожидаемые отношения между переменными в терминах дифференциальных уравнений (ДУ). Проще говоря, такие уравнения выражают взаимосвязь между значениями переменных и скоростью изменения этих значений. Это хорошо согласуется с тем, как мы обычно думаем о системах. Решением множества ДУ является траектория переменных во времени. Таким образом, ДЭ часто естественным образом возникают при анализе динамических данных.

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

Начнем с простой модели DE. Далее мы узнаем, как получать численные решения в Р . Наконец, мы встраиваем численное интегрирование в контекст нелинейного метода наименьших квадратов, чтобы получить оценки параметров и выполнить вывод.

Модель SIR

Классическая компартментальная модель SIR делит популяцию хозяев на три класса: восприимчивые, инфицированные, выздоровевшие (см. диаграмму ниже). Модель описывает, как количество людей в каждом из этих классов меняется со временем. Рождения моделируются как потоки «из ниоткуда» в восприимчивый класс; смерти моделируются как потоки из отсека S, I или R в «никуда».

Схема отсековой модели SIR.

Популяция хозяев делится на три класса в зависимости от их инфекционного статуса:

S, восприимчивые хозяева;
I, инфицированные (и заразные) носители;
Р, хозяева удалены из заразной популяции.

распространенность

Решение ОДУ в R

Поскольку эти уравнения нелинейны, неудивительно, что их нельзя решить аналитически. Однако мы можем вычислить траекторий модели с непрерывным временем, такой как эта, путем численного интегрирования уравнений. Точное выполнение этого требует большого количества вычислений, и есть умные и не очень умные способы сделать это. Эта очень распространенная проблема очень тщательно изучалась численными аналитиками на протяжении поколений. Результатом всей этой работы является то, что, когда уравнения представляют собой гладкие функции с хорошим поведением, превосходные алгоритмы численного интегрирования легко доступны для вычисления приближенных решений с высокой точностью. В частности, R имеет несколько сложных решателей дифференциальных уравнений, которые (для многих задач) дают очень точные решения. Эти алгоритмы гибки, автоматически выполняют проверки, выдают информативные ошибки и предупреждения. Чтобы использовать пакет решателя дифференциальных уравнений, мы загружаем deSolve package

 library (deSolve)

. Решатель ODE, который мы будем использовать, называется ode . Взгляните на страницу помощи для этой функции.

Страница справки сообщает нам, что вызов ode выглядит как

 ode(y, times, func, parms,
    method = c("lsoda", "lsode", "lsodes", "lsodar", "vode", "daspk",
    "эйлер", "rk4", "ode23", "ode45", "radau",
    "bdf", "bdf_d", "adams", "impAdams", "impAdams_d", "итерация"),
    .  {}}} & = \ gamma \, I \\
\end{выровнено}
\] Чтобы закодировать эти уравнения в форме, подходящей для использования в качестве  func  аргумент для  ode , нам нужно написать функцию. Например: 
  closed.sir.model <- function (t, x, params) {
  ## сначала извлеките переменные состояния
  S <- х[1]
  я <- х[2]
  R <- х[3]
  ## теперь извлеките параметры
  бета <- параметры["бета"]
  гамма <- параметры["гамма"]
  Н <- С+И+Р
  ## теперь закодируйте уравнения модели
  dSdt <- -бета*S*I/N
  dIdt <- бета*S*I/N-гамма*I
  dRdt <- гамма*I
  ## объединить результаты в один вектор
  dxdt <- c(dSdt,dIdt,dRdt)
  ## вернуть результат в виде списка!
  список (dxdt)
} 
 Обратите внимание, что порядок и тип аргументов и выходных данных этой функции должны точно соответствовать ожиданиям  или . Так, например, переменная времени  t  должна быть первым аргументом, даже если, как здесь, в функции от времени ничего не зависит. [Когда правая часть ОДУ не зависит от времени, мы говорим, что ОДУ являются  автономными  ; когда они явно зависят от времени, говорят, что они являются  неавтономными  . ] Заметим также, что  ода  ожидает, что значения ODE RHS будут первым элементом списка   .
 Теперь мы можем вызывать  или  для вычисления траекторий модели при условии, что у нас есть значения параметров и  начальных условий  , т. е. значения переменных состояния \(S\), \(I\), и \(R\) в некоторый начальный момент времени. Если мы думаем о заболевании, например кори, и измеряем время в годах, мы можем использовать что-то вроде:
 пармс <- c(бета=400,гамма=365/13) 
 Теперь укажем моменты времени, когда нам нужны решения, и зададим начальные условия:
 раз <- seq(from=0,to=60/365,by=1/365/4)
xstart <- c(S=999,I=1,R=0) 
 Затем мы вычисляем траекторию модели с помощью команды  ode  и сохраняем результат в кадре данных:
 библиотека (tidyverse)
ода(
  func=closed.sir.model,
  у=хстарт,
  раз = раз,
  пармс=пармс
) %>%
  as.data.frame() -> out 
 и постройте результаты с помощью команд:
 выход %>%
  собрать (переменная, значение, время) %>%
  ggplot(aes(x=время,y=значение,цвет=переменная))+
  geom_line (размер = 2) +
  тема_классика()+
  labs(x='время (год)',y='количество людей') 
 Упражнение 
 Предположим, вы предпочитаете измерять время в днях. Измените параметры соответствующим образом и проверьте ваши изменения.
 Давайте изучим, как кривая эпидемии зависит от скорости передачи \(\beta\) и инфекционного периода. В частности, мы исследуем, как меняется кривая эпидемии при изменении \(\beta\) от 20 до 500 и инфекционного периода от 5 до 30 дней.
 бета-значения <- c(20,50,500)
ipvals <- c(5,10,30)/365
gammavals <- 1/ipvals 
 expand.grid(beta=betavals,gamma=gammavals)%>%
  group_by(бета,гамма) %>%
  делать(
    {
      ode(func=closed.sir.model,y=xstart,times=times,
        parms=c(beta=.$beta,gamma=.$gamma)) %>%
        as.data.frame()
    }
  ) %>%
  ggplot(aes(x=время,y=I))+
  геометрическая_линия()+
  facet_grid(beta~gamma,scales='free_y',labeller=label_both)+
  theme_bw() 
 Динамика SIR в открытой популяции 9{}}} &= \gamma\,I-\mu\,R\\
\end{выровнено}
\] 
 Мы должны соответствующим образом изменить функцию ОДУ:
 open.sir.model <- function (t, x, params) {
  B <- параметры["B"]
  бета <- параметры["бета"]
  мю <- параметры["му"]
  гамма <- параметры["гамма"]
  N <- х[1]+х[2]+х[3]
  dSdt <- B - бета*x[1]*x[2]/N - mu*x[1]
  dIdt <- бета*x[1]*x[2]/N - (мю+гамма)*x[2]
  dRdt <- gamma*x[2] - mu*x[3]
  список (с (dSdt, dIdt, dRdt))
} 
 Нам нужно указать уровень рождаемости/смертности в дополнение к двум параметрам, которые мы указали ранее:
 parms <- c(B=20,mu=1/50,beta=400,gamma=365/13) 
 Интегрируем уравнения, как и раньше:
 ode(
  func=open. sir.model,
  у=хстарт,
  раз = последовательность (от = 0, до = 25, по = 1/365),
  пармс=пармс
) %>%
  as.data.frame() -> out 
 Мы можем построить график зависимости каждой из переменных состояния от времени и \(I\) от \(S\), используя команды:
 out %>%
  собрать (переменная, значение, время) %>%
  ggplot(aes(x=время,y=значение,цвет=переменная))+geom_line()+
  facet_wrap(~переменная,масштабы="free_y",ncol=1)+
  масштаб_y_log10()+
  направляющие(цвет=ЛОЖЬ)+theme_bw() 
 выход %>%
  ggplot(aes(x=S,y=I))+geom_path()+
  scale_x_log10()+scale_y_log10()+theme_bw() 
 Упражнение 
 Исследуйте динамику системы для различных значений параметров \(\beta\) и \(\gamma\) путем моделирования и построения траекторий в виде временных рядов и в фазовом пространстве (например, \(I\) против \(S\)). Используйте те же значения \(B\) и \(\mu\), которые мы рассмотрели выше. Как значение \(R_0=\tfrac{\beta}{\gamma+\mu}\) влияет на результаты? Возможно, вам будет полезно обернуть вышеуказанные вызовы в  ода  в  манипулировать  .
 Beyond SIR 
 Стратегию сегментарного моделирования можно использовать для моделирования огромного количества инфекций. Этот раздел и более сложные упражнения в нем делают первые шаги в этом направлении.
 Латентность: модель SEIR 
 Люди, инфицированные большинством патогенов, часто не сразу становятся заметно заразными. Большинство инфекций имеют латентный период, в течение которого передача маловероятна или невозможна. Мы можем уловить эту сложность, введя латентную фазу в нашу модель SIR. На приведенной ниже диаграмме показаны отсеки и потоки модели SEIR. 9{}}} &= \gamma\,I.\\
\end{выровнено}
\] Следующее кодирует правую часть этих уравнений в форме, подходящей для использования с  deSolve::ode  .
 closed.seir.model <- function (t, x, params) {
  infec <- params["beta"]*x[1]*x[3]/sum(x)
  symp <- params["сигма"]*x[2]
  recov <- параметры["gamma"]*x[3]
  список (c (-infec, infec-symp, symp-recov, recov))
} 
 Упражнение 
 Адаптируйте приведенные выше коды модели, чтобы включить рождаемость и смертность, как в открытой модели эпидемии SIR.
No related posts.