Численное решение уравнений биокинетики в курсах «Общая экология» и «Моделирование биологических процессов»

О системе MATLAB

О численных методах

Список используемых сокращений и соглашение об обозначениях

О терминах «граничные», «краевые» и «начальные условия»

О системе MATLAB

Сферы применения персональных компьютеров поистине неисчерпаемы и охватывают подготовку и обработку текстов, финансово-экономические расчеты, ведение баз данных, подготовку проектных документов, работу в глобальной сети Internet и т.д. Среди этого обилия применений видное место занимают математические и научно-технические расчеты – та область, ради которой компьютеры (от слова computer – вычислитель) были изначально созданы [Дьяконов и Абраменкова, 1999]. Для естественно-научных специальностей знание Фортрана или Бейсика длительное время представлялось вполне достаточным. Однако стремительное расширение сферы применения компьютеров неизбежно породило много новых систем программирования, гораздо лучше приспособленных к новым обстоятельствам и возможностям их практического использования. Так, для управления совместной работой компьютера с контрольно-измерительной аппаратурой используются Ассемблер и С; разработаны самые различные системы для поиска, анализа и визуализации информации [Кондрашов и Королев, 2002]. Научное программирование также претерпевает серьезную трансформацию: развиваются интегрированные среды, основанные на алгоритмических языках, растет применение универсальных математических систем (Maple, Mathematica, MATLAB, MatCad и др.) [Говорухин и Цибулин, 2001], среди которых MATLAB – наиболее развитая система для научно-технических расчетов [Кондрашов и Королев, 2002].

Систему MATLAB (сокращение от MATrix LABoratory – матричная лаборатория [Кондрашов и Королев, 2002]) разработал C.B. Moler, и с конца 70-х годов она широко использовалась на больших ЭВМ. В начале 80-х годов John Litl из фирмы The MathWorks, Inc. (США, г. Нейтик, шт. Массачусетс), разработал версию системы РС MATLAB для IBM PC [Потемкин, 1997]. К расширению системы были привлечены крупнейшие научные школы мира в области математики, программирования и естествознания [Дьяконов, 2001]. MATLAB является интерактивной системой для выполнения инженерных или научных расчетов, ориентированной на работу с массивами данных [Потемкин, 1997], которая интег-
рирует в единое целое [Потемкин, 2001]:

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

Экономисты используют MATLAB для исследования, анализа, планирования и прог-
нозирования экономической динамики [Цисарь, 2002]. MATLAB лидирует в военно-промышленном комплексе, аэрокосмической индустрии и в автомобилестроении [Дьяконов и Абраменкова, 2002]. С равным успехом она применима к задачам математики, физики, биологии 2, электротехники, радиотехники и т.д. [Дьяконов и Абраменкова, 1999].

MATLAB представляет собой одну из старейших, тщательно проработанных и апробированных (а потому – надежных) систем. Это уникальная коллекция реализаций численных методов, созданных за последние десятилетия, вобравшая в себя также опыт, правила и методы математических вычислений, накопленные за тысячи лет развития математики [Дьяконов и Абраменкова, 1999]. В 1998 г. систему использовали более 500 000 легально зарегистрированных пользователей [Дьяконов, 2001]. А к 1999 г. MATLAB стал инструментом почти всех инженеров и прикладных математиков [Мэтьюз и Финк, 2001].

Система MATLAB – это одновременно и операционная среда и язык программирования [Потемкин, 1998]. Важным достоинством системы является ее открытость и расширяемость. Большинство команд и функций системы реализованы в виде так называемых М-файлов (с расширением .m) и файлов на языке С, причем файлы доступны для модификации. Пользователю предоставлена возможность создавать не только отдельные файлы, но и библиотеки файлов для реализации специфических задач. Благодаря текстовому формату М-файлов пользователь может ввести в систему любую новую команду, оператор или функцию и пользоваться затем ими столь же просто, как и заведомо встроенными операторами или функциями [Дьяконов и Абраменкова, 1999].

Язык, используемый в системе MATLAB, можно сравнить с языком BASIC по простоте его применения [Потемкин, 1997]. По сути, он напоминает BASIC с примесью Фортрана и Паскаля. Но, в то же время, в части программирования математических вычислений язык системы MATLAB намного богаче любого универсального языка программирования высокого уровня [Потемкин, 1998]. MATLAB – это язык четвертого поколения, а теперь – в связи с возможностями создания GUI-приложений – и язык, который можно отнести к пятому поколению. Он весьма удобен для проведения анализа, для моделирования и для построения полноценных самодостаточных конечных программных продуктов [Потемкин, 2001].

Основным объектом является массив, для которого не требуется указывать размерность явно. Это позволяет решать многие вычислительные задачи, связанные с векторно-матричными формулировками, существенно сокращая время, которое понадобилось бы для программирования на языках типа С или Fortran [Потемкин, 1998] (опыт разработки финансовых приложений показывает, что использование системы MATLAB экономит до 90% времени по сравнению с программированием на языке С [Потемкин, 2001]). Также система работает с многомерными массивами, массивами записей и ячеек, с матрицами, в том числе и разреженными, позволяет проводить анализ и обработку данных, включая аппроксимацию и интерполяцию, численное интегрирование, решение систем обыкновенных дифференциальных уравнений, вычисление минимумов и нулей функций, преобразование Фурье, свертку и фильтрацию и т.д. и т.п. [Васильков и Василькова, 1999].

MATLAB прекрасно интегрируется с Microsoft Word и Excel. Надстройка Excel Link, поставляемая вместе с MATLAB, существенно расширяет возможности Excel, обеспечивая доступ пользователя к функциям MATLAB. Подготовка данных осуществляется непосредственно в электронных таблицах, а обращение к функциям производится либо из ячеек рабочего листа, либо в модуле, написанном на Visual Basic [Ануфриев, 2002].

Возможности MATLAB весьма обширны, а по скорости выполнения задач эта система превосходит другие подобные системы (при прочих равных условиях скорость почти на порядок выше, чем у системы Mathcad) [Дьяконов и Абраменкова, 1999]. Однако возможности системы MATLAB еще сильнее расширяются при использовании дополнительных пакетов программ, так называемых тулбоксов (TOOLBOX – ящик с инструментами) [Егоренков с соавт., 1994].

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

Значимым недостатком системы MATLAB можно назвать то, что это эта система является коммерческим продуктом, а её приобретение стоит довольно дорого, особенно в российских реалиях. Впрочем, «студенческие» версии программы с ограниченным функционалом может позволить себе приобрести практически каждый. Главными свободно распространяемыми аналогами MATLAB являются R и Python. Поскольку программы на этих языках, реализующие различные численные методы, развиваются широким программистским сообществом, а не одной коммерческой корпорацией, в некоторых случаях доступный набор возможных опций решения задач на них может оказаться шире, чем в MATLAB. Однако качество работы таких решателей не всегда может быть достаточно высоким.

О численных методах

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

Аналитические решения позволяют получить представление об общем решении, однако для лучшего понимания эволюции решений исследуемого дифференциального уравнения полезно нарисовать графики некоторых его частных решений [Шампайн и др., 2009, с. 12]. Часто эти решения содержат специальные функции – в качестве примера можем привести классическое уравнение ферментативной кинетики – уравнение Михаэлиса-Ментен

dS/dt = –Vm × S/(KM + S),

где S – концентрация вещества, на которое действует фермент; t – время; (тогда dS/dt – скорость ферментативной реакции); Vm – ее максимальная скорость; KM – константа Михаэлиса (равная такой S, при которой скорость реакции достигает ровно половину максимальной). Решение этого уравнения (т.е. S как функция t):

S(t) = KM × W(KM–1 × e–[t + C1] × Vm/KM),

здесь W(х) – W-функция Ламберта 3, C1 – константа, удовлетворяющая уравнению

Sо = KM × W(KM–1 × e–C1 × Vm/KM),

где Sо = S(0) – начальная концентрация вещества.

Для получения значений специальных функций, присутствующих в записи аналитического решения, необходимо применить существующие вычислительные процедуры [Шампайн и др., 2009, с. 12]. Более того, пусть нам пришлось решать какое-нибудь совсем простое уравнение, решение которого выражается не через специальные, а через элементарные функции. Но если это будут элементарные трансцендентные функции (т.е. типа экспоненты или синуса), то и здесь ситуация аналогичная – для получения значений элементарных трансцендентных функций тоже необходимо применить численные процедуры 4. Действительно, пусть решение выражается через cos(x). Как нам посчитать, например, cos(0.12)? Вы скажете, что какой-нибудь эксель или даже просто калькулятор выдаст ответ почти мгновенно? Но это лишь потому, что кто-то ранее написал соответствующую программу, позволяющую при помощи некоторого численного метода определить данное значение.

Однако если в подобных ситуациях приходится использовать численные методы, почему бы нам не отказаться от предварительного получения решения в аналитической форме и всегда находить его численно [Шампайн и др., 2009, с. 12]? Справедливости ради отметим, что у оппонентов есть на это возражения, впрочем, как сможет увидеть вдумчивый читатель – достаточно нелепые. Состоят они в следующем.

Первое возражение. Непосредственное использование численных методов для получения решений может быть приемлемым подходом при исследовании частных случаев, однако в общем случае применение численных методов имеет некоторые ограничения. Численное определение значений специальных функций (например, функций Эйри и Бесселя) может вызвать определенные трудности, т.к. эти функции могут иметь сингулярности или осциллировать с высокой частотой. Если решение исследуемого уравнения характеризуется подобным поведением, или если мы заинтересованы в исследовании поведения этого решения при t→¥, мы не сможем воспользоваться численными методами для непосредственного нахождения этого решения. Получение решения в аналитической форме позволяет изолировать указанные проблемы, в результате чего вычисление точного решения будет зависеть лишь от качества программ, предназначенных для нахождения значений специальных функций [Шампайн и др., 2009, с. 13]. Наши контраргументы таковы: во-первых, при решении реально возникающих на практике задач мы вряд ли будем заинтересованы в исследовании поведения решения при t→¥ (хотя максимальное значение интересующего нас t может быть весьма велико); во-вторых, если математическая модель правильно описывает биологический процесс, то никаких сингулярностей в решении быть не должно, потому что в реальных биологических процессах мы эти самые сингулярности не наблюдаем; в-третьих, решение иных проблем (типа осцилляции решения с высокой частотой) будет зависеть лишь от качества программ, предназначенных для нахождения численного решения дифференциальных уравнений. И если предполагается, что высокого качества программ можно добиться в случае численного определения значений специальных функций, то почему такого же качества нельзя будет добиться в случае численного решения уравнений?

Второе возражение. Поскольку численные методы позволяют найти лишь частное решение, не всегда можно ответить на вопрос о том, как это решение зависит от параметров [Шампайн и др., 2009, с. 13]. Но почему же? Что мешает прорешать задачу многократно – при разных значениях параметров!?

Список используемых сокращений
и соглашение об обозначениях

ДиРУ – дифференциально-разностные уравнения;

ЗаСГУ – задача с граничными 5 условиями;

ЗаСНУ – задача с начальными условиями (задача Коши);

МаМ – математическое моделирование;

МЭНЯ – метод Эйлера (неявный);

МЭЯ – метод Эйлера (явный);

НеЗ – нежесткие задачи;

ОДУ – обыкновенные дифференциальные уравнения;

ПроЖИ – программы жесткого интегрирования;

САДУ – система алгебро-дифференциальных уравнений;

СЛАУ – система линейных алгебраических уравнений;

УЧП – уравнения в частных производных

ФОС – функция обработки событий;

BDF – «формулы дифференцирования назад» (backward differentiation formulas).

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

• скалярные величины – латинские буквы (чаще всего 6 строчные), выделенные курсивом (например, t, х);
• элементы векторов будут обозначаться так же, но с добавлением подстрочного индекса – обычно i или j (например, bi , cj );
• сами вектора будем обозначать соответствующими строчными латинскими буквами, выделенными полужирным шрифтом (например, b, c);
• матрицы (прямоугольные таблицы чисел) будем обозначать прописными латинскими буквами;
• операторы обычно будут помечаться «крышечкой» (например, Â, Ĉ), за исключением широко распространенных давно используемых обозначений (типа ∇-оператора Гамильтона).

О терминах «граничные», «краевые»
и «начальные условия»

При решении задач из какой-либо области науки математическими методами необходимо написать:

а) уравнение, которому удовлетворяет искомая функция, описывающая исследуемое явление;

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

Дело в том, что интересующие нас явления обычно имеют однозначный характер, в то время как описывающие их уравнения имеют множество решений. Дополнительные условия позволяют вычленить лишь одно решение, описывающее конкретное явление [Соболев, 1954, с. 28–29; Кошляков и др., 1970, с. 11–12; Арсенин, 1984, с. 32–33; Владимиров, 1988, с. 71–72].

При этом выделяют граничные условия, характеризующие особенности взаимодействия граничной поверхности тела с окружающей средой, и начальные условия, характеризующие состояние тела в исходный момент времени. Эти условия в совокупности называются краевыми условиями: начальное условие является временным краевым условием, а граничные – пространственным краевым условием. Дифференциальное уравнение вместе с краевыми условиями составляет краевую задачу. Разумеется, для установившегося процесса необходимость задавать начальное условие отпадает, и в этом случае краевая задача будет состоять из уравнения и граничных условий [Карташов, 1985, с. 49]. Из числа классических учебников, придерживающихся такой классификации, можно указать [Владимиров, 1988, с. 71–72]. Эта терминология нашла отражение и в прикладных работах в области МаМ почвенных процессов. В частности, ей следуют Д.С. Орлов и др. [1987, с. 146–147].

Однако вышеприведенная стройная терминологическая система присутствует лишь в части литературы, причем, как представляется, в не слишком большой ее части. Иногда начальные условия даже в пределах одной и той же книги то считаются краевыми, то – нет, как это имеет место, например, в [Бицадзе, 1982, с. 30, 39, 177; Самарский и др., 1987, с. 18, 50]. Тем более это может происходить в разных книгах одного и того же автора – ср., например, превращение «начальных» условий в «краевые», происходящее в учебнике [Бицадзе, 1982] с четким разграничением между этими условиями в [Бицадзе, 1981, с. 37]. Впрочем, наше внимание сейчас занимает не столько сумятица в различении начальных и краевых условий (мы упомянули об этом лишь для того, чтобы показать, что даже по такому краеугольному вопросу «в товарищах согласья нет»), сколько краевые условия сами по себе и в их отношении к условиям граничным.

К сожалению, даже в классических учебниках по математической физике 7 царит полный разнобой в применении вышеперечисленных терминов. Видимо, это явилось следствием того, что терминология менялась на протяжении десятилетий. Так, в начале 20-х гг. ХХ в. В.А. Стеклов [1983, c. 57–58, 62] выделял «начальные» условия (по времени), тогда как дополнительные условия на пространственной границе называл «предельными» (изредка – «поверхностными»). Позднее, С.Л. Соболев [1954, с. 31–32], а вслед за ним Н.С. Кошляков и др. [1970, с. 254, 449] также говорили о начальных и предельных условиях, но последние называли одновременно еще и «граничными», и «краевыми» (давая синонимичные термины в скобках после основного). Позднее о том же самом писал В.Я. Арсенин [1984, с. 33], хотя у него термин «предельные условия» уже исчез: осталась только эквивалентность условий граничных и краевых вкупе с их отличием от начальных. И это также находит отражение в работах, посвященных МаМ в биологии. В частности, то, что условия на пространственной границе можно одновременно называть и «граничными», и «краевыми», явным образом указывается в [Ризниченко, 2011, с. 291; Плохотников, 2017, с. 58]. А Глаголев и Смагин [2005, с. 91, 103], хотя и не обращают явным образом внимание читателя на это, но исповедуют ту же эклектику: один из авторов в написанном им примере называет дополнительные условия на пространственной границе «краевыми», тогда как другой (в своем примере) – «граничными».

В результате вышеописанной терминологической чехарды, в современной литературе, использующей МаМ в виде УЧП, условия на пространственной границе называют как «граничными» (см., например, [Рубин и др., 1987, с. 75; Дружинин и Шишкин, 1989, с. 56; Братусь и др., 2010, с. 353; Глаголев и др., 2016]), так и «краевыми» – [Говорухин и Цибулин, 2001, с. 410; Братусь и др., 2010, с. 319; Акимов и др., 2017, с. 223]. При работе с ОДУ различные авторы также иногда называют условия на этой границе «краевыми» (при этом саму задачу то называют «краевой» – например, [Говорухин и Цибулин, 2001, с. 407–408; Глаголев и Смагин, 2005, с. 73; Акимов и др., 2017, с. 207; Плохотников, 2017, с. 48–49], то никак ее не называют [Братусь и др., 2010, с. 324–325]; впрочем, в старой литературе бывала и обратная ситуация: данную задачу называли «краевой», но сами условия никак не называли – [Хемминг, 1972, с. 229–232; Каханер и др., 1998, с. 387, 399]). Иногда же эти условия называют «граничными», а задачу, соответственно, либо тоже «граничной» [Кондрашов и Королев, 2002, с. 115; Шампайн и др., 2009, с. 536], либо – «краевой» [Ортега и Пул, 1986, с. 68; Хайрер и др., 1990, с. 109–110; Мэтьюз и Финк, 2001, с. 536], либо – никак [Вавилин и Васильев, 1979, с. 28, 31; Рубин и др., 1987, с. 76–77; Мюррей, 2009, с. 621].

Учитывая,

• во-первых, вышеизложенное (в настоящее время, фактически, нет общепринятого стандарта различия в употреблении терминов «граничные» и «краевые» условия);
• во-вторых, то, что при написании данной книги нам пришлось пользоваться многочисленными разнородными литературными источниками;
• и, в-третьих (главное!), что в данном томе мы почти везде будем рассматривать лишь ЗаСГУ для ОДУ, а в этом случае граничные условия действительно эквивалентны краевым даже со строгой точки зрения работ [Карташов, 1985, с. 49; Владимиров, 1988, с. 71–72],

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

---