Жесткие дифференциальные уравнения (интегрирование в среде Mathcad)
готовящейся к публикации в издательстве МЭИ
The stiff
Differential equations (integration in Mathcad)
Table of contents:
Introduction
Method
rkfixed. Instability of the solution
Method
rkadapt
Method
stiffr. The solution of a stiff model problem.
Method
stiffr. The solution of the chemical kinetics equations
Explicit
and implicit methods
Literature
The
elementary model of control system for some object Y is considered which provides
change Y according to the law F (X): Y = F(X). By virtue of the
different reasons (for example, random exterior actions or special initial
state of object) there can be considerable diversions from the ordered law
Y - F(X) ¹ 0.
The
control should ensure rapid returning to an equilibrium trajectory F (X) (see
figure).
There
are many examples of the such systems:
q
The robot pilot provides the given small height of flight h above a
relief [F(X)‑h], where X is the route on a pilot chart. The diversions
from given height should be promptly eliminated, catastrophe otherwise is
possible. It is necessary to ensure rapid returning of Y to the ordered value F
(X), hence parameter a should be great.
q
Temperature of object Y should hold by means of thermostat at the
necessary level F (X) where X is a time. The intensity a of heat exchange
between object and fluid in thermostat should be great, that the temperature
balance between object and thermostat, Y = F (X), was promptly
renewed.
q
q
Etc.
The
behavior of such systems by means of numerical integration by various methods
is explored:
ü
by explicit Runge-Kutta
method rkfixed with a fixed step,
ü
by explicit Runge-Kutta
method rkadapt with step adaptation,
ü
by special implicit method stiffr for the stiff differential
equations.
During
numerical experiments the characteristics for stiff systems are found out. It is
shown, that the special method stiffr has the important advantage before
others. Two basic schemes of a numerical integration - explicit and implicit -
is considered and is shown, that the stiff equations should be solved by
implicit methods, such as stiffr.
The
stiff equations frequently meet at mathematical model operation of technical
systems. However, despite of practical importance of this question, it is
difficult for the student or engineer to find the accessible description of
this phenomenon. The publication fills somewhat this blank. On substantial
examples, by means of Mathcad, the effective expedients of the solution of such problems
are demonstrated.
The
materials of the publication could be useful to Mathcad users, because its help
system almost nothing contains on a discussed problem.
____________________________________________________________________________________________________________
Содержание
Жесткие дифференциальные уравнения
(интегрирование в системе Mathcad)
Метод rkfixed. Неустойчивость решения
Mетод
stiffr. Решение модельной жесткой задачи.
Метод stiffr.
Решение уравнений химической кинетики.
Рассмотрим простейшую модель управления некоторым объектом Y. Пусть система управления должна обеспечить изменение Y согласно закону F(X):
(1)
Предположим, что в силу разных причин (например, случайных внешних воздействий или особого начального состояния объекта) возможно возникновение значительных отклонений от предписанного закона
Необходимо обеспечить быстрое возвращение на равновесную траекторию F(X). Поэтому скорость возвращения, т.е. производная dY/dX, должна быть пропорциональна отклонению и противоположна ему по знаку, а коэффициент пропорциональности должен быть большой величиной. Такая модель представляется дифференциальным уравнением (ДУ) с большим параметром a:
Поясним смысл параметра a,
оценив по порядку величины период DX
быстрого изменения функции в окрестности некоторой точки X0.
Для этого просто заменим в ДУ (2) производную
ее конечно-разностной аппроксимацией, а величины в правой части возьмем при X=X0
(см. Рис. 1 и выражение (3)).
Рис. 1. Быстрое уменьшение отклонения
при большом значении параметра a
При операциях оценки примем во внимание, что функция управления F(X) есть медленно меняющаяся на больших промежутках X функция. В нижеследующей записи выражения на Mathcad’е нам пришлось воспользоваться логическим знаком равенства, чтобы сформировать уравнение, хотя по смыслу приближенной операции оценки более подошел бы знак ~:
Итак, если параметр a сделать большой величиной, то ликвидация отклонений будет происходить очень быстро, на малых промежутках DX (Рис. 1), протяженность которых по порядку величины равна обратному значению большого параметра a .
Существует много примеров систем, описываемых дифференциальными уравнениями вида (2):
q Автопилот обеспечивает заданную малую высоту h полета над рельефом [F(X)‑h], где X – маршрут полета на карте. Возникшие по той или иной причине отклонения высоты должны быть быстро устранены, иначе возможна катастрофа. Необходимо обеспечить быстрое возвращение величины Y к предписанному значению F(X), следовательно параметр a должен быть большим.
q Температура объекта Y устанавливается на нужном уровне F(X) в зависимости от времени X с помощью термостата. Интенсивность теплообмена a между объектом и жидкостью в термостате должна быть большой, чтобы быстро устанавливалось тепловое равновесие между объектом и термостатом, Y=F(X).
q Концентрация Y компонента смеси в химическом реакторе быстро стремится к равновесному значению F вследствие высокой скорости реакции.
и т.п.
Мы исследуем далее поведение таких систем посредством численного интегрирования различными методами, а именно:
ü явным методом Рунге-Кутта с фиксированным шагом rkfixed,
ü явным методом Рунге-Кутта с адаптацией шага rkadapt,
ü специальным (неявным) методом stiffr для жестких дифференциальных уравнений.
В процессе исследования мы обнаружим характерные признаки, определяющие собственно свойство жесткости. Мы увидим, что специальный метод stiffr имеет важное преимущество перед двумя другими. Наконец, на простейшем примере мы рассмотрим две основные численные схемы - явную и неявную - и покажем, что жесткие уравнения следует решать неявными методами (stiffr - такой метод).
В качестве закона управления выберем медленно меняющуюся периодическую функцию с волновым числом k, k »1:
(4)
Зададим основной параметр задачи a . Это должна быть большая величина, чтобы
обеспечить быстрое возвращение величины Y к предписанному (равновесному) значению F(X).
Запишем правую часть модельного дифференциального уравнения (2):
(5)
Получим несколько вариантов решения при различных значениях шага интегрирования h и при различных начальных условиях. Смысл параметров встроенной функции rkfixed понятен из обозначений, которые мы старались сделать мнемоническими.
Рассчитаем эволюцию системы из начального состояния, заметно отличающегося от равновесного:
Область интегрирования зададим протяженной по сравнению с периодом быстрого изменения (Xend>>1/a):
Число шагов NofSteps и соответственно шаг интегрирования h выберем как
Это относительно малый шаг интегрирования. Если попытаться увеличить шаг, то возникнут неприятности (см. пункт «Вариант 2»).
Обратимся к встроенной функции интегрирования с фиксированным шагом rkfixed:
Полученная функция Y(X) представлена на Рис. 2. Обсуждение результатов будет проведено ниже.
Пусть в промежуточной точке Xbeg > 0
возникает значительное отклонение:
Шаг интегрирования сохраним прежний:
Рассчитаем эволюцию отклонения, возникшего в промежуточной точке:
Результаты численного интегрирования для двух рассчитанных выше вариантов представлены далее на Рис. 2 как функции Y(X) и Y1(X1). Там же показано равновесное решение F(X).
Рис. 2. Эволюция отклонений,
возникающих в различных точках X.
Как видно из графиков, начальное отклонение, заданное в точке X=0, быстро затухает благодаря большому значению a. Далее изменение Y происходит медленно в соответствии с заданной функцией управления F(X).
Такое поведение решения сохраняется для всей области интегрирования. Отклонения, возникающие в любой точке X, эволюционируют аналогичным образом, что иллюстрируется ветвью решения Y1.
Итак, решение модельной задачи характеризуется сменой периодов быстрого и медленного изменения искомой величины. Отклонения, возникшие в любой точке области интегрирования, порождают быстроменяющиеся решения.
Эти признаки составляют определение жесткой задачи.
Мы продолжим далее исследование свойств модельного жесткого уравнения (2) и рассмотрим возможности регулирования шага интегрирования.
В связи с характерным поведением искомой функции - сменой периодов быстрого и медленного изменения - появляется идея интегрировать с разным шагом: малым в области быстрого изменения и большим - в области медленных изменений (с очевидной целью - сэкономить время и память).
Расчеты в следующем пункте и рис. 3 показывают, что при этом произойдет.
Попытаемся увеличить шаг вдвое в области слабых изменений, начиная с точки i_beg (см. Рис. 2):
Стартуем со значения Y(Xbeg), используя решение, полученное выше в пункте «Исходный вариант решения». Это решение сохранено в виде массива Y, и нужно просто выбрать элемент с нужным индексом i_beg :
Увеличим шаг вдвое, разбив оставшийся интервал (Xend -Xbeg) на вдвое меньшее число интервалов, чем в рассмотренных выше примерах:
Итак, до точки i_beg включительно решение получено в пункте при счете с малым шагом, а далее мы пытаемся продолжить решение с увеличенным вдвое шагом:
Результат показан на Рис. 3.
Рис. 3. Неустойчивость решения при
интегрировании жесткой задачи методом rkfixed c увеличенным шагом.
Кривая Y (это по-прежнему результат счета с постоянным малым шагом) представляет правильное решение задачи.
Решение Y2 стартует с увеличенным в два раза шагом после того, как быстрые изменения закончились. Видно, что эта попытка неудачна. Отклонения катастрофически нарастают, и получившаяся функция Y2 не имеет ничего общего с истинным решением, которое мы наблюдаем как функцию Y на этом же рисунке.
Описанное поведение характерно для жестких систем, если использовать для численного интегрирования явные методы с фиксированным шагом, такие как rkfixed.
Потенциально жесткие системы вида (2) даже в периоды медленных изменений имеют быстроменяющуюся составляющую. Как только возникают заметные отклонения от равновесного состояния F(X) - в силу случайных воздействий в реальном объекте или из-за вычислительных погрешностей в математической модели - появляется необходимость адекватно описывать очень быстрые изменения, т.е. вести счет с малым шагом. Фактически интегрирование нужно вести все время с самым малым шагом, даже несмотря на то, что решение изменяется очень медленно. Это парадокс жестких систем.
Рассмотрим еще два варианта, иллюстрирующих полную неудачу при решении, если шаг выбран неправильно. Увеличим параметр a вдвое, оставив шаг таким же как в исходном варианте.
Параметр a, определяющий скорость релаксации, увеличен вдвое. Остальные параметры такие же как в примере «Исходный вариант решения».
Используется новое, увеличенное значение параметра a (см. «Вариант 3») :
Другие параметры задачи такие же как в пункте «Вариант 1» для решения, стартующего в промежуточной точке:
Результат расчетов в пунктах «Вариант 3» и «Вариант 4» показаны на рис. 4. Видно, что с самого начала возникают катастрофические отклонения от истинного решения. Чтобы показать рост отклонений, то есть меру неустойчивости решения, пришлось выбрать очень большой масштаб по вертикальной оси.
Рис. 4. Неустойчивость решения при больших значениях параметра a
Ясно, что величина шага должна быть значительно уменьшена. Вопросы устойчивости численных схем интегрирования обсуждаются ниже в разделе «Явные и неявные методы».
Сделаем еще одно замечание относительно выбора шага при использовании rkfixed в общем случае. Априори, до решения задачи, трудно выбрать такой шаг интегрирования, который обеспечивал бы необходимую точность численного решения. Поэтому на практике поступают следующим образом. Проводят интегрирование при некотором разумно выбранном шаге и запоминают полученные результаты. Затем уменьшают шаг вдвое и снова проводят интегрирование. Разность между новыми и старыми значениями искомой функции является хорошей оценкой погрешности.
В целом приведенные выше примеры показывают, что метод rkfixed мало пригоден для решения жестких задач, и мы перейдем далее к анализу метода с адаптирующимся шагом.
Идея этого метода состоит в том, чтобы автоматически изменять, подстраивать величину шага к характеру изменения искомой функции, как мы пытались это сделать "вручную" в разделе «Метод rkfixed. Неустойчивость решения».
Список параметров процедуры таков:
где дополнительные (по сравнению с rkfixed) три параметра имеют следующий смысл:
· acc (accuraccy level) - точность интегрирования; если необходимо, процедура уменьшает шаг, чтобы ошибка на каждом шаге не превосходила эту величину,
· Nmax - максимальное число точек, в которых должно быть выведено решение; если задать Nmax небольшим целым числом, то выведенные точки будут сосредоточены у левого конца отрезка.
· save_int -минимальный шаг, с которым выводятся результаты интегрирования.
Если положить save_int=0 и Nmax="очень большое число", то выходной массив будет содержать все точки, в которых было вычислено решение, даже очень близко расположенные друг к другу.
Как и ранее, примем
Рис. 5. Интегрирование жесткой задачи методом rkadapt
Число шагов интегрирования - 29.
Видно, что адаптирующийся алгоритм справляется с этой жесткой задачей (Рис. 5). Обработав быстропеременное решение вблизи начала координат, программа пытается увеличить шаг, но подчиняясь требованию обеспечить нужную точность, оставляет шаг относительно малым. Благодаря этому не происходит срыва решения из-за численной неустойчивости явного метода, как это случилось в примере с грубым ручным увеличением шага (см. «Метод rkfixed. Неустойчивость решения»).
Однако здесь же проявляется недостаток метода rkadapt при работе с жесткой задачей: в области медленного изменения шаг должен оставаться малым, практически таким же как для переходного процесса в начале координат. Если представить себе, что область интегрирования составит не один период (2*p в нашем примере), а сотни или тысячи периодов медленного изменения, то значение этого недостатка станет очевидным.
Покажем теперь, насколько эффективнее работает rkadapt c обычной, не жесткой задачей. Чтобы сравнение было убедительным, проинтегрируем дифференциальное уравнение второго порядка
очевидным решением которого является та же функция, что и для жесткой задачи, а именно:
Представим дифференциальное уравнение второго порядка как систему двух дифференциальных уравнений первого порядка и обратимся к методу rkadapt :
Рис. 6. Интегрирование обычной, не
жесткой задачи методом rkadapt
Число шагов интегрирования – 7.
Сопоставьте Рис. 5 и Рис. 6. В обычной задаче потребовалось гораздо меньше шагов интегрирования (7 против 29). Если в жесткой задаче параметр a сделать еще большим, например, равным 40, результат будет 7 против 80, то есть еще более показательным.
Сопоставим теперь работу явного метода с адаптацией rkadapt и специального метода для жестких уравнений stiffr при интегрировании модельного дифференциального уравнения (2). Условия сравнения будут действительно жесткими - при очень больших значениях параметра a:
Для обращения к процедуре stiffr требуется составить матрицу Якоби (см. пояснения в разделе «Явные и неявные методы»):
Смысл других параметров вызова процедур уже обсуждался.
Рис. 7. Интегрирование жесткой
задачи с очень большим значением параметра a методом rkadapt
Число шагов интегрирования для rkadapt очень большое (1685), поэтому точки на графике сливаются в жирную линию.
Рис. 8. Интегрирование жесткой
задачи с очень большим значением параметра a методом stiffr
Число шагов интегрирования для stiffr относительно невелико (73).
Результаты вычислений представлены на Рис. 7- для метода rkadapt и на Рис. 8 - для stiffr.
Оба метода справились с обработкой быстрого переходного процесса в начале координат, но эффективность счета в области медленных изменений разительно отличается - сопоставьте приведенные выше цифры для числа шагов интегрирования (1685 и 73).
Метод rkadapt бесполезно тратит машинное время и память для расчета медленных изменений с очень малым шагом, в то время как метод stiffr оказывается действительно способным подстраиваться к характеру изменения искомой функции.
В разделе «Явные и неявные методы» мы обсудим сущность этих различий, учитывая, что сравниваемые процедуры находятся по разные стороны главного водораздела в численном анализе - между явными и неявными методами.
Приведем еще один пример эффективной работы метода stiffr - теперь для жесткой системы дифференциальных уравнений [1]:
Параметры принимают значения в соответствии с неравенствами a1>>1, a0<<1.
Этот пример взят из области химической кинетики. Уравнение (7) описывает медленное изменение концентрации 0-компонента за счет автокатализа (первое слагаемое справа) и участия в реакции с 1-компонентом. Согласно уравнению (8), концентрация 1-компонента уменьшается в процессе быстрого распада (первое слагаемое справа) и возрастает благодаря бимолекулярной реакции.
Запишем правую часть системы уравнений в матричной форме, обозначая независимую переменную (время) через T:
Матрица Якоби и начальные условия:
Рис. 9.Интегрирование жесткой задачи
химической кинетики методом stiffr
Изменение концентрации компонентов во времени представлено на Рис. 9. Видно , что решение развивается на большом интервале времени, но имеется также стадия быстрого изменения, когда концентрация компонента "1" резко убывает на малом промежутке времени порядка 1/a1. Метод stiffr воспроизводит это быстрое, практически скачкообразное изменение Y1 вблизи начала координат благодаря счету с очень малым шагом, а затем обеспечивает экономные вычисления с большим шагом в области относительно медленного изменения концентраций.
Если стадия быстрого изменения Y1 не представляет интереса, то математическую модель можно серьезно упростить. Вместо дифференциального уравнения (8) , отвечающего за изменение концентрации компонента "1", следует использовать так называемое адиабатическое (равновесное) приближение, приравняв правую часть нулю:
откуда
На Рис. 9 видно совпадение "точного" решения (10) (ромбы) и адиабатического приближения (11) (сплошная линия) везде, за исключением очень короткого начального периода. Метод адиабатического приближения эффективно используется при анализе сложных нелинейных систем [1].
Рассмотренный пример с системой дифференциальных уравнений позволяет коротко коснуться общего определения понятия жесткости [2]. Говорят, что система дифференциальных уравнений является жесткой, если собственные числа матрицы Якоби сильно различаются одно от другого.
В связи с этим определением мы имеем хороший повод вычислить эти собственные значения, применив встроенную функцию Mathcad eigenvals
Выпишем матрицу Якоби, определенную через производные правой части по зависимым переменным (см. (9) ):
Найдем локальные значения собственных чисел и построим график (Рис. 10) их изменения на интервале интегрирования.
Рис. 10. Собственные числа матрицы
Якоби для жесткой задачи химической кинетики
Видно, что на всем интервале решения собственные числа различаются более чем на порядок, следовательно это жесткая задача. Несложно показать, что по крайней мере локально (т.е. на небольших интервалах независимой переменной) малые и большие собственные значения определяют соответственно медленное и быстрое решения в виде экспонент: exp(EV0 ×T) и exp(EV1 ×T).
В нашей задаче наибольшее по модулю собственное значение равно примерно параметру a1: EV1 » -a1.
Рассмотрим простейшее модельное уравнение для жесткой задачи:
с начальным условием y(0) = 1 и очевидным решением y = exp(-ax). Свойства жесткости будут проявляться, если a>>1 и область интегрирования достаточно протяженная по сравнению с областью быстрых изменений 1/a.
Выпишем конечно-разностное представление этого уравнения на интервале
двумя различными способами:
Явная схема:
Неявная схема:
Величина yj подлежит определению, а значение yj-1 известно из начального условия или с предыдущего шага решения.
Принципиальное отличие двух численных схем состоит в способе вычисления правой части:
v в явной (explicit) схеме используется известное значение yj-1 (значение из "прошлого")
v в неявной (implicit) схеме - неизвестное, подлежащее определению значение yj (значение из "будущего")
Теперь необходимо получить из этих аппроксимаций расчетные формулы для искомой величины yj.
В случае явной формулы (14) никаких проблем не возникает, так как неизвестное присутствует только в левой части. Схема и называется явной потому, что фактически мы уже имеем раскрытое относительно yj выражение. После элементарных преобразований получим:
В случае неявной схемы (15)
Хотя благодаря линейности исходного дифференциального уравнения (13) выражение для искомой величины в нашем простом примере легко получить, и не обращаясь к решателю solve, для неявной схемы действительно необходимо решать алгебраическое уравнение или систему уравнений (если интегрируется система дифференциальных уравнений), что превращается в сложную проблему, если уравнения нелинейны.
Однако этот недостаток неявной схемы искупается ее замечательным свойством - абсолютной устойчивостью при счете.
Мы продемонстрируем сказанное, решая жесткое уравнение по явной и неявной схеме:
Рис. 11. Решение жесткого уравнения
по явной и неявной схеме
Как видно из графика, неявная схема обеспечивает удовлетворительное (учитывая выбранный большой шаг интегрирования) согласие с точным решением, в то время как явная дает совершенно нефизичное осциллирующее решение с нарастающей амплитудой. Это наглядное проявление численной неустойчивости явной схемы при выбранном шаге интегрирования. Если существенно уменьшить шаг, то неустойчивости можно избежать.
Из конечно-разностной аппроксимации (16) следует условие устойчивости явной схемы:
Таким образом, расчет быстропеременных (a>>1) решений требует работы с очень малым шагом (1/a), иначе результат будет не просто неточным, но не имеющим ничего общего с действительностью (см. пилообразное решение на Рис. 11). Для жестких задач счет с таким малым шагом должен продолжаться и в области слабых изменений (для x >>1/a), как мы видели это на приведенных выше примерах с явными методами rkfixed, rkadapt.
Большие затраты машинного времени, опасность возникновения трудно диагностируемых паразитных осцилляций, эстетически неприятная необходимость считать с очень малым шагом медленно меняющиеся функции - все это делает явные методы мало пригодными для решения жестких задач.
Неявная схема, как это непосредственно видно из (17), устойчива при любых (положительных) значениях шага. Поэтому жесткие задачи предпочитают решать по различным неявным схемам.
Для исследовательской работы, когда расчеты приходится проводить, еще не располагая информацией об особенностях поведения моделируемой системы, важным преимуществом неявных схем является получение разумных, физически правильных результатов даже при грубом счете с большим шагом.
Эффективность неявных схем объясняется тем, что при прогнозировании поведения функции используется информация из "будущего" - правая часть дифференциального уравнения выражается через еще неизвестные значения зависимой переменной на следующем шаге интегрирования (см. (15)). Неявный метод связан с попыткой заглянуть в будущее. Этим определяются и достоинства метода, и сложности при его реализации.
Нам осталось обсудить еще одну проблему, возникающую при обращении к встроенным функциям stiffr, stiffb для интегрирования жестких уравнений, а именно, необходимость задавать выражение для матрицы Якоби J(x,y) от правой части f(x,y) дифференциального уравнения y' = f(x,y) (см. уравнения (6), (9), (12)).
Это связано с обсуждавшейся выше необходимостью решать алгебраические уравнения на каждом шаге интегрирования по неявной схеме. Если дифференциальное уравнение нелинейное, то проводят линеаризацию правой части посредством ее разложения в ряд в окрестности начальной точки отрезка интегрирования:
Поскольку все величины в точке (x0, y0) известны, получается линейное уравнение относительно неизвестной величины y, как в нашем простейшем примере (17).
Если решается система дифференциальных уравнений порядка (N+1), то функцию f, приращение у, производную по х следует рассматривать как векторы-столбцы:
а производную по у как матрицу:
В справочной системе Mathcad матрица Якоби рассматривается как прямоугольная матрица, получающаяся при объединении (augment) вектора-столбца производных от правой части по x и матрицы производных по у. Такую матрицу необходимо сформировать при обращении к функциям stiffr, stiffb.
В математической литературе матрица Якоби чаще понимается как квадратная матрица производных правой части ДУ по зависимой переменной y :
Разложение правой части приспосабливают следующим образом к итерационной процедуре отыскания решения y :
где y(x) - значение y в правой точке отрезка интегрирования, взятое с предыдущей итерации, у - отыскиваемое новое значение.
Жесткие дифференциальные уравнения часто встречаются при математическом моделировании технических систем (см., например [3-5]).
Практически любое сколько-нибудь сложное устройство или установка характеризуется несколькими сильно различающимися временными постоянным. Другими словами, имеют место быстропротекающие переходные процессы, которые возникают в пусковые периоды или в любой момент времени вследствие изменений режима, случайных возмущений, ошибок и т.п. С другой стороны, в базовом режиме основные параметры системы обычно стационарны или медленно изменяются на больших промежутках времени.
Однако, несмотря на практическую важность вопроса, студенту или инженеру трудно найти доступное описание сущности этого феномена. Публикация восполняет в некоторой степени этот пробел и на простых примерах, с помощью средств, предоставляемых Mathcad' ом, демонстрирует эффективные способы решения таких задач.
Материалы публикации могут быть полезны для пользователей Mathcad'а, поскольку его справочная система почти ничего не содержит по обсуждаемой проблеме.
1. Хакен. Г. Синергетика. Иерархии неустойчивостей в самоорганизующихся системах и устройствах. M: Мир. 1985. 423 с.
2. Современные численные методы решения обыкновенных дифференциальных уравнений. Редакторы Дж. Холл и Дж. Уатт. Москва, изд-во Мир. 1979. С.312 (Modern Numerical Methods for Ordinary Differential Equations. Edited by G. Hall and J.M. Watt. Clarendon Press. Oxford. 1976/
3. Солодов А.П. Расчетные модели теплообмена при контактной
конденсации //Теплоэнергетика. 1990. № 10.
4. Солодов А.П. Неустойчивость и межфазная турбулентность в системе жидкость-пар.//Двухфазные течения: Труды Первой Российской национальной конференции по теплообмену. М.: Изд-во МЭИ. 1994.
5. Закревский С.Л, Солодов А.П. Динамическая модель подогревателя низкого давления смешивающего типа // Теплоэнергетика. 1998. №7. С. 48-51.