Дифференциальные уравнения

Типичный пример — численное решение в среде Mathcad системы алгебраических уравнений (см. главу 2): сначала задается начальное приближение к корню (ведь корней может быть много), потом за ключевым словом Given (дано) в естественном виде пишется сама система уравнений, после которой помещается функция Find, возвращающая значения своих аргументов, превращающих ранее записанные уравнения (до слова Given) в тождества. Функция Find также аномальна: она возвращает значения, зависящие не только от значений аргументов, но и от того, что вблизи нее находится. Это, конечно, грубейшее нарушение правил построения функциональных зависимостей, и это уже отмечено в разд. 1.2.1. Но разработчикам Mathcad пришлось пойти на это ради того, чтобы алгебраические уравнения в Mathcad-документе хранились в естественной форме. При создании математических пакетов типа Mathcad приходится решать не только чисто математические проблемы (разработка алгоритмов решения и реализация их средствами вычислительной математики), но и проблемы единства формы представления задачи (интерфейс пользователя) и методов ее решения (содержание задачи).

Функция sbval не решает краевую задачу, а только находит недостающие значения на краю отрезка. После этого краевая задача переходит в задачу с начальными условиями (задача Коши), которая для проверки и решается так, как показано внизу рис. 5.11.

Для решения краевой задачи в среде Mathcad есть еще одна функция — bvaifit. Она используется в тех случаях, когда нет всей необходимой информации для функции sbval, но известно решение задачи в промежуточной точке. Функция bvaifit решает задачу с двумя граничными точками, начиная с конечных точек и следуя траекториям решения и его производным в промежуточных точках:

bvaifit(xl, х2, tнач, tкон, tf, f, loadl, load2, score)

где:

- x1 — вектор предполагаемых начальных значений, не определенных в точке xl;

- х2 — то же для точки х2;

- tнач, tкон — конечные точки интервала, в котором должно быть вычислено решение дифференциального уравнения;

tf — точка между tнач и tкон, в которой траектории решений, начатые с tнач, и траектории, начатые с tкон, должны совпадать;

- f(t, x) — векторная функция, состоящая из п элементов и содержащая первые производные неизвестной функции;

- loadl (tнач, x1) — векторная функция, п элементов которой соответствуют значениям п неизвестных функций в точке x1. Некоторые из этих значений будут постоянными, заданными начальными условиями. Другие будут неизвестны вначале, но будут найдены. Если значение неизвестно, то следует использовать соответствующее предполагаемое значение из вектора x1;

- load2 (tкон, x2) — аналог load1, но для значений п неизвестных функций в точке х2;

- score (tf, x)— векторная функция, состоящая из п элементов. Эта функция применяется для того, чтобы задать, как решения должны совпадать в точке tf. Обычно нужно определить score (tf, x) :=x, чтобы все решения неизвестных функций совпадали в точке tf.

Рис. 5.12. Решение краевой задачи об эпидемии средствами программирования

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

На рис. 5.12 зафиксирован еще один подход к решению краевой задачи "Эпидемия", без использования сложных и довольно туманных встроенных функций sbval И bvalfit.

На рис. 5.12 с помощью средств программирования Mathcad автоматизированы последовательные приближения при реализации метода стрельбы и метода половинного деления, "ручной" вариант которых был показан на рис. 5.10.

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

Казалось бы, что подстраиваемая (см. рис. 5.8) функция odesolve, введенная в Mathcad 2000 и расширенная в Mathcad 2001 и 14, должна "похоронить" все остальные функции, оставив их только для совместимости — для того, чтобы Mathcad-документ, написанный в 6-й версии, например, работал и в 14-й. Но это не совсем так. Основной недостаток функции odesolve в том, что ее нельзя вставить в программы. Этому мешает ключевое слово Given. И если возникнет необходимость программно реализовать метод стрельбы (половинного деления) для решения краевой задачи, то тут нужно будет вернуться к функции rkfixed (см. рис. 5.12), например, несмотря на все другие вышеописанные преимущества функции odesolve.

Но и "первородная" функция rkfixed и другие подобные функции имеют свои недостатки. Основной из них— невозможность создания функции пользователя для многократного вызова этой функции при формировании различного вида зависимостей— новых пользовательских функций. Такую задачу можно решить только за счет некоего обмана, пакета Mathcad. Так на рис. 5.14 (первая его половина) сделана попытка ввода в функцию f, хранящую правые части системы двух дифференциальных уравнений развития эпидемии15, дополнительного аргумента Пр— f(t, x, Пр). Но ничего из этого не вышло. Было выдано сообщение об ошибке "This function needs more arguments" ("Эта функция требует больше аргументов"), хотя здесь аргументов должно быть наоборот меньше.

Рис. 5.13. Фазовый портрет задачи об эпидемии

Положение спасло добавление в функцию-вектор F(t, x) дополнительного не аргумента, а элемента (компонента вектора) с нулевым (пустым) значением. После этого матрицу решения задачи z стало возможным заменить на функцию z (Пр) и построить "одним махом", например, график изменения максимального количества больных (вершина "горы", показанной на рис. 5.2 и 5.7) от значения коэффициента профилактики (см. рис. 5.14 внизу).

В среде Mathcad 12/13/14 появилась возможность создания локальной функции, что позволяет нам более естественно решить задачу о создании функции пользователя на базе встроенной rkfixed (рис. 5.15).

Рис. 5.14. Создание функции пользователя на базе встроенной функции rkfixed

Встроенные инструменты Mathcad имеют еще один общий недостаток, связанный с тем, что при возникновении ошибки счета не дается полной информации о причинах сбоя. Из-за этого часто приходится отказываться от встроенных средств и возвращаться к пользовательским, "доморощенным" функциям. Так "в пику" встроенным в Mathcad функциям Find и Minimize, мы описали функции поиска корней систем уравнений и минимумов функции средствами программирования (см. рис. 3.27, например), которые при всей их примитивности имеют одно преимущество — они выдают траекторию пути к решению.

Рис. 5.15. Создание функции пользователя на базе встроенной функции rkfixed с использованием локальной функции

Встроенные в Mathcad средства решения дифференциальных уравнений также обладают этим недостатком. Очень часто решения систем дифференциальных уравнений прерываются сообщением об ошибке, связанной с переполнением предела счета (10^307). Такое случилось с нашей задачей об эпидемии, когда значения коэффициента профилактики пр заменили с 0.0001 на 0.00885, например. При этом не ясно, какое именно уравнение системы "зашкалило", как нужно его подправить или как следует изменить пределы интегрирования, чтобы получить хотя бы упрощенное или укороченное решение для его последующего анализа. Такая ситуация отображена на рис. 5.16, где (первая половина) функция rkfixed дает сбой.

Рис. 5.16. Замена встроенной функции rkfixed на пользовательскую

Во второй половине рис. 5.16 встроенная функция rkfixed заменена на программу, реализующую метод Рунге — Кутты и возвращающую при ошибке частичное решение задачи за счет операторов обработки ошибок on error, вставленных в программу. Это позволяет вести интегрирование не до конечной точки (у нас это 3), а до последней допустимой. Встроенная же функция rkfixed излишне категорична — она возвращает "все или ничего".

В среде Mathcad 13/14 встроенные средства отладки позволяют провести трассировку и встроенных средств решения ОДУ (рис. 5.17).

Рис. 5.17. Трассировка встроенной функции odesolve