Воронеж 2006

«ГОСУДАРСТВЕННОЕ ОБРАЗОВАТЕЛЬНОЕ УЧРЕЖДЕНИЕ

ВЫСШЕГО ПРОФЕССИОНАЛЬНОГО ОБРАЗОВАНИЯ

ВОРОНЕЖСКИЙ ГОСУДАРСТВЕННЫЙ ТЕХНИЧЕСКИЙ УНИВЕРСИТЕТ»

Кафедра автоматизированных и вычислительных систем

ЗАДАНИЕ

на курсовую работу

по дисциплине «Пакеты прикладных программ»

Тема «Решение систем дифференциальных уравнений в MathCad»

Студент группы ВМ – 052___________________Тютин Алексей Анатольевич

Перечень вопросов, подлежащих разработке:

-изучить основы MathCad;

-определить инструменты MathCad, наиболее подходящие для достижения поставленных целей;

-изучить механизм работы инструментов

В качестве методовисследования выступили изучение специализированной литературы, а также сбор информации во Всемирной Сети, как в наиболее динамично развивающейся информационной среде, чрезвычайно популярной среди пользователей MathCad.

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

Срок защиты курсовой работы________________________________________

Руководитель ________________ М.П. Носачева__

Подпись,

дата

Задание принял студент ________________ А.А. Тютин

Подпись,

дата

ЗАМЕЧАНИЯ РУКОВОДИТЕЛЯ

Содержание

Задание на курсовую работу…………………………………………………….2

Замечания руководителя…………………………………………………………3

Введение…………………………………………………………………………..5

1 Функции для решения дифференциальных уравнений……………………..6

2 Функция odesolve………………………………………………………………7

3 Функции rkfixed, Rkadapt, Bulstoer…………………………………………..11

4 Функции Stiffb, Stiffr………………………………………………………….15

5 Функций rkadapt, bulstoer, stiffb, stiffr……………………………………….17

6 Функции sbval, bvalfit……………………………………………………..…..19

Заключение……………………………………………………………………….23

Список литературы……………………………………………………………....24

Введение

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

Функции для решения дифференциальных уравнений

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

Некоторые из функций предназначены для быстрого нахождения конечных значений искомых переменных на правом конце выбранного интервала интегрирования. В таблице 1 и далее использованы сокращения: ДУ - дифференциальные уравнения, ОДУ - обыкновенные дифференциальные уравнения.

Перед описанием функций введем обозначения основных параметров, которые используются в качестве аргументов большинства встроенных функций:

у - (n х 1) - вектор результирующих переменных (n 1);

у0 - (nх 1) - вектор начальных значений переменных;

х, х1, х2 - аргумент, левая и правая границы его диапазона соответственно;

т - число точек, в которых находится решение внутри интервала (xl, х2);

D(x,y) - (n х 1) - вектор правых частей системы дифференциальных уравнений первого порядка, соответствующий первым производным вектора у. Этот вектор должен быть предварительно, до использования какой-либо функции, введен в виде выражения D(х,у): = (правые части уравнений).

При п = 1 решение ищется для одного дифференциального уравнения.

Для решения обыкновенных дифференциальных уравнений в системах MathCAD предусмотрены четыре функции, обеспечивающие:

odesolve(x, x2,[m]) - решение одного обыкновенного дифференциального уравнения методом Рунге-Кутта с постоянным (по умолчанию) или адаптивно вычисляемым системой шагом интегрирования;

rkfixed(y0,xl,x2,m,D) - решение уравнений методом Рунге-Кутта с постоянным шагом интегрирования, равным (х2 - xl)/m;

Rkadapt(yO,xl,x2,m,D) - решение уравнений методом Рунге-Кутта с шагом интегрирования, адаптивно выбираемым в зависимости от характера изменения у(х);

Bulstoer(yO,xl,x2,m,D) - решение уравнений методом Булирш-Штоера.

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

Функция odesolve

Встроенная функция odesolve введена в системы MathCAD, начиная с восьмой версии. Эта функция предназначена для решения отдельных обыкновенных дифференциальных уравнений, порядок которых n 1. Функция odesolve используется в составе вычислительных блоков для интегрирования дифференциальных уравнений, которые могут рассматриваться в качестве частного случая рассмотренных выше вычислительных блоков для решения алгебраических уравнений и задач оптимизации.

Вычислительные блоки для интегрирования дифференциальных уравнений (ВБДУ) имеют несколько отличительных особенностей:

• ВБДУ не нуждаются в задании начальных приближений;

• ВБДУ также начинается с ключевого слова Given, но за ним следуют) обыкновенное дифференциальное уравнение и его начальные (граничные для краевых задач) условия^-, последние выступают в качестве ограничений;

• ВБДУ завершается функцией odesolve, которая (в отличие от ранее рассмотренных вычислительных блоков) не является исполнимой и которой следует присваивать имя;

• результат вычислений в ВБДУ непосредственно можно использовать лишь для построения графика и вывода значений получаемой функции при фиксированных значениях аргумента.

При использовании ВБДУ следует учитывать ряд правил, выступающих в качестве необходимых условий успешного решения задач:

1. все математические выражения, расположенные между ключевым словом Given и завершающей функцией odesolve (дифференциальное уравнение и начальные условия), должны иметь так называемые "жирные знаки равенства" - операторы отношения, вставка которых осуществляется клавишами Ctrl + = или с помощью математической палитры;

2. формат функции odesolve предполагает указание в числе ее аргументов переменной интегрирования, ее правого граничного значения и числа точек, в которых будет найдено решение, т. е. odesolve(x, x2,[m]), где квадратные скобки означают, что указание числа точек, в которых ищется решение, не обязательно;

3. при обозначении - производных внутри ВБДУ можно пользоваться двумя операторами - знаков дифференциала, которые вводятся с помощью математической палитры или с клавиатуры, а также надстрочными штрихами, каждый из которых вводится клавишами Ctrl + F7;

4. число ограничений при решении дифференциального уравнения порядка n 1 должно равняться п (начальное условие для функции и начальные условия ее (n - 1) производных). При нарушении равенства числа ограничений порядку ДУ, системой выдается сообщение о вариантах этого нарушения ("Too few initial conditions" при слишком, малом и "Too many initial conditions" - при слишком большом количестве начальных условий). Формат задания начальных условий для функции и ее производных единый, например, для функции: y(x1) = 5, где x1- левая граница интервала интегрирования. Параметры ДУ и ограничений могут быть заданы в символьной форме; в этом случае их числовые значения должны быть введены предварительно, до ВБДУ;

5. дифференциальное уравнение должно быть линейно относительно старшей производной.

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

На рисунке 1 представлено решение линейного ОДУ третьего порядка, описывающего поведение устойчивого динамического звена с постоянным входным воздействием. На рисунке 1,а ВБДУ используется в стандартном варианте - для построения графика и вывода отдельных значений функции. Перед ВБДУ задается значение входного воздействия, равное трем, и параметры ДУ. Фрагмент, изображенный на рисунке 1,б,содержит следующую часть рассматриваемого примера - формирование массивов данных по результатам использования ВБДУ. Целью представленных (рисунок 1,6) преобразований служит получение вектора моментов времени х с шагом 0.5 и соответствующего ему вектора Z значений функции; число элементов этих векторов равно 61. Матрица W, полученная с помощью встроенной функции augment, объединяет векторы х и Z. Отдельные значения функции -результата интегрирования ДУ, полученные ранее (рисунок 1,а), выводятся на рис. 3.20,б как соответствующие элементы второго столбца матрицы W. В левой нижней части рисунока 1,б выводится массив значений функции z для каждого момента времени .

Рисунок 1,а. Пример стандартного применения ВБДУ

Рисунок 1,б. Формирование массивов данных по результатам работы ВБДУ

Введение вектора Z и матрицы W позволяет использовать результаты, полученные ВБДУ, в последующих вычислениях. Приведенный фрагмент (рисунок 1,б) может быть использован для получения вектора значений функции для моментов времени, изменяющихся с произвольным шагом.

Пример решения нелинейного ДУ четвертого порядка с помощью ВБДУ, приведен на рисунке 2.

Рисунок 2. Решение нелинейного ДУ с применением ВБДУ

Приведенные примеры иллюстрируют простоту и удобствоприменения ВБДУ при интегрировании ОДУ. Использование встроенной функции odesolve по сравнению с другими встроенными функциями для решения задач интегрирования ОДУ, обладает заметными преимуществами в части решения ДУ высокого порядка. Другие встроенные функции (rkfixed, Rkadapt и проч.) предназначены для интегрирования систем ДУ и требуют, поэтому, предварительного преобразования ДУ высоких порядков в систему ДУ первого порядка (форма Коши). Кроме того, применение этих функций для решения одного ДУ первого порядка несколько сложнее, чем с помощью функции odesolve. Это будет проиллюстрировано ниже при рассмотрении функции rkfixed. Однако при решении задач с исходным заданием ДУ в форме Коши предпочтительнее использовать rkfixed, Rkadapt, Bulstoer. В таких задачах функцию odesolve можно применить лишь для ограниченного класса систем линейных дифференциальных уравнений, которые могут быть преобразованы в эквивалентное ДУ высокого порядка с константой в правой части. Поэтому при решении конкретной задачи следует предварительно проанализировать ее особенности и, после этого, выбрать наиболее рациональную встроенную функцию. Это, достаточно общая, рекомендация вызвана функциональной избыточностью средств, предоставляемых системами MathCAD, и относится к любым встроенным функциям, в том числе - и к функциям рассматриваемой категории.

Функции rkfixed, Rkadapt, Bulstoer

Встроенные функции rkfixed, Rkadapt, Bulstoer имеют однотипный состав аргументов и аналогичные условия применения. Эти функции предназначены для решения систем дифференциальных уравнений. Их можно использовать и для интегрирования одиночных ДУ (n=1), но в этом случае их применение имеет некоторые особенности, рассматриваемые ниже, формат этих функций приводился выше. При их использовании следует учесть, что предварительно в документе MathCAD должны быть введены: векторная функция D(x,y) с указанием ее аргументов и начальный вектор состояний. В аргументах функций rkfixed, Rkadapt, Bulstoer векторная функция D(x,y) указывается лишь своим первым символом.

Функции rkfixed и Rkadapt могут быть использованы в большинстве задач интегрирования систем п (n 1) дифференциальных уравнений. Функция rkfixed реализует метод Рунге-Кутта с постоянным шагом интегрирования, который задается в соответствии с динамикой изменения переменных системы. Величина шага зависит от заданных значений аргументов и составляет = (х2 - xl)/m, где x1, х2, и т - аргументы функции rkfixed. Функция Rkadapt также реализует метод Рунге-Кутта, но шаг интегрирования выбирается системой адаптивно, в зависимости от текущей изменчивости исследуемых переменных. При наличии в составе системы дифференциальных уравнений переменных с резко различными динамическими характеристиками или при существовании в переходных процессах, являющихся результатом интегрирования, участков с медленным и очень быстрым изменением переменных, функция rkfixed может привести к сильному увеличению времени вычислений. Это объясняется тем, что данная функция работает с постоянным шагом интегрирования, который в этом случае будет выбран очень мелким, исходя из требуемой точности аппроксимации быстрых составляющих. Этот шаг для медленно меняющихся составляющих в этих условиях будет неоправданно малым, что приведет к резкому увеличению числа шагов при попытке проанализировать весь переходный процесс. Для таких случаев, при наличии в составе системы переменных с резко различными динамическими характеристиками, более целесообразным является использование функции Rkadapt. Эта функция для участка быстро меняющихся переменных система выберет мелкий шаг интегрирования, а вычисление переменных на участках с медленным их изменением будет осуществляться крупными шагами, что сократит общее время решения. Независимо от выбранного функцией шага интегрирования результаты решения с применением функции Rkadapt выводятся в точках, число которых в интервале (xl, х2)задается значением т.

Функцию Bulstoer предпочтительнее использовать для интегрирования обыкновенных дифференциальных уравнений с медленным изменением переменных без разрывов.

Результаты решения задач интегрирования систем дифференциальных уравнений с использованием функций rkfixed, Rkadapt, Bulstoer формируются системами MathCAD в виде (т+ 1) * (п + 1) - матрицы (таблицы), первый столбец которой содержит значения аргументов от x1 до x2, а остальные п ее столбцов образуются значениями элементов вектора у переменных состояний, исследуемой системы. Таким образом, число элементов каждого из столбцов результирующей матрицы определяется параметром т, введенным в качестве аргумента соответствующей функции.

На рисунке 3 приведены результаты решения одиночного нелинейного дифференциального уравнения с применением rkfixed. Для сопоставления, в представленном файле приводится решение этого же уравнения с использованием функции odesolve.

По условиям применения рассматриваемой группы однотипных функций начальное состояние и искомые переменные должны быть векторами, поэтому для скалярной задачи (рисунок 3) приходится прибегать к искусственному приему. Он заключается в указании системе, что начальное состояние ее первый (ORIGIN = 0) элемент некоего вектора r, а в правых частях уравнения также используется первый элемент некоего вектора результирующих переменных у. Результирующая матрица s в данном случае имеет всего два столбца, первый из которых есть значения аргумента t, а второй значения искомой переменной у. Результаты представлены в виде графика y(t). Привлечение искусственного приема вызвано тем, что рассматриваема задача не является типичной для области рационального использования функций rkfixed, Rkadapt, Bulstoer. Применение odesolve в той же задаче (рисунок 3) иллюстрирует простоту практического использования ВБДУ при решении задач области приложения функции odesolve. Результаты ВБДУ формируются в виде вектора у (рисунок 3) и практически совпадают с результатами, полученными использованием rkfixed.

Рисунок 3. Решение ДУ первого порядка с применением rkfixed и odesolve

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

Рисунок 4. Пример решения системы нелинейных уравнений

Результирующая матрица решений, обозначенная на рисунке 4 через S, имеет 1001 строку и 4 столбца. Вектор искомых переменных имеет три составляющих, первые две из которых выведены на график.

Достаточно большой класс прикладных задач, возникающих, например, при анализе переходных процессов линейных динамических систем, допускает векторно-матричную запись исходных ДУ вида dx/dt = Ах + Вu, z = Нх, где x - n-мерный вектор состояний, u - вектор входных воздействий, z- вектор выходных (наблюдаемых переменных), А, В, Н - матрицы согласованных размерностей. Пример решения подобной системы линейных дифференциальных уравнений, представленной в компактной форме, приведен на рисунке 5. Эта задача соответствует моделированию реакции z(t) динамической системы третьего порядка на единичный скачок. Исходная система уравнений 'эквивалентна дифференциальному уравнению третьего порядка, решение которого с помощью функции odesolve был приведен на рисунке 1. В примере (рисунок 5) применяются функции Bulstoer и rkfixed. Начальный вектор состояний задается нулевым; для этого последнему элементу вектора х присваивается нулевое значение. Выходная переменная z (в данном случае первый элемент вектора х) выводится с помощью вспомогательной матрицы H1, которая учитывает, что первый столбец результирующей матрицы у, как отмечалось, содержит значения аргумента (в нашем случае - времени). На рисунке 5 построен график зависимости выходной переменной от времени. Шаг интегрирования при выбранных Параметрах равен 0.5 и хорошо согласуется с динамическими свойствами моделируемой системы.

Рисунок 5. Пример решения системы линейных ДУ

Для уточнения состава и структуры результирующей матрицы, дополнительно выведены результаты для первых 4 шагов на интервале времени (0..2) с тем же интервалом дискретности, равным 0.5. Результирующая матрица yl получена с помощью rkfixed и выведена в форме матрицы. В первом ее столбце, в соответствии с порядком перечислений аргументов в функции правых частей D, расположены значения времени t с шагом 0.5: Остальные столбцы матрицы yl содержат значения элементов вектора состояний. Выходная переменная zl (рисунок 5) также получена с привлечением вспомогательной матрицы H1. Введение вектора v (рисунок 5) иллюстрирует возможность использования рассматриваемых функций в составе выражений. В данном случае вектор v получен с помощью функции Bulstoer и представляет выходную переменную. Сопоставление данных (рисунки 1 и 5), полученных для одной и той же задачи с применением различных функций, демонстрирует полную идентичность результатов, характерную при корректном выборе интервалов дискретности.

Функции Stiffb, Stiffr

Для интегрирования "жестких" дифференциальных уравнений и систем таких уравнений в MathCAD предусмотрены встроенные функции:

Stiffb(y0,xl,x2,m,D,J) - решение уравнений методом Булирш-Штоера;

Stiffr(y0,xl,x2,m,D,J) - решение уравнений методом Розенброка,

Их применение совершенно аналогично. По сравнению с ранее рассмотренными функциями, эта группа функций имеет дополнительный аргумент J - [n х (n + 1)]-матричную функцию, первый столбец которой равен вектору производных D(x,y) по х, а остальную часть составляет матрица Якоби для вектора D(x,у). Формат матричной функции J, который необходимо учитывать при ее вводе, имеет вид J(x,y).

Вывод результата функциями Stiffb, Stiffr производится аналогично функциям rkfixed, Rkadapt, Bulstoer, т.е. в виде (m + 1) х (п + 1) - матрицы (таблицы), первый столбец которой содержит значения аргументов от xl до x2, а остальные п ее столбцов образуются значениями элементов вектора у переменных состояний исследуемой системы.

При формировании матрицы J и вычислении якобиана (определителя матрицы Якоби) удобно использовать возможности символьных преобразований.

Методы, реализованные функциями rkfixed, Rkadapt, Bulstoer, практически неприменимы к жестким дифференциальным уравнениям, поскольку их решение требует исключительно малого шага интегрирования. Применение функций Stiffb, Stiffr по сравнению с ранее рассмотренными функциями характеризуется большой точностью интегрирования жестких дифференциальных уравнений при существенно (иногда на порядки) увеличенном шаге интегрирования. Пример использования функции Stiffb для системы дифференциальных уравнений второго порядка приведен на рисунке 6. Для сопоставления, интегрирование исходной системы осуществлено также и с использованием функций rkfixed для двух вариантов шага интегрирования.

Вспомогательная матрица H1 использована для выведения первой составляющей вектора состояний рассматриваемой системы. Значения этой составляющей, полученные с применением указанных функций, образуют векторы: zl - функция rkfixed с шагом 0.05, z2 - функция rkfixed с шагом 0.01, Z0 - функция Stiffb с шагом 0.05. Функция Stiffb (рисунок 6) обеспечила точное решение на интервале (0, 10) с шагом 0.05. При таком же шаге решение, полученное с применением функции rkfixed, имеет заметное искажение. Для достижения сравнимой точности, интегрирования с помощью функции rkfixed рассматриваемой системы потребовало уменьшения шага в пять раз, т.е. число шагов увеличилось в пять раз.

Рисунок 6. Пример применения функции Stiffb

Функций rkadapt, bulstoer, stiffb, stiffr

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

rkadapt(y0,xl,x2,a,D,k,s) - решение обыкновенных уравнений методом Рунге-Кутта с переменным шагом интегрирования;

bulstoer(y0,xl,x2,a,D,k,s) - решение уравнений методом Булирш-Штоера;

stiffb(y0,xl,x2,a,D,J,k,s) - решение "жестких" систем методом Булирш-Штоера;

stiffr(y0,xl,x2,a,D,J,k,s) - решение "жестких" систем методом Розенброка.

Здесь к - число промежуточных точек, в которых функции будут находить решение в интервале (xl, х2). Значение к не должно быть меньше двух; s -минимально допустимая величина шага; а - параметр, управляющий точностью вычислений (обычно а = 0.001). Остальные аргументы описаны мной ранее, при рассмотрении соответствующих функций, начинающихся с прописных букв.

Правила, которые необходимо учитывать при использовании приведенных _/ функций, в целом аналогичны правилам применения ранее рассмотренных одноименных функций. Отличия касаются лишь формы выдачи результата и отсутствия необходимости задавать шаг интегрирования. Рассматриваемые функции формируют результат в виде

[кх (n + 1)] - матрицы (таблицы), первый столбец которой содержит значения аргументов от x1 до x2, а остальные n ее столбцов образуются значениями элементов вектора у переменных состояний исследуемой системы. Значение может быть задано равным двум. При к = 2 результирующая матрица имеет всего две строки - для начального и конечного состояний. В случае задания к > 2, функции сами выбирают местоположение точек в заданном интервале интегрирования.

Применение одной из рассматриваемых функций (bulstoer) для получения финальных значений переменных проиллюстрировано на рисунке 7. В качестве исходной системы была принята система дифференциальных уравнений, приведенная на рисунке 5.

В примере (рисунок 7) задан ненулевой вектор начальных значений. Для выбранного числа точек к = 2, в которых находится решение, результирующая матрица s, как отмечалось, содержит две строки. Первая строка соответствует; начальному состоянию системы, а вторая - ее состоянию на правом конце выбранного диапазона (при х2 = 25). Справа приведены данные теоретического значения установившегося режима выходной переменной и её расчетное значение при х2 = 25, полученное применением функции bulstoer. В этом примере выбранный интервал практически соответствует времени затухания переходных процессов системы, поэтому полученное расчетное значение незначительно отличается от теоретического.

Рисунок 7. Анализ финального участка переходного процесса

Применение остальных функций этой группы аналогично применению одноименных функций, которые имеют прописной начальный символ.

Функции sbval, bvalfit

Встроенные функции sbval и bvalfit предназначены для решения (совместно с другими функциями) двухточечных краевых задач. В задачах этого класса дифференциальные уравнения могут иметь ограничения в начале, в середине и конце интервала интегрирования. Часть этих ограничений обычно неизвестно. Функции sbval и bvalfit позволяют определить недостающие граничные условия, которые используются затем для получения решения на всем интервале интегрирования с помощью другой (например, rkfixed) функции. При решении задач доопределения граничных условий функции sbval и bvalfit реализуют так называемый алгоритм, пристрелки (shooting method), обеспечивающий подбор недостающих граничных условий с опорой на известную часть ограничений. В наиболее простом случае (sbval)-этим алгоритмом осуществляется серия решений ДУ с различными начальными условиями, имеющая цель определить такие начальные значения, при которых решение выходит в заданную на правом конце точку.

Исходные условия для двухточечных краевых задач могут быть различными. В наиболее распространенном случае для ДУ задаются граничные условия на левом и правом концах интервала интегрирования. Определение недостающих граничных, условий в таких задачах осуществляется с помощью функции sbval. При задании граничных условиях не только на концах интервала, но и внутри его, используется функция bvalfit.

Функция sbval(v,xl,x2,D,L,s) имеет следующие аргументы:

v - r-мерный вектор начальных приближений искомых недостающих граничных условий;

xl, х2 - левая и правая границы интервала интегрирования (область определения аргумента х);

D - n-мерная векторная функция двух аргументов (х, Y), где Y - n-мерный вектор переменных ДУ; D содержит первые производные вектора Y (правые части формы Коши). Формат этой векторной функции (ее изображение в документах MathCAD) имеет вид: D (х, Y);

L - n-мерная векторная функция двух аргументов (xl, v). Это функция начальных условий в точке xl; вектор L содержит (n - r)известных значений (числа) и г неизвестных начальных значений (символы соответствующих элементов вектора v). Формат: L (xl, v);

s - r-мерная векторная функция двух аргументов (х2, Y), содержащая разности между известными граничными условиями и получаемыми решениями на правом конце интервала интегрирования (в точке х2). Формат: s (x2,Y).

Результат вычислений с помощью функции sbval формируется MathCAD в виде r-мерного вектора-столбца недостающих граничных условий.

Пример использования функции sbval для решения двухточечной краевой задачи приведен на рисунке 8. При составлении примера использованы дифференциальное уравнение четвертого порядка ранее рассмотренной задачи Коши (рисунок 2). Краевая задача содержит это ДУ, два левых (у"(0) = - 0.1, у"'(0) = 0) и два правых (у(6) = 50, у'(6) = 50) граничных условия. Таким образом, для решения задачи Коши недостает двух левых граничных условий, которые определяются функцией sbval из условия "попадания" в заданные правые условия;

Рисунок 8. Пример применения функции sbval

Условия применения функции sbval требуют, чтобы ДУ высокого порядка было предварительно представлено в виде эквивалентной системы ДУ первого; порядка (стандартная форма или форма Коши). Для рассматриваемого дифференциального уравнения у""(х)+у²(х) = 2.5х² эта система имеет вид:

Вектор состояний Y полученной системы ДУ имеет четыре составляющих, обозначенные (рисунок 8.) как элементы вектора (ORIGIN = 0), причем переменная y(x) исходного уравнения соответствует первой составляющей ( ) вектора Y. Правые части этих уравнений нашли свое отражение в векторной функции D(x,Y). Вектор v (рисунок 8) содержит начальные приближения недостающих граничных условий. Из представленного файла видно, что вектор L играет роль начального вектора состояний; два последних элемента этого вектора заданы и, потому, известны, а два первых элемента неизвестны ивходят в состав L как элементы вектора v. Вектор-функция s содержит целевую установку решаемой задачи - невязки решения на правом конце; минимизацией этих невязок и определяются недостающие начальные значения. После задания всех аргументов функции sbval, ее использование для получения вектора-столбца недостающих граничных условий не составляет труда. Из представленного решения (рисунок 8) ясно, что вектор начальных условий, рассматриваемой системы, который обеспечивает "выход" решения на заданные справа значения переменной у(х) и ее производной при х = 6, равен: | -0.042, 1.874, -0.1, 0| .

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

На рисунке 9 приведены результаты решения задачи Коши для рассмотренного (рисунок8) примера. Интегрирование осуществляется с применением функции rkfixed; шаг интегрирования принят равным 0.1.

Рисунок 9.Решение задачи при найденных начальных условиях

Фрагмент (рисунок 9) представляет собой продолжение файла, приведенного ранее (рисунок 8) поэтому два первых элемента в начальном векторе состояний Y обозначены как элементы вектора V найденных недостающих граничных условий. С помощью вспомогательных матриц H1 и Н2 из результирующей матрицы yl формируются векторы z1 и z2 значений двух первых элементов вектора состояний рассматриваемой системы, для которых условиями задачи были заданы правые граничные значения, равные 50. Полученные в результате интегрирования финальные значения этих элементов выведены в нижней части рисунка 9 и представлены в составе шестидесятой строки результирующей матрицы. График иллюстрирует характер процесса перехода двух первых элементов вектора состояний от начальных значений к заданным конечным. Достигнутая точность выхода переходных процессов на заданные правые граничные условия свидетельствует о высоком качестве решения двухточечной задачи.

Функция bvalfit(vl,v2,xl,x2,xf,D,Ll,L2,sf) имеет следующие аргументы:

vl, v2 - векторы начальных приближений недостающих левых и правых граничных условий;

xl, х2- левая и правая границы интервала интегрирования (область определения аргумента х);

xf- промежуточная точка внутри интервала интегрирования,

D - n-мерная векторная функция двух аргументов (х, Y), где Y - n-мерный вектор переменных ДУ; D содержит первые производные вектора Y (правые части формы Коши). Формат этой векторной функции (ее изображение в документах MathCAD) имеет вид: D (x,Y);

L1 - n-мерная векторная функция двух аргументов (xl, vl). Это функция начальных условий в точке xl; вектор L1 содержит как заданные значения (числа), так и неизвестные начальные значения (символы соответствующих элементов вектора vl). Формат: L (xl, vl);

L2 - n-мерная векторная функция двух аргументов (х2, v2). Это функция конечных условий в точке х2; вектор L2 содержит как заданные значения (числа), так и неизвестные конечные значения (символы соответствующих элементов вектора v2). Формат: L (х2, v2);

sf - векторная функция двух аргументов (xf, Y), содержащая значение вектора Y в промежуточной точке xf. Формат: sf (xf, Y).

Результат вычислений с помощью функции bvalfit формируется MathCAD в виде вектора-строки недостающих граничных условий.

Дополнительное условие Y(xf) в промежуточной точке используется функцией bvalfit для определения недостающих граничных условий с помощью алгоритма двусторонней пристрелки, имеющего цель обеспечить выход пременных на заданные промежуточные значения. Чаще всего условие Y(xf) вводится в задачу, когда первые производные вектора Y имеют разрыв в точке xf. Для сопряжения решений слева и справа от xf, векторная функция sf задается в виде: sf (xf, Y) := Y. В центре ресурсов и справочной базе MathCAD имеются примеры