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

Министерство науки и высшего образования Российской Федерации 

НОВОСИБИРСКИЙ ГОСУДАРСТВЕННЫЙ ТЕХНИЧЕСКИЙ УНИВЕРСИТЕТ 
__________________________________________________________________________ 
 
 
 
 
Д.В. ВАГИН 
 
 
 
ЧИСЛЕННОЕ МОДЕЛИРОВАНИЕ  
ДИНАМИЧЕСКИХ СИСТЕМ,  
ОПИСЫВАЕМЫХ  
ОБЫКНОВЕННЫМИ  
ДИФФЕРЕНЦИАЛЬНЫМИ  
УРАВНЕНИЯМИ 
 
 
Утверждено Редакционно-издательским советом университета 
в качестве учебного пособия 
 
 
 
 
 
 
 
 
 
 
 
НОВОСИБИРСК 
2019 

УДК 519.62(075.8) 
В 124 
 
Рецензенты: 
канд. техн. наук, ведущий инженер-программист  
АО «Завод Экран» А.В. Волкова 
д-р техн. наук, профессор М.Э. Рояк 
 
 
 
Работа подготовлена на кафедре прикладной математики НГТУ 
 
 
 
Вагин Д.В. 
В 124  
Численное моделирование динамических систем, описываемых обыкновенными дифференциальными уравнениями: учебное пособие / Д.В. Вагин. – Новосибирск: Изд-во НГТУ, 2019. – 
63 с. 

ISBN 978-5-7782-3941-8 

В данном учебном пособии рассмотрены элементы теории из раздела численных методов решения систем обыкновенных дифференциальных уравнений. Пособие может быть рекомендовано как для самостоятельного изучения курса «Численное моделирование динамических систем, описываемых обыкновенными дифференциальными 
уравнениями», так и для подготовки к выполнению практических заданий. 
 
 
 
 
 
 
 
УДК 519.62(075.8) 
 
ISBN 978-5-7782-3941-8 
© Вагин Д.В., 2019 
 
© Новосибирский государственный 
 
технический университет, 2019 

 
 
 
 
 
 
 
ВВЕДЕНИЕ 
 
Для изучения поведения (эволюции) динамических систем необходимо решать описывающие их системы дифференциальных уравнений. 
Сложные динамические системы описываются большими системами 
нелинейных дифференциальных уравнений разного порядка. Для 
большинства практических задач эти системы могут быть сведены к 
системам обыкновенных дифференциальных уравнений. Для решения 
подобных систем используются численные методы. 
В данном учебном пособии рассматриваются основные элементы 
теории численного решения систем обыкновенных дифференциальных 
уравнений. Среди них: 
 Явные и неявные схемы; 
 Аппроксимация, сходимость и устойчивость численных схем; 
 Одношаговые и многошаговые схемы; 
 Методы Эйлера, Рунге–Кутты и Адамса; 
 Методы типа «прогноз-коррекция»; 
 Системы обыкновенных дифференциальных уравнений; 
 Жесткие задачи. 
Помимо этого рассматриваются основные принципы написания 
программ для имитационного моделирования сложных технических 
систем. 
 
 

1. ДИСКРЕТНАЯ ЗАДАЧА КОШИ 
 
Рассмотрим задачу Коши: 

 
( )
( , ( ))
y t
f t y t


, 
(1.1) 

 
0
0
( )
y t
y

,    
0
[ , ]
t
t
T

. 
(1.2) 

Ее решением является функция 
( )
y t . Процесс решения уравнения (1.1) с начальным условием (1.2) (нахождение функции 
( ))
y t
 также называют интегрированием уравнения (1.1). Когда аналитическое 
решение (1.1) невозможно или чрезвычайно сложно, прибегают к численному интегрированию. 
Для использования численного метода необходимо перейти от непрерывной задачи, решением которой является функция 
( )
y t , определенная для всех 
0
[ , ]
t
t
T

, к дискретной задаче, в которой функция 
( )
y t  заменена набором значений  


0
1
2
, , 
,..., 
i
N
y
y
y
y
y

, ассоциированных со значениями  


0
1
2
...
i
N
t
t
t
t
t
T






. То есть если решения непрерывной и дискретной задачи совпадают, то будет выполнено: 
1
1
( )
y t
y

, 
2
2
( )
y t
y

, ..., 
(
)
N
N
y t
y

. Отметим, что 
0
0
( )
y t
y

 
выполняется всегда, так как является исходным данным для обеих задач. Набор значений  
it
 называется сеткой, а значения it  – узлами 
сетки. Величины 
1
n
n
n
h
t
t 


, 
1, 2,..., 
,
n
N

 называются шагами 
сетки. Набор значений  
iy
 называется дискретным аналогом функции ( )
y t  (или сеточной функцией). Если 
0
(
) /
,
n
N
h
h
t
t
N



 для всех 
1, 2,..., 
,
n
N

 то сетка  
it
 называется равномерной и обозначается 

как ht , а сеточная функция 
iy
 – как 
h
y  [1]. 

Для того чтобы сформулировать дискретную задачу, необходимо 
левую и правую части уравнения (1.1) заменить на их дискретные аналоги. Для упрощения выкладок ограничимся рассмотрением равномерных сеток. Отметим, что данное упрощение не нарушает общности 
излагаемой теории. 
Аппроксимируем 
( )
y t

, 
1
[ , 
]
n
n
t
t
t 

 с помощью некоторой разностной схемы, включающей в себя 
1
k   значение сеточной функции 
h
y : 

 
0
1
1
1
1
( )
(
...
)
n
n
k
n k
y t
y
y
y
h

 



 

 
. 
(1.3) 

Полагается, что 
0
0
 
. При этом аппроксимация функции 
( , ( ))
f t y t
 в общем виде может быть записана, как 

 
1
1
( , ( ))
( , 
,..., 
, 
, )
n
n k
n
n
f t y t
t
y
y
y
h
 

 
. 
(1.4) 

Таким образом, дискретным аналогом уравнения (1) является уравнение 

 
1
1
1
0

1
(
, 
,..., 
, 
, )
k

j n
j
n
n k
n
n
j
y
t
y
y
y
h
h
 
 



 

. 
(1.5) 

Из уравнения (1.5) можно найти значение 
1
n
y  , если известны k  

предыдущих значений 
1
n k
y   , ..., 
n
y  сеточной функции 
h
y . Поэтому 
основанные на решении уравнения (1.5) методы получили название  
k-шаговых. Очевидно, что для того чтобы использовать k-шаговый 
метод, необходимо задать k  начальных значений: 

 
0
0
1
1
1
1
(
)
,  
( )
, ...,  
(
)
h
h
h
k
k
y
t
y
y
t
y
y
t
y





. 
(1.6) 

Уравнение (1.5) с начальными условиями (1.6) будем называть дискретной задачей Коши. Так как ее решением является сеточная функция 
h
y , соответствующие методы нахождения 
h
y  также называются 
сеточными методами. Принято считать, что уравнение (1.5) задает 
численный метод, который также называется численной схемой. Далее 

будем отождествлять свойства уравнения (1.5), численного метода 
(схемы) и соответствующей дискретной задачи Коши. 
Отметим также, что если в правой части (1.5) функция   не зависит от 
1
n
y  , то соответствующая численная схема называется явной. 
В противном случае (функция   зависит от 
1
n
y  ) соответствующая 
численная схема называется неявной. Очевидно, что для нахождения 
1
n
y   с помощью явной схемы достаточно подставить в (1.5) известные из начальных условий или уже вычисленные на предыдущих 
шагах интегрирования значения искомой сеточной функции 
h
y . Поэтому явные схемы очень просты с точки зрения их программной реализации. Для нахождения 
1
n
y   в (1.5) с помощью неявной схемы в 
общем случае приходится решать нелинейное уравнение. Решение 
нелинейного уравнения, даже таким быстро сходящимися методом 
как, например, метод Ньютона, требует существенных вычислительных затрат относительно простого вычисления правой части в явном 
методе [2]. Тем не менее положительные свойства неявных схем, как 
правило, позволяют использовать более грубый шаг интегрирования, 
что может частично или полностью компенсировать временные затраты на решение нелинейного уравнения. К тому же в случае, если 
решаемая задача является жесткой, то использование явной схемы 
зачастую вообще не позволяет получить сколько-нибудь приемлемое 
решение.  

АППРОКСИМАЦИЯ 

Итак, мы сформулировали дискретную задачу Коши, но сразу возникает вопрос о том, как сильно она отличается от исходной непрерывной задачи. 
Заменим в (1.5) значения 
1
n j
y    сеточной функции 
h
y  соответ
ствующими значениями 
1
(
)
n j
y t  
 функции 
( )
y t . После этой подстановки левая и правая части дискретного уравнения могут начать отличаться. Рассмотрим их разность: 

 


1
1
1
0

1
(
)
, (
),..., (
), (
), 
k

n
j
n
j
n
n k
n
n
j
y t
t
y t
y t
y t
h
h
 
 





 

. (1.7) 

Сеточная функция 
n
  называется погрешностью аппроксимации 
дискретного уравнения (1.5). Говорят, что дискретное уравнение (1.5) 
аппроксимирует дифференциальное уравнение (1.1), если max
0
n


 
при 
0
h 
, 
1
k
n
N
 

. Также логично предположить, что в этом 
случае 

 
1
0

1
(
)
( )
k

j
n
j
j
y t
y t
h
 






; 
(1.8) 

 
1
1
( , (
),..., ( ), (
), )
( , ( ))
n
n k
n
n
t
y t
y t
y t
h
f t y t
 



. 
(1.9) 

Из (1.8) в частности следует, что 

 
0
0
k

j
j

 

. 
(1.10) 

Качество аппроксимации непрерывного уравнения дискретным 
оценивают с помощью порядка аппроксимации, выражающего порядок малости погрешности аппроксимации относительно шага h. Говорят, что дискретное уравнение (1.5) аппроксимирует дифференциальное уравнение (1.1) с p-м порядком, если 

 
max
p
n
Ch


,   
0
p 
,   
1
k
n
N
 

. 
(1.11) 

Помимо погрешности аппроксимации для оценки качества одношаговых методов также используют другую величину – локальную 
погрешность. Пусть 
1
n
y   – значение, найденное из уравнения 

 
1
1
(
)
(
, (
), 
, )
n
n
n
n
n
y
y t
t
y t
y
h
h




 
. 
(1.12) 

Уравнение (1.12) – это уравнение (1.5) для одношагового метода,  
в которое вместо 
n
y  подставлено точное значение 
( )
n
y t
. Тогда разность 
1
1
(
)
n
n
n
l
y t
y




 называется локальной погрешностью (погрешностью на шаге) метода. Нетрудно показать, что для явных схем локальная погрешность и погрешность аппроксимации связаны простым 
соотношением n
n
l
h
 
. Хотя для неявных схем локальная погреш

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

УСТОЙЧИВОСТЬ 

Порядок аппроксимации является важной характеристикой дискретной задачи, однако еще более важно, чтобы решение дискретной 
задачи было устойчивым к малым возмущениям начальных значений и 
правой части уравнения. Отсутствие этого свойства делает соответствующую численную схему непригодной для практического использования. 
Внесем в правую часть уравнения (1.5) и начальные условия (1.6) 
произвольные малые возмущения 
n
  и 
0
 , 
1 , 
2
 ,..., 
1
k

 соответ
ственно. Положим 
*
0
0
0
y
y

  , 
*
1
1
1
y
y

  , ..., 
*
1
1
1
k
k
k
y
y




 
. 

Обозначим за 
*h
y
 решение соответствующей возмущенной задачи 

 


*
*
*
*
1
1
1
0

1
, 
,..., 
, 
, 
k

j n
j
n
n k
n
n
n
j
y
t
y
y
y
h
h
 
 



 
 

; 
(1.13) 

 
*
*
*
*
*
*
0
0
1
1
1
1
(
)
,  
( )
, ..., 
(
)
h
h
h
k
k
y
t
y
y
t
y
y
t
y





. 
(1.14) 

Будем называть дискретную задачу Коши и соответствующий численный метод устойчивыми на конечном отрезке (или просто устойчивыми), если при всех 
0
h
h

 (где 
0
h  достаточно мало) справедливо 
неравенство 

 

1
*

0
0
1
1
max
( )
max
N

n
n
n
n
n N
n k
n k
y
y
K T
h


 
  
 













, 
(1.15) 

где величина 
 
K T  не зависит от h, n
  и 
n
 , входящих в (1.15). Выражение (1.15) по сути означает, что ошибка решения, вызванная возмущением начальных значений и правой части уравнения, может быть 
сделана сколь угодно малой, если соответствующим образом уменьшить породившее ее возмущение. Последнее предложение и есть 
наиболее часто используемое определение понятия устойчивой задачи. 

НУЛЬ-УСТОЙЧИВОСТЬ 

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

 
( )
0
y t


. 
(1.16) 

Для всякого практически пригодного метода (решение непрерывным образом зависит от начальных значений) ошибку 
*
n
y  можно оценить как 

 
*

0
0
1
max
( ) max
h
n
n
n
n N
n k
y
y
K
T
 
  




. 
(1.17) 

Согласно (1.17) ошибка, внесенная в начальные значения, возрастает в 
( )
h
K
T  раз. То есть величина 
( )
h
K
T  играет роль числа обусловленно
сти метода. Величина 
( )
h
K
T  зависит от h  и от T , и ее рост с ро
стом T  в общем случае допустим. Однако если величина 
( )
h
K
T  может неограниченно возрастать с уменьшением h, то уменьшение шага 
интегрирования приведет не к уточнению, а наоборот, к неограниченному росту погрешности. Соответственно для устойчивости метода 
естественно потребовать, чтобы для достаточно малых h K  не зависела от h, т. е. чтобы выполнялось 

 
*

0
0
1
max
( ) max
n
n
n
n N
n k
y
y
K T
 
  




. 
(1.18) 

В (1.16) 
0
f 
. Следовательно, согласно (1.5) (при фиксированном 
h и 
0
 
), получаем 

 
1
0
0
k

j n j
j
y  





. 
(1.19) 

Уравнение (1.19) относится к классу линейных однородных разностных уравнений k-го порядка с постоянными коэффициентами. 
Если 
*h
y
 – решение этого уравнения, соответствующее возмущенным 

начальным значениям 
*
0
y , 
*
1y ,..., 
*
1
k
y  , то в силу линейности уравне
ния (1.19) погрешность 
*
n
n
n
y
y



 также будет решением (1.19): 

 
1
0
0
k

j n
j
j
 


 


. 
(1.20) 

Будем искать частное решение (1.20) в виде 
n
n
q


. Подставив его 

в (1.20) и разделив на общий множитель 
1
n k
q   , получим уравнение 

 
0
0
k
k
j
j
j
q 






, 
(1.21) 

которое называют характеристическим уравнением, соответствующим 
методу (1.5). Левая часть (1.21) называется характеристическим многочленом. 
Из теории линейных разностных уравнений можно показать, что 
уравнение (1.18) для задачи (1.16) будет выполнено (метод будет обладать нуль-устойчивостью), если корни характеристического уравнения лежат внутри или на границе единичного круга комплексной плоскости, причем среди корней, лежащих на границе единичного круга, 
нет кратных. Сформулированное условие также называют корневым 
условием. 

АБСОЛЮТНАЯ УСТОЙЧИВОСТЬ 

Нуль-устойчивость гарантирует устойчивое развитие погрешности 
при 
0
h 
 только в том случае, когда отрезок интегрирования 

0, 
t
T  
фиксирован. При этом наличие нуль-устойчивости вовсе не исключает 
того, что сколь угодно малая погрешность в начальных значениях при 
T   может приводить к сколь угодно большой погрешности решения. Часть методов, непригодных для практического применения при 
решении задачи Коши на больших временных отрезках, можно отбросить, исследовав результат их применения к решению модельной задачи  

 
y
y
   ,    
0
0
( )
y t
y

. 
(1.22)