Матричные вычисления в Mathcad
Что такое жесткие ОДУ?
9.4.1. Что такое жесткие ОДУ?Исторически интерес к жестким системам возник в середине XX в. при изучении уравнений химической кинетики с одновременным присутствием очень медленно и очень быстро протекающих химических реакций (оно будет рассмотрено в разд. 9.4.3). Тогда неожиданно оказалось, что считавшиеся исключительно надежными методы Рунге—Кутты стали давать сбой при расчете этих задач.
Строгого общепринятого математического определения жестких ОДУ нет. В рамках этой книги будем считать, что жесткие системы — это те уравнения, решение которых получить намного проще с помощью определенных неявных методов, чем с помощью явных методов (типа тех, что были рассмотрены в предыдущих разделах). Примерно такое определение было предложено в 1950-х гг. классиками в этой области Кертиссом и Хиршфельдером. Начнем разговор о жестких ОДУ с примера нежесткого уравнения (листинг 9.8), решение которого удается получить численным методом Рунге— Кутты (Рисунок 9.12).
Модели динамики биологических популяций
9.5.1. Модели динамики биологических популяцийМодель взаимодействия "хищник—жертва" независимо предложили в 1925— 1927 гг. Лотка и Вольтерра. Два дифференциальных уравнения (листинг 9.13) моделируют временную динамику численности двух биологических популяций жертвы уо и хищника уь Предполагается, что жертвы размножаются с постоянной скоростью с, а их численность убывает вследствие поедания хищниками. Хищники же размножаются со скоростью, пропорциональной количеству пищи (с коэффициентом г), и умирают естественным образом (смертность определяется константой о). В листинге рассчитываются три решения о, с, р для разных начальных условий.
Встроенные функции для решения систем ОДУ
9.3.1. Встроенные функции для решения систем ОДУВ Mathcad имеется несколько встроенных функций, которые позволяют решать задачу Коши различными численными методами. Для "хороших" нежестких систем ОДУ применяются следующие функции:
ВНИМАНИЕ!
Соблюдайте регистр первой буквы рассматриваемых функций, поскольку это влияет на выбор алгоритма счета, в отличие от многих других встроенных функций Mathcad, например, Find=f ind (см. разд. 9.3.2).
Каждая из приведенных функций выдает решение в виде матрицы размера (M+1)х(N+1). В ее левом столбце находятся значения аргумента t, делящие интервал на равномерные шаги, а в остальных N столбцах — значения искомых функций y0(t) ,y1(t), ... ,yN-i(t), рассчитанные для этих значений аргумента (Рисунок 9.5). Поскольку всего точек (помимо начальной) м, то строк в матрице решения будет всего M+1.
В листинге 9.3 приведен пример решения той же самой системы ОДУ осциллятора с затуханием (см. разд. 9.1 и 9.2) при помощи первой из функций rkfixed. Результат расчетов представлен на Рисунок 9.5 как содержимое матрицы и на Рисунок 9.6 в виде графика. Точки, в которых получено решение, отмечены на графике Рисунок 9.6 кружками. Чтобы использовать другой численный алгоритм, достаточно поменять имя функции rkfixed в последней строке листинга на другое (на практике как раз более эффективны функции Rkadapt и Bulstoer).
Задачи Коши для ОДУ
9.1.1. Задачи Коши для ОДУДифференциальные уравнения — это уравнения, в которых неизвестными являются не переменные (т. е. числа), а функции одной или нескольких переменных. Эти уравнения (или системы) включают соотношения между искомыми функциями и их производными. Если в уравнения входят производные только по одной переменной, то они называются обыкновенными дифференциальными. В противном случае говорят об уравнениях в частных производных (см. главу II). Таким образом, решить (иногда говорят проинтегрировать) дифференциальное уравнение — значит, определить неизвестную функцию на определенном интервале изменения ее переменных.
ОДУ с неизвестной функцией у (t), в которое входят производные этой функции вплоть до y(N) (t), называется ОДУN-го порядка. В частности, уравнение первого порядка может по определению содержать помимо самой искомой функции y(t) только ее первую производную у'(t), второго порядка— у' (t) и у'' (t) и т. д. В подавляющем большинстве практических случаев дифференциальное уравнение можно записать в стандартной форме (форме Коши): у'(t)=f (у (t) ,t). Уравнение второго порядка может содержать, помимо самой функции, ее первую и вторую производные и т. д.
Автоколебания
9.5.2. АвтоколебанияРассмотрим решение уравнения Ван дер Поля, описывающего электрические колебания в замкнутом контуре, состоящем из соединенных последовательно конденсатора, индуктивности, нелинейного сопротивления и элементов, обеспечивающих подкачку энергии извне (листинг 9.14). Неизвестная функция времени y(t) имеет смысл электрического тока, а в параметре ц заложены количественные соотношения между составляющими электрической цепи, в том числе и нелинейной компонентой сопротивления.
Фазовый портрет динамической системы
9.1.2. Фазовый портрет динамической системыМодели, основанные на задачах Коши для ОДУ, часто называют динамическими системами, подчеркивая, что, как правило, они содержат производные по времени t и описывают динамику некоторых параметров. Проблемы, связанные с динамическими системами, на самом деле весьма разнообразны и зачастую не сводятся к простому интегрированию ОДУ. Некоторые из них мы обозначим в данном разделе, отметив, что для изучения динамических систем центральным моментом является анализ фазовых портретов, т. е. решений, получающихся при выборе всевозможных начальных условий.
Решение ОДУ часто удобнее изображать не в виде графика у0 (t), y1(t), ..., как на Рисунок 9.1, а в фазовом пространстве, по каждой из осей которого откладываются значения каждой из найденных функций. При таком построении графика аргумент t будет присутствовать на нем лишь параметрически. В рассматриваемом случае двух ОДУ (мы свели к ним в предыдущем разделе дифференциальное уравнение осциллятора второго порядка) фазовое пространство является координатной плоскостью, а решение представляет собой кривую, или, по-другому, траекторию, выходящую из точки, координаты которой равны начальным условиям (Рисунок 9.2). В общем случае, если система состоит из N ОДУ, то фазовое пространство является N-мерным. При N>3 наглядность теряется, и для визуализации фазового пространства приходится строить его различные проекции или прибегать к другим специальным приемам (например, отображению Пуанкаре).
Функции для решения жестких ОДУ
9.4.2. Функции для решения жестких ОДУРешение жестких систем дифференциальных уравнений можно осуществить только с помощью встроенных функций, аналогичных по действию семейству рассмотренных выше функций для обычных ОДУ:
Как вы можете заметить, для двух последних функций серьезным отличием от функций, решающих нежесткие системы, является добавление к стандартному набору параметров дополнительной матричной функции, задающей якобиан системы ОДУ. Решение выдается в виде матрицы, по форме идентичной аналогичным функциям решения нежестких задач Коши.
Примечание 1
Примечание 1
Встроенная функция Radau, которая не требует явного задания якобиана системы уравнений, появилась в версии Mathcad 2001I, а остальные две — в Mathcad 2001.
Решение жесткой задачи из предыдущего раздела при помощи функции Radau приведено в листинге 9.9. Результат показан в виде графика на Рисунок 9.14 вместе с графиком решения менее жесткой задачи (для которого применялся листинг 9.8). Как вы видите, хватило всего пяти точек разбиения интервала интегрирования жесткого ОДУ, чтобы метод с ним справился. Специфика применения других встроенных функций, требующих дополнительного задания якобиана, будет рассмотрена в следующем разделе на примере уравнения химической кинетики.
Решение одного уравнения (N=1)
9.3.2. Решение одного уравнения (N=1)Метод решения ОДУ при помощи встроенных функций rkfixed, Rkadapt или Bulstoer (в противоположность вычислительному блоку Given/ odesoive) сохранился с прежних версий Mathcad (до 2000-й). В большинстве случаев лучше использовать вычислительный блок Given/odesolve, который выигрывает в простоте и в наглядности, однако иногда предпочтительнее решать ОДУ первого порядка с помощью второго способа, например, при следующих обстоятельствах:
Поскольку решение вторым способом одного ОДУ не отличается от решения систем ОДУ (см. предыдущий разд.), приведем пример его использования (листинг 9.4) практически без комментариев. Отметим лишь в случае одного ОДУ, что как само уравнение, так и начальное условие можно задавать не в векторной, а в скалярной форме. Результат выдается в виде матрицы размерности мх2, которая состоит из двух столбцов: в одном находятся значения аргумента t (от t0 до t1 включительно), а в другом соответствующие значения искомой функции у (t).
Пример химическая кинетика
9.4.3. Пример: химическая кинетикаРассмотрим классическую модель химической кинетики (Робертсон, 1966), которая как нельзя лучше передает смысл понятия жесткости ОДУ.
Попытка решения стандартными методами
Рассмотрим составную схему химического взаимодействия трех веществ. Пусть вещество "О" медленно превращается в "1": "0"->"1" (со скоростью 0.1), вещество "1" при каталитическом воздействии самого себя превращается очень быстро в вещество "2": "1"+"1"->"2"+"1" (103). И, наконец, подобным образом (но со средней скоростью) реагируют вещества "2" и "1": "1"+"2"->"0" + "2" (102). Система ОДУ, описывающая динамику концентрации реагентов, с попыткой решения методом Рунге—Кутты, приведена, в символической форме в первой строке листинга 9.10.
Решение систем ОДУ в одной заданной точке
9.3.3. Решение систем ОДУ в одной заданной точкеЗачастую при решении дифференциальных уравнений требуется определить значения искомых функций не на всем интервале (t0.t1), а только в одной его последней точке. Например, весьма распространены задачи поиска аттракторов динамических систем. Известно, что для широкого класса ОДУ одна и та же система при разных (или даже любых, как рассмотренный выше пример осциллятора с затуханием) начальных условиях при t->? приходит в одну и ту же точку (аттрактор). Поэтому часто нужно определить именно эту точку.
Такая задача требует меньше ресурсов компьютера, чем решение системы ОДУ на всем интервале, поэтому в Mathcad имеются модификации встроенных функций Rkadapt и Buistoer. Они имеют несколько другой набор параметров и работают значительно быстрее своих аналогов:
Как легко заметить, вместо числа шагов на интервале интегрирования ОДУ в этих функциях необходимо задать точность расчета численным методом значения функций в последней точке. В этом смысле параметр асе похож на константу TOL, которая влияет на большинство встроенных численных алгоритмов Mathcad. Количество шагов и их расположение определяются численным методом автоматически, чтобы обеспечить эту точность. Два последних параметра нужны для того, чтобы пользователь мог искусственно повлиять на разбиение интервала на шаги. Параметр к служит для того, чтобы шагов не было чрезмерно много, причем нельзя сделать k>1000. Параметр s — для того, чтобы ни один шаг не был слишком малым для появления больших погрешностей при разностной аппроксимации дифференциальных уравнений внутри алгоритма. Эти параметры следует задавать явно, исходя из свойств конкретной системы ОДУ. Как правило, проведя ряд тестовых расчетов, можно подобрать их оптимальный набор для каждого конкретного случая.
Пример использования функции buistoer для той же модели линейного осциллятора приведен в листинге 9.5. В его первых двух строках, как обычно, определяется система уравнений и начальные условия; в следующей строке матрице и присваивается решение, полученное с помощью buistoer. Структура этой матрицы в точности такая же, как и в случае решения системы ОДУ посредством уже рассмотренных нами ранее встроенных функций (см. разд. 9.3.1). Однако в данном случае нам интересна только последняя точка интервала. Поскольку сделанное численным методом количество шагов, т. е. размер матрицы и, заранее неизвестно, то его необходимо предварительно определить. Это сделано в четвертой строке листинга, присваивающей это число переменной м, в этой же строке оно выведено на экран. В предпоследней строке листинга осуществлен вывод решения системы ОДУ на конце интервала, т. е. в точке t=s0 в виде вектора. В последней строке для примера еще раз выводится искомое значение первой функции из системы ОДУ (сравните его с соответствующим местом вектора из предыдущей строки). Чтобы попробовать альтернативный численный метод, достаточно в листинге 9.5 заменить имя функции buistoer на rkadapt.
Странный аттрактор
9.5.3. Странный аттракторОдна из самых знаменитых динамических систем предложена в 1963 г. Лоренцем в качестве упрощенной модели конвективных турбулентных движений жидкости в нагреваемом сосуде тороидальной формы. Система состоит из трех ОДУ и имеет три параметра модели (листинг 9.15). Поскольку неизвестных функций три, то фазовый портрет системы должен определяться не на плоскости, а в трехмерном пространстве.
Брюсселятор
9.5.4. БрюсселяторДо сих пор в этой главе в качестве примеров расчета динамических систем мы приводили графики 1—3 траекторий на фазовой плоскости. Однако для надежного исследования фазового портрета необходимо решить систему ОДУ большое количество раз с самыми разными начальными условиями (и, возможно, с разным набором параметров модели), чтобы посмотреть, к каким аттракторам сходятся различные траектории. В Mathcad можно реализовать эту задачу, например, в форме алгоритма, приведенного в листинre 9.16 для решения системы уравнений автокаталитической химической реакции с диффузией. Эта модель, называемая моделью брюсселятора, предложена в 1968 г. Лефевром и Пригожиным. Неизвестные функции отражают динамику концентрации промежуточных продуктов некоторой реальной химической реакции. Параметр модели в равен исходной концентрации катализатора.
О численных методах
9.3.4. О численных методахВ завершение раздела сделаем несколько важных замечаний относительно выбора численного алгоритма решения ОДУ и задания его параметров. Они не претендуют на общность, но, надеемся, будут весьма полезны читателю, особенно в случае возникновения проблем.
Рекомендации по выбору численного алгоритма
Все численные методы решения ОДУ основаны на аппроксимации дифференциальных уравнений разностными аналогами. В зависимости от конкретной формы аппроксимации получаются алгоритмы различной точности и быстродействия. В Mathcad использован наиболее популярный алгоритм Рунге—Кутты четвертого порядка, описанный в большинстве книг по методам вычислений. Он обеспечивает малую погрешность для широкого класса систем ОДУ за исключением жестких систем. Поэтому в большинстве случаев стоит применять функцию rkfixed. Если по различным причинам время расчетов становится критичным или точность неудовлетворительна, стоит попробовать вместо rkfixed другие функции, тем более, что сделать это очень просто благодаря одинаковому набору параметров. Для этого нужно только поменять имя функции в программе.
Функция Rkadapt может быть полезна в случае, когда известно, что решение на рассматриваемом интервале меняется слабо либо существуют участки медленных и быстрых его изменений. Метод Рунге—Кутты с переменным шагом разбивает интервал не на равномерные шаги, а более оптимальным способом. Там, где решение меняется слабо, шаги выбираются более редкими, а в областях его сильных изменений — частыми. В результате для достижения одинаковой точности требуется меньшее число шагов, чем для rkfixed. Метод Булирша— Штера Buistoer часто оказывается более эффективным для поиска гладких решений.
Имея в виду сделанные замечания, приведем краткую сводку алгоритмов решения задач Коши для ОДУ, отмечая, какие из встроенных функций следует использовать в конкретных случаях.
Мы не будем в этой книге приводить описание численных алгоритмов, поскольку им посвящено огромное количество книг, и читатель, несомненно, найдет соответствующую литературу с подробным объяснением их работы. Тем не менее, разберем несколько важных моментов, касающихся выбора алгоритма и его параметров.
Число шагов
Все встроенные функции требуют задания для численного метода количества шагов, разбивающих интервал интегрирования ОДУ. Это очень важный параметр, непосредственно влияющий на погрешность результата и скорость расчетов. При решении систем ОДУ многие проблемы могут быть устранены при помощи простой попытки увеличить число шагов численного метода. В частности, сделайте так при неожиданном возникновении ошибки "Found a number with a magnitude greater than 10^307" (Найдено число, превышающее значение 10307). Данная ошибка не обязательно говорит о том, что решение в действительности расходится, а возникает из-за недостатка шагов для корректной работы численного алгоритма.
В качестве иллюстрации приведем результаты простых расчетов (листинг 9.6), в которых определяется пользовательская функция и(м>, представляющая решение ОДУ в зависимости от числа шагов м. На Рисунок 9.8 показано несколько графиков решения системы ОДУ затухающего линейного осциллятора для нескольких значений м.
Аттрактор Лоренца на фазовой плоскости (продолжение листинга 9 15)
Рисунок 9.22. Аттрактор Лоренца на фазовой плоскости (продолжение листинга 9.15)
Дифференциальное уравнение Nго порядка
9.2. Дифференциальное уравнение N-го порядкаДля решения ОДУ порядка N>1 в Mathcad предусмотрены две возможности:
В данном разделе рассматривается первый из вариантов, а второй будет описан в следующем, применительно к общему случаю системы N уравнений. Заметим лишь, что первый путь предпочтительнее из соображений наглядности представления задачи и результатов, а второй дает пользователю больше рычагов воздействия на параметры численного метода.
Вычислительный блок для решения ОДУ, реализующий численный метод Рунге—Кутты, состоит из трех частей:
Примечание 1
Примечание 1
На запись ОДУ в пределах вычислительного блока накладывается несколько ограничений. Во-первых, ОДУ должно быть линейно относительно старшей производной, т. е. фактически должно быть поставлено в стандартной форме. Во-вторых, начальные условия должны иметь форму у (t0) =b или y(N) (t0) =b, а не более сложную (как, например, встречающаяся в некоторых математических приложениях форма у (t0)+y' (t0)=b).
Примечание 2
Примечание 2
Допустимо, и даже часто предпочтительнее, задание функции Odesoive (t,t1, step) с тремя параметрами, где step— внутренний параметр численного метода, определяющий количество шагов в которых метод Рунге—Кутты будет рассчитывать решение дифференциального уравнения. Чем больше step, тем с лучшей точностью будет получен результат, но тем больше времени будет затрачено на его поиск. Помните, что подбором этого параметра можно заметно (в несколько раз) ускорить расчеты без существенного ухудшения их точности.
Пример решения задачи Коши для ОДУ второго порядка, являющегося усложнением модели из листинга 9.1 на нелинейный случай (с параметром квадратичной нелинейности у)> приведен в листинге 9.2, а результат — на Рисунок 9.4. Учтите, что символ производной допускается вводить как средствами панели Calculus (Вычисления), как это сделано в листинге 9.1, так и в виде штриха, набрав его с помощью сочетания клавиш
Фазовый портрет брюсселятора при В=0 5 (продолжение листинга 9 16)
Рисунок 9.24. Фазовый портрет брюсселятора при В=0.5 (продолжение листинга 9.16)
Как видно из Рисунок 9.24, все траектории, вышедшие из разных точек, асимптотически стремятся к одному и тому же аттрактору (1,0.5). Из теории динамических систем нам известно, что такой аттрактор называется узлом (с узлом мы уже встречались в примере разд. 9.4.1). Конечно, в общем случае при анализе фазового портрета желательно "прощупать" большее число траекторий, задавая более широкий диапазон начальных условий. Не исключено, что в других областях фазовой плоскости траектории будут сходиться к другим аттракторам.
Эволюцию фазового портрета брюсселятора можно наблюдать, проводя расчеты с различным параметром в. При его увеличении узел будет сначала постепенно смещаться в точку с координатами (1, в), пока не достигнет бифуркационного значения B=2. В этой точке происходит качественная перестройка портрета, выражающаяся в рождении предельного цикла. При дальнейшем увеличении в происходит лишь количественное изменение параметров этого цикла. Решение, полученное при B=2.5, показано на Рисунок 9.25.
Примечание 1
Примечание 1
Чтобы найти аттракторы динамической системы, как известно, нужно решить систему алгебраических уравнений, получающуюся из системы ОДУ заменой нулями их левых частей. Эти задачи также удобно решать средствами Mathcad (см. главу 5). В частности, исследование зависимости фазового портрета от параметров системы ОДУ и поиск бифуркаций можно проводить методами продолжения (с/и. разд. 5.3.3).
с расчетом динамических систем, несомненно,
Рисунок 9.25. Фазовый портрет брюсселятора при В=2.5
Читатели, сталкивающиеся с расчетом динамических систем, несомненно, оценят возможности Mathcad по построению фазовых портретов и исследованию бифуркаций. Возможно также, что они найдут лучшие программные решения этой задачи, чем алгоритм, предложенный в данном разделе автором.
Фазовый портрет полученный bulstoer (a) и rkadapt (б) (продолжение листинга 9 5)
Рисунок 9.10. Фазовый портрет, полученный bulstoer (a) и rkadapt (б) (продолжение листинга 9.5)
Поэтому для экономии времени расчетов в функции buistoer можно выбирать и большие асе.
Чтобы обеспечить заданную точность, данные алгоритмы могут изменять как количество шагов, разбивающих интервал (t0,t1), так и их расположение вдоль интервала. Чтобы выяснить, на сколько шагов разбивался интервал при расчетах, воспользуемся простой программой, представленной в листинге 9.7. В ней определены пользовательские функции R(е) и B(е), вычисляющие для конкретной задачи (методов Рунге—Кутты и Булирша—Штера соответственно) зависимость числа шагов (т. е. числа строк получающейся матрицы) от параметра е=асс. Графики функций R(e) и B(е) показаны на Рисунок 9.11. Как видно, методу Рунге—Кутты при увеличении требуемой точности требуется все возрастающее количество шагов (и, соответственно, времени на расчеты), в противоположность методу Булирша—Штера, демонстрирующему большую экономичность.
Примечание 1
Примечание 1
Максимальное количество шагов алгоритма ограничивается еще одним параметром k встроенных функций buistoer и rkadapt (с/и. рази 9.3.3).
Фазовый портрет решения системы ОДУ при М=40 и М=200 (продолжение листинга 9 6)
Рисунок 9.9. Фазовый портрет решения системы ОДУ при М=40 и М=200 (продолжение листинга 9.6)
Как можно заметить, несмотря на, в целом, правильное решение ОДУ в узлах (оно показано кружками и квадратами), профиль решения при малых м получается неверным. Еще лучше это видно на Рисунок 9.9, представляющем фазовый портрет системы с разным числом шагов. Видно, что для м=200 решение получается более гладким (конечно, при этом увеличивается и время расчетов).
Обыкновенные дифференциальные уравнения динамические системы
Глава 9. Обыкновенные дифференциальные уравнения: динамические системы| Содержание |
График решения системы ОДУ осциллятора (продолжение листинга 9 3)
Рисунок 9.6. График решения системы ОДУ осциллятора (продолжение листинга 9.3)
ВНИМАНИЕ!
Обратите внимание на некоторое разночтение в обозначении индексов вектора начальных условий и матрицы решения. В ее первом столбце собраны значения нулевой компоненты искомого вектора, во втором столбце — первой компоненты и т. д.
График решения (слева) и фазовый портрет (справа) модели конкуренции популяций
Рисунок 9.18. График решения (слева) и фазовый портрет (справа) модели конкуренции популяций
График решения (слева) и фазовый портрет (справа) уравнения Ван дер Поля (продолжение листинга 9 14)
Рисунок 9.19. График решения (слева) и фазовый портрет (справа) уравнения Ван дер Поля (продолжение листинга 9.14)
Решением уравнения Ван дер Поля являются колебания, вид которых для µ=1 показан на Рисунок 9.19. Они называются автоколебаниями и принципиально отличаются от рассмотренных нами ранее (например, колебаний маятника в модели осциллятора или численности популяций в модели Вольтерpa) тем, что их характеристики (амплитуда, частота, спектр) не зависят от начальных условий, а определяются исключительно свойствами самой динамической системы. Через некоторое время расчетов после выхода из начальной точки решение выходит на один и тот же цикл колебаний, называемый предельным циклом. Аттрактор типа предельного цикла является замкнутой кривой на фазовой плоскости. К нему асимптотически притягиваются все окрестные траектории, выходящие из различных начальных точек, как изнутри (Рисунок 9.19), так и снаружи (Рисунок 9.20) предельного цикла.
Примечание 1
Примечание 1
Если компьютер у вас не самый производительный, то расчет в Mathcad фазового портрета с Рисунок 9.19, 9.20 может занять относительно продолжительное время, что связано с численным определением сначала решения y(t), а потом его производной. Время расчетов можно было бы существенно сократить, если использовать вместо вычислительного блока Given/odesolve одну из встроенных функций, которые выдают решение в виде матрицы, например, rkfixed.
Примечание 2
Примечание 2
По мере возрастания параметра µ модель Ван дер Поля становится все более жесткой. Например, при ц=5000 решение уже придется искать при помощи специфических методов решения жестких задач (см. разд. 9.3). Это еще раз доказывает, что одна и та же система ОДУ с различными коэффициентами может быть жесткой в разной степени.
График решения (слева) и фазовый
Рисунок 9.17. График решения (слева) и фазовый портрет (справа) системы "хищник—жертва" (продолжение листинга 9.13)
Рассмотренную модель динамики двух популяций легко можно модифицировать, изменив тип взаимодействия "хищник—жертва" на тип конкуренции. Для этого надо учесть, что рост численности каждой популяции тормозит, во-первых, межвидовая, и, во-вторых, внутривидовая конкуренция.
В результате система (во второй строке листинга) запишется в виде:
где матрица г задает коэффициенты убывания численности вследствие конкурентной борьбы (диагональные элементы соответствуют внутри-, а недиагональные — межвидовой конкуренции).
График решения (для разных начальных условий) и фазовый портрет для описанной системы ОДУ показаны на Рисунок 9.18. Как видно, конкурентная борьба приводит к установлению некоторого стационарного состояния, выражающего равновесие видов. Особая точка, к которой стремится решение системы ОДУ подобным образом, называется узлом.
Примечание 2
Примечание 2
Соответствующий файл с программой вы найдете на компакт-диске, вместе с более простой моделью динамики одной популяции с внутривидовой конкуренцией (называемой логистической моделью).
демонстрирует решение
Листинг 9.1 демонстрирует решение простого ОДУ второго порядка, описывающего модель затухающего гармонического осциллятора. Само уравнение приведено во второй строке листинга, после задания параметров модели, а вычисленный результат, т. е. искомая функция у (t), показан на Рисунок 9.1.Примечание 1
Примечание 1
Модель гармонического осциллятора описывает, в частности, колебания маятника: y(t) описывает изменения угла его отклонения от вертикали, у' (t) — угловую скорость маятника, у" (t) — ускорение, а начальные условия, соответственно, начальное отклонение маятника у (0) =1.0 и начальную скорость у' (0)=0. Важно отметить, что модель является линейной, т.е. неизвестная функция (и ее производные) входят в уравнение в первой степени.
Решение задачи Коши
Листинг 9.1. Решение задачи Коши для ОДУ второго порядка (модель затухающего гармонического осциллятора)
Как показывает листинг 9.1, помимо самого уравнения потребовалось определить два начальных условия (третья и четвертая строки листинга) — начальные значения y(t) и у' (t) при t=0. Вообще говоря, ОДУ (или система ОДУ) имеет единственное решение, если помимо уравнения определенным образом заданы начальные или граничные условия. В соответствующих курсах высшей математики доказываются теоремы о существовании и единственности решения в зависимости от тех или иных условий.
Жесткая система ОДУ химической кинетики
Листинг 9.10. Жесткая система ОДУ химической кинетики
Бросается в глаза сильно различающийся порядок коэффициентов при разных слагаемых. Именно степень этого различия чаще всего и определяет жесткость системы ОДУ. В качестве соответствующей характеристики выбирают матрицу Якоби (якобиан) векторной функции F(t,y), т.е. функциональную матрицу, составленную из производных F(t,y) (см. разд 3.4.3). Чем сильнее вырождена матрица Якоби, тем жестче система уравнений. В приведенном примере определитель якобиана и вовсе равен нулю при любых значениях у0, y1 и у2 (листинг 9.11, вторая строка). В первой строке листинга 9.11 приведено напоминание способа вычисления якобиана средствами Mathcad на примере определения элементов его первой строки.
Якобиан рассматриваемой системы ОДУ химической кинетики
Листинг 9.11. Якобиан рассматриваемой системы ОДУ химической кинетики
Для примера, приведенного в листинге 9.10, стандартным методом Рунге— Кутты все-таки удается найти решение (оно показано на Рисунок 9.15). Однако для этого требуется очень большое число шагов, M=20000, что делает расчеты очень медленными. При меньшем числе шагов численному алгоритму не удается найти решение. В процессе работы алгоритма оно расходится, и Mathcad вместо результата выдает ошибку о превышении предельно большого числа.
Еще один факт, на который стоит обратить внимание, — это различие в порядке величины получающегося решения. Как видно из Рисунок 9.15, концентрация первого реагента уг существенно (в тысячи раз) превышает концентрацию остальных. Это свойство также очень характерно для жестких систем.
Примечание 1
Примечание 1
В принципе, можно было бы снизить жесткость системы "вручную", применяя масштабирование. Для этого нужно искусственно уменьшить искомую функцию у1 к примеру, в тысячу раз, разделив все слагаемые в системе ОДУ, содержащие y1, на 1000. После масштабирования для решения полученной системы методом Рунге—Кутты будет достаточно взять всего м=20 шагов. Соответствующий пример вы найдете на компакт-диске.
Решение жесткой системы ОДУ химической кинетики
Листинг 9.12. Решение жесткой системы ОДУ химической кинетики
Расчеты показывают, что для получения того же результата (который был изображен на Рисунок 9.15) оказалось достаточно в тысячу раз меньшего количества шагов численного алгоритма, чем для стандартного метода Рунге— Купы ! Примерно во столько же раз требуется меньше компьютерного времени на проведение расчетов. Стоит ли говорить, что, если вы имеете дело с жесткими (в той или иной степени) системами, применение описанных специальных алгоритмов просто необходимо.
Важно заметить, что до сих пор мы имели дело с примером не очень жесткой системы. Попробуйте вместо скоростей упомянутых химических реакций 0.1, 103 и 102 взять другие числа, например, 0.05, 104 и 107 соответственно. Заметим, что такое соотношение скоростей часто встречается в прикладных задачах химической кинетики и определяет куда более жесткую систему ОДУ. Ее уже никак не удается решить стандартными методами, поскольку число шагов численного метода должно быть просто гигантским. А между тем алгоритмы для жестких ОДУ справляются с этой задачей с легкостью (Рисунок 9.16), причем практически при тех же значениях шага, что были взяты в листинге 9.12. Обратите внимание, что порядки величины решений для концентраций различных веществ на Рисунок 9.16 различаются еще сильнее, чем в предыдущем (менее жестком) примере.
Модель "хищникжертва"
Листинг 9.13. Модель "хищник-жертва"
Модель замечательна тем, что в такой системе наблюдаются циклическое увеличение и уменьшение численности и хищника (Рисунок 9.17), и жертвы, так часто наблюдаемое в природе. Фазовый портрет системы представляет собой концентрические замкнутые кривые, окружающие одну стационарную точку, называемую центром. Как видно, модельные колебания численности обеих популяций существенно зависят от начальных условий — после каждого периода колебаний система возвращается в ту же точку. Динамические системы с таким поведением называют негрубыми.
Примечание 1
Примечание 1
Пример негрубой системы (модель осциллятора без затухания) с особой точкой типа "центр" был нами рассмотрен ранее (см. разд. 9.1).
Модель Ван дер Поля
Листинг 9.14. Модель Ван дер Поля
Модель Лоренца
Листинг 9.15. Модель Лоренца
Решением системы Лоренца при определенном сочетании параметров (Рисунок 9.21 и 9.22) является странный аттрактор (или аттрактор Лоренца) — притягивающее множество траекторий на фазовом пространстве, которое по виду идентично случайному процессу. В некотором смысле аттрактор Лоренца является стохастическими автоколебаниями, которые поддерживаются в динамической системе за счет внешнего источника.
Решение в виде странного аттрактора появляется только при некоторых сочетаниях параметров. В качестве примера на Рисунок 9.23 приведен результат для r=10 и тех же значений остальных параметров. Как видно, аттрактором в этом случае является фокус. Перестройка типа фазового портрета происходит в области промежуточных г. Критическое сочетание параметров, при которых фазовый портрет системы качественно меняется, называется в теории динамических систем точкой бифуркации. Физический смысл бифуркации в модели Лоренца, согласно современным представлениям, описывает переход ламинарного движения жидкости к турбулентному.
Построение фазового портрета для модели брюсселятора
Листинг 9.16. Построение фазового портрета для модели брюсселятора
Предложенный алгоритм формирует из отдельных матриц решений системы ОДУ с разными начальными условиями объединенную матрицу U. Пары начальных условий задаются в первой строке листинга в виде матрицы v размера 2х10. Последнее означает построение десяти траекторий. Для того чтобы поменять количество траекторий, измените соответствующим образом размер этой матрицы. Затем (Рисунок 9.24) элементы матрицы и выводятся на график в виде отдельных точек. Отсутствие соединения точек линиями является недостатком алгоритма, но это плата за возможность представить в Mathcad несложным образом сразу большое количество траекторий на фазовой плоскости.
решение задачи Коши для ОДУ второго порядка (модель нелинейного осциллятора)
Листинг 9.2. решение задачи Коши для ОДУ второго порядка (модель нелинейного осциллятора)
Пользователь имеет возможность выбирать между двумя модификациями численного метода Рунге—Кутты. Для смены метода необходимо нажатием правой кнопки мыши на области функции odesoive вызвать контекстное меню и выбрать в нем один из трех пунктов: Fixed (С фиксированным шагом), Adaptive (Адаптивный) или Stiff (Для жестких ОДУ). Подробнее о различиях этих методов сказано в разд. 9.3.4 и 9.4.
Решение системы двух ОДУ(модель осциллятора)
Листинг 9.3. Решение системы двух ОДУ(модель осциллятора)
Решение задачи Коши для ОДУ первого порядка
Листинг 9.4. Решение задачи Коши для ОДУ первого порядка.
Поиск аттрактора системы двух ОДУ модели осциллятора
Листинг 9.5. Поиск аттрактора системы двух ОДУ модели осциллятора
ВНИМАНИЕ!
Функции bulstoer и rkadapt (те, что пишутся со строчной буквы) не предназначены для нахождения решения в промежуточных точках интервала, хотя они и выдают их в матрице-результате.
Решение системы ОДУ как функции числа шагов
Листинг 9.6. Решение системы ОДУ как функции числа шагов
Погрешность алгоритмов решения ОДУ в точке
Обратимся теперь к функциям buistoer и rkadapt (тем, что пишутся со строчной буквы), которые предназначены для нахождения решения в определенной точке интервала интегрирования системы ОДУ. В первую очередь, подчеркнем еще раз, что их нельзя применять для получения решения в промежуточных точках интервала, несмотря на их вывод в матрице-результате. На Рисунок 9.10 показан фазовый портрет системы ОДУ динамической модели осциллятора, полученный при помощи функции buistoer в листинге 9.5 (см. разд. 9.3.3) и с помощью rkadapt (при соответствующей замене третьей строки листинга 9.5). Видно, что, несмотря на высокую точность, равную 10-5, и верный результат на конце интервала, график, полученный при помощи функции buistoer (Рисунок 9.10, а), мало напоминает правильный фазовый портрет (см. Рисунок 9.9), начиная быть приемлемым только при предельно допустимом для обсуждаемых функций значении асс=10-16. В данном случае конкретной системы ОДУ функция rkadapt дает качественно верное решение (Рисунок 9.10, б), однако его необходимая точность обеспечивается только в последней точке временного интервала.
В заключение остановимся на влиянии выбора параметра асе на расчеты. Как правило, разные численные методы работают несколько по-разному. Алгоритм Рунге—Кутты дает результат тем ближе к истинному, чем меньше выбирается параметр асе, а Булирша—Штера часто демонстрирует менее естественную зависимость у(е): даже при относительно больших е реальная погрешность остается хорошей (намного лучше метода Рунге—Кутты).
Число шагов адаптивного численного алгоритма в зависимости от погрешности acc
Листинг 9.7. Число шагов адаптивного численного алгоритма в зависимости от погрешности acc
Таким образом, проводя тестовые расчеты для различных задач и подбирая наилучший набор параметров, можно существенно сэкономить ресурсы компьютера. Конечно, проводить подобный анализ стоит в случаях, когда время расчетов для вас критично.
Решение нежесткого ОДУ
Листинг 9.8. Решение нежесткого ОДУ
Решение жесткого ОДУ алгоритмом RADAUS
Листинг 9.9. Решение жесткого ОДУ алгоритмом RADAUS
Неверное решение более жесткого ОДУ методом Рунге—Кутты (листин 9 8 с коэффициентом 50)
Рисунок 9.13. Неверное решение более жесткого ОДУ методом Рунге—Кутты (листин 9.8 с коэффициентом -50)
Таким образом, во-первых, мы выяснили, что одни и те же уравнения с разными параметрами могут быть как жесткими, так и нежесткими. Во-вторых, чем жестче уравнение, тем больше шагов в обычных численных методах требуется для его устойчивого решения. С классическим примером ОДУ из листинга 9.8 все получилось хорошо, т. к. оно было не очень жестким, и небольшое увеличение числа шагов разрешило все проблемы. Для решения обычными методами более жестких уравнений требуются миллионы, миллиарды и даже большее число шагов.
Примечание 1
Примечание 1
Некоторые ученые замечают, что в последние годы методы Рунге—Кутты, считавшиеся ранее незыблемыми, стали уступать свое главенствующее положение среди алгоритмов решения ОДУ методам, способным решать жесткие задачи.
О постановке задач
9.1. О постановке задачНачнем разговор об обыкновенных дифференциальных уравнениях с формулировки типичных задач, сопровождая вопросы их постановки конкретными примерами. При этом мы будем, несколько забегая вперед, приводить их решение, не обсуждая особенностей его отыскания средствами Mathcad.
Обыкновенные дифференциальные уравнения динамические системы
Обыкновенные дифференциальные уравнения: динамические системыВ этой главе рассматриваются численные методы решений задач с начальными условиями (называемых задачами Коши) для обыкновенных дифференциальных уравнений (далее используется сокращение ОДУ). Такие задачи требуют нахождения функции (или нескольких функций) одной переменной, если, во-первых, определено дифференциальное уравнение (или система уравнений), содержащее производную функции, и, во-вторых, необходимое количество дополнительных условий, задающих значение функции в некоторой начальной точке.
Решение задач Коши для ОДУ — давно и детально разработанная технология. С "хорошими" ОДУ вообще никаких вычислительных проблем обычно не возникает (чаще всего они решаются при помощи алгоритма Рунге—Купы), а для ОДУ особого типа, называемых жесткими, необходимо применять специальные методы. Все эти возможности заложены в Mathcad, причем пользователю позволено выбирать конкретный алгоритм решения ОДУ.
ы классические динамические системы
9.5. Примеры: классические динамические системыВ предыдущих разделах было использовано в качестве примера в основном линейное уравнение осциллятора (оно содержало только первую степень неизвестных функций и их производных). Между тем многие нелинейные уравнения демонстрируют совершенно удивительные свойства, причем решение большинства из них можно получить лишь численно.
Рассмотрим несколько наиболее известных классических примеров динамических систем, имея в виду, что читателю они могут пригодиться как в познавательных, так и в практических целях. Это модели динамики популяций (Вольтерры), генератора автоколебаний (Ван дер Поля), турбулентной конвекции (Лоренца) и химической реакции с диффузией (Пригожина). Как уже было сказано, для изучения динамических систем разработана специальная теория, центральным моментом которой является анализ фазовых портретов, т. е. решений, получающихся при выборе всевозможных начальных условий.
Примечание 1
Примечание 1
В большинстве примеров, изложенных ниже, для построения схемы фазового портрета рассчитывается несколько решений для разных начальных условий. О, том, как проделать такие расчеты в Mathcad, будет рассказано ниже на примере модели брюсселятора (см. разд. 9.5.4).
Ограничимся в дальнейшем минимальными комментариями и приведем листинги и графики решений без подробного обсуждения.
Решение более жесткой системы ОДУ химической кинетики методом Розенброка
Рисунок 9.16. Решение более жесткой системы ОДУ химической кинетики методом Розенброка
Для решения очень жестких систем особенно подходит функция Radau, которая соблазнительна еще и тем, что избавляет от необходимости предварительного расчета якобиана. В случае очень жесткой задачи химической кинетики результат, выдаваемый функцией Radau, практически тот же, что и в случае использования алгоритма Розенброка (Рисунок 9.16). Любопытно, что решение менее жесткой системы (из предыдущего листинга) при помощи функции Radau удается получить только для первой половины интервала, примерно до t1=20. Для больших интервалов она дает сбой, выводя вместо решения сообщение об ошибке.
Решение нежесткого ОДУ методом Рунге—Кутты (продолжение листинга 9 8)
Рисунок 9.12. Решение нежесткого ОДУ методом Рунге—Кутты (продолжение листинга 9.8)
Обратите внимание на значение коэффициента -10 во второй строке листинга. Если изменить его на -50, то решение меняется до неузнаваемости (Рисунок 9.13). Вас, несомненно, должен насторожить результат, выданный Mathcad. Характерная "разболтка" решения говорит о неустойчивости алгоритма. Первое, что можно сделать, — увеличить количество шагов в методе Рунге—Куттты. Для этого достаточно добавить третий параметр step в функцию odesoive (t, i,step). После нескольких экспериментов можно подобрать такое значение step, которое будет обеспечивать устойчивость решения. Читатель может самостоятельно убедиться, что при step>50 "разболтка" пропадает, и решение принимает вид графика, показанного ниже на Рисунок 9.14.
что решение подобных нелинейных динамических
Рисунок 9.23. Решение системы Лоренца с измененным параметром г=10
Замечательно, что решение подобных нелинейных динамических систем можно получить только численно, поэтому их изучение стало бурно развиваться с ростом возможностей вычислительной техники в последние полвека.
Примечание 1
Примечание 1
На компакт-диске вы найдете решение в виде странного аттрактора еще одной классической динамической системы (палеомагнетизм Рикитаке), моделирующей глобальное движение магнитных полюсов Земли.
Примечание 2
Примечание 2
Также на компакт-диске находится программа с реализацией алгоритма поиска отображения Пуанкаре (на примере модели Лоренца) — эффективного приема теории динамических систем.
Решение системы ОДУ в зависимости от числа узлов М (продолжение листинга 9 6)
Рисунок 9.8. Решение системы ОДУ в зависимости от числа узлов М (продолжение листинга 9.6)
Решение уравнения со2у' '+у=0 для различных начальных условий (коллаж графиков)
Рисунок 9.3. Решение уравнения со2-у' '+у=0 для различных начальных условий (коллаж графиков)
Резюмируя содержание вводного раздела главы, перечислим еще раз типичные постановки задач, характерные для динамических систем:
В дальнейших разделах этой главы при рассказе о возможностях Mathcad мы будем в первую очередь описывать решение первой (базовой) задачи, для которой предусмотрен целый арсенал средств. А именно: вычислительный блок для решения ОДУ (см. разд. 9.2), несколько встроенных функций для решения систем ОДУ (см. разд. 9.3), в том числе жестких, которые не поддаются решению стандартными методами (см. разд. 9.4). Приемы, которые автор рекомендует в качестве идиом решения остальных задач, будут рассмотрены эпизодически, на конкретных примерах классических динамических систем вычислительной физики, химии и биологии (см. разд. 9.5 и 9.4.3), связанных с динамическими системами. В частности, программа для визуализации фазового портрета рассмотрена в конце главы, на примере модели брюсселятора (см. разд. 9.5.4). Сводка алгоритмов с рекомендациями по их применению в зависимости от типа задачи приведена конспективно, без детального разбора (см. разд. 9.3.4).
Решение уравнения у'=1у2 (продолжение листинга 9 4)
Рисунок 9.7. Решение уравнения у'=1-у2 (продолжение листинга 9.4)
Построение графика (Рисунок 9.7) осуществляется так же, как и в рассмотренном в предыдущем разделе случае N уравнений, при помощи выделения столбцов из матрицы решения посредством оператора <>
Решение уравнения Ван дер Поля
Рисунок 9.20. Решение уравнения Ван дер Поля при других начальных условиях у=-2, у' =-3
Решение уравнения w2у' '+?у'+у=0 на фазовой плоскости (продолжение листинга 9 1)
Рисунок 9.2. Решение уравнения w2у' '+?у'+у=0 на фазовой плоскости (продолжение листинга 9.1)
Как правило, решение задач Коши для ОДУ и их систем — задача хорошо разработанная и с вычислительной точки зрения довольно простая. На практике чаще встречаются другие, более сложные задачи, в частности, исследование поведения динамической системы в зависимости от начальных условий. При этом в большинстве случаев бывает необходимым изучить только асимптотическое решение ОДУ, т.е. y(t->?), называемое аттрактором. Очень наглядным образом можно визуализировать такую информацию на фазовой плоскости, во многом благодаря тому, что существует всего несколько типов аттракторов, и для них можно построить четкую классификацию.
С одной стороны, каждое решение будет выходить из точки, координаты которой являются начальными условиями, но, оказывается, для большинства ОДУ целые семейства траекторий будут заканчиваться в одних и тех же аттракторах (стационарных точках или предельных циклах). Множество решений, вычисленное для всевозможных начальных условий, образует фазовый портрет динамической системы. С вычислительной точки зрения задача исследования фазового портрета часто сводится к обычному сканированию семейств решений ОДУ при разных начальных условиях.
Примечание 1
Примечание 1
Для рассматриваемого примера модели гармонического осциллятора имеется единственная стационарная точка (аттрактор), на которую "накручивается" решение, из каких бы начальных условий оно ни выходило. В теории динамических систем аттрактор такого типа называется фокусом.
Дальнейшее усложнение задач анализа фазовых портретов связано с их зависимостью от параметров, входящих в систему ОДУ. В частности, при плавном изменении параметра модели может меняться расположение аттракторов на фазовой плоскости, а также могут возникать новые аттракторы и прекращать свое существование старые. В первом случае, при отсутствии особенностей, будет происходить простое перемещение аттракторов по фазовой плоскости (без изменения их типов и количества), а во втором — фазовый портрет динамической системы будет коренным образом перестраиваться. Критическое сочетание параметров, при которых фазовый портрет системы качественно меняется, называется в теории динамических систем точкой бифуркации.
Поясним понятие бифуркации на примере той же модели осциллятора, которая зависит от двух параметров (ш и р). При р>о существует единственная стационарная точка типа фокуса (см. Рисунок 9.2), которая в точке бифуркации Р=о вырождается в аттрактор типа центр, характеризующийся тем, что решения ОДУ представляют собой циклы, совершаемые вокруг этой точки с амплитудой, которая существенно зависит от начальных условий (Рисунок 9.3). Для надежного исследования фазового портрета практически всегда необходимо решить систему ОДУ большое количество раз с самыми разными начальными условиями (и, возможно, с разным набором параметров модели), чтобы посмотреть, к каким аттракторам сходятся различные траектории.
Решение уравнения w2у"+?у'+у=0 (продолжение листинга 9 1)
Рисунок 9.1. Решение уравнения w2у"+?у'+у=0 (продолжение листинга 9.1)
Имеются два типа задач, которые возможно решать с помощью Mathcad:
Сказанное относится как к отдельным дифференциальным уравнениям (см. разд. 9.2), так и к системам ОДУ (см. разд. 9.3). Важно подчеркнуть, что Mathcad умеет решать только такие системы дифференциальных уравнений, которые могут быть представлены в стандартной форме (форме Коши). Для неизвестных функций у0 (t), y1(t), ..., yN-1(t) система ОДУ должна быть записана в форме:
что эквивалентно следующему векторному представлению:
Y' (t)=F(Y(t) ,t). (9.2)
Здесь Y и Y' — соответствующие неизвестные векторные функции переменной t размера Nx1, a F — векторная функция того же размера и количества переменных (N+1) (N компонент вектора и, возможно, t). Именно векторное представление (9.2) используется для ввода системы ОДУ в среде Mathcad.
Для того чтобы определить задачу Коши для системы из N ОДУ первого порядка, следует определить еще ровно N начальных условий, задающих значение каждой из функций yi(t0) в начальной точке интегрирования системы t0. В векторной форме они могут быть записаны в виде
Y(t0)=B, (9.3)
где B— вектор начальных условий размера Nx1, составленный из yi(t0). Обратите внимание на необходимость векторной записи как самого уравнения, так и начального условия. В случае одного ОДУ первого порядка соответствующие векторы имеют только один элемент, а в случае системы N>1 уравнений — N.
Стандартные процедуры Mathcad применимы для систем ОДУ первого порядка, записанных в форме (9.2)—(9.3). Тем не менее, если в систему входят и уравнения высших порядков, то ее можно свести к системе большего числа уравнений первого порядка. Рассмотрим в качестве примера уравнение второго порядка модели осциллятора w2у'+?у'+у=0, решенное в листинге 9.1. Если ввести обозначение y0(t)=y(t), yi(t)=y' (t), то уравнение сведется к эквивалентной системе:
y0(t)=y1(t);
w2у1'+?у'+у=0,
форма которой удовлетворяет форме (9.2) и уравнение может быть решено средствами Mathcad, предназначенными для систем ОДУ. Именно эта система решена ниже в листинге 9.3 (см. разд. 9.3.1). Отметим, что на Рисунок 9.1 показаны как зависимость у (t), так и у1 (t)=y' (t).
Решение уравнения w2у''+?у '+у+?у2=0 (продолжение листинга 9 2)
Рисунок 9.4. Решение уравнения w2у''+?у '+у+?у2=0 (продолжение листинга 9.2)
Еще раз подчеркнем, что результатом применения блока Given/odesolve является функция y(t), определенная на промежутке (to,ti). Следует воспользоваться обычными средствами Mathcad, чтобы построить ее график или получить значение функции в какой-либо точке указанного интервала, например: у(10)=0.048.
Решение в виде аттрактора Лоренца (продолжение листинга 9 15)
Рисунок 9.21. Решение в виде аттрактора Лоренца (продолжение листинга 9.15)
Решение жесткого ОДУ методом RADAUS (продолжение листингов 9 8 и 9 9)
Рисунок 9.14. Решение жесткого ОДУ методом RADAUS (продолжение листингов 9.8 и 9.9)
В заключение приведем соответствующие встроенные функции, которые применяются для решения жестких систем ОДУ не на всем интервале, а только в одной заданной точке t1.
Имена этих функций пишутся со строчной буквы, а их действие и набор параметров полностью аналогичны рассмотренным нами ранее для функций, относящихся к решению в заданной точке нежестких систем (см. разд. 9.3.2). Отличие заключается в специфике применяемого алгоритма и необходимости задания матричной функции якобиана J(t,y) (для двух последних функций).
Решение жесткой системы ОДУ химической кинетики методом Рунге—Кутты (продолжение листинга 9 10)
Рисунок 9.15. Решение жесткой системы ОДУ химической кинетики методом Рунге—Кутты (продолжение листинга 9.10)
Применение специфических алгоритмов
Покажем действие этих алгоритмов на том же примере жесткой системы ОДУ химической кинетики (листинг 9.12). Обратите внимание, как следует представлять в данном случае якобиан, сравнив задание матричной функции в предпоследней строке листинга 9.12 с заданием якобиана из листинга 9.11.
Результат выдаваемый встроенной функцией в качестве решения системы ОДУ (продолжение листинга 9 3}
Рисунок 9.5. Результат, выдаваемый встроенной функцией в качестве решения системы ОДУ (продолжение листинга 9.3}
Первая строка листинга представляет задание параметров модели, вторая — начального условия задачи Коши, а в третьей строке листинга определено число шагов, на которых рассчитывается решение. Самая важная — это предпоследняя строка листинга, в которой, собственно, определяется система ОДУ. Последняя строка присваивает матричной переменной и результат действия функции rkfixed. Решение системы ОДУ будет осуществлено на промежутке (0,40).
Сравните рассматриваемую систему (см. разд. 9.1.Т), записанную в стандартной форме с формальной ее записью в Mathcad, чтобы не делать впоследствии ошибок. Во-первых, функция D, входящая в число параметров встроенных функций для решения ОДУ, должна быть функцией обязательно двух аргументов. Во-вторых, второй ее аргумент должен быть вектором того же размера, что и сама функция о. В-третьих, точно такой же размер должен быть и у вектора начальных значений у0. Не забывайте, что векторную функцию D(t,y) следует определять через компоненты вектора у с помощью кнопки Subscript (Нижний индекс) с наборной панели Calculator (Калькулятор) или нажатием клавиши <[>.
Матрица, представляющая решение, показана на Рисунок 9.5. Размер полученной матрицы будет равен (M+i)x(N+l), т. е. 51x3. Просмотреть все компоненты матрицы и, которые не помещаются на экране, можно с помощью вертикальной полосы прокрутки. Как нетрудно сообразить, на этом рисунке отмечены выделением расчетные значения искомых векторов на 10-м шаге. Это соответствует, с математической точки зрения, найденным значениям уо(8.0)=о. 307 и yi (8.0) =0.204. Для вычисления элементов решения в последней точке интервала используйте вывод элемента ин, ь Для построения графика решения надо отложить соответствующие компоненты матрицы решения по координатным осям: значения аргумента и<0> — вдоль оси х, а u<1> и u<2> — вдоль оси Y (Рисунок 9.6).
Система N дифференциальных уравнений
9.3. Система N дифференциальных уравненийПри помощи Mathcad можно решать системы N>1 ОДУ первого порядка, если они записаны в стандартной форме (Коши) в виде векторного соотношения: Y' (t)=F(Y(t),t) (см. разд. 9.1. 1).
Зависимость числа шагов от параметра асе численных методов (продолжение листинга 9 7)
Рисунок 9.11. Зависимость числа шагов от параметра асе численных методов (продолжение листинга 9.7)
Жесткие системы ОДУ
9.4. Жесткие системы ОДУДо сих пор мы имели дело с "хорошими" уравнениями, которые надежно решались численными методами Рунге—Кутты. Однако имеется класс так называемых жестких (stiff) систем ОДУ, для которых стандартные методы практически неприменимы, поскольку их решение требует исключительно малого значения шага численного метода. Некоторые из специальных алгоритмов, разработанных для этих систем, реализованы в Mathcad.
Матричные вычисления в Mathcad
Разностном методе
10.4.1.0 разностном методеРазберем идею разностного метода решения краевых задач на примере взаимодействия световых пучков (см. Рисунок 10.1), переобозначив в системе (10.1) интенсивность излучения вправо на Y, а интенсивность излучения влево на у (просто в целях удобства, чтобы не писать в дальнейшем индекс). Суть метода заключается в покрытии расчетного интервала сеткой из N точек. Тем самым определяются (м-l) шагов (Рисунок 10.7). Затем надо заменить дифференциальные уравнения исходной краевой задачи аппроксимирующими их уравнениями в конечных разностях, выписав соответствующие разностные уравнения для каждого 1-го шага. В нашем случае достаточно просто заменить первые производные из (10.1) их разностными аналогами (такой метод называется еще методом Эйлера):
Алгоритм стрельбы
10.2.1. Алгоритм стрельбыСуть метода стрельбы заключается в пробном задании недостающих граничных условий на левой границе интервала и решении затем полученной задачи Коши хорошо известными методами (см. главу 11). В нашем примере не хватает начального условия для y1(0), поэтому сначала зададим ему произвольное значение, например, у1(0)=10. Конечно, такой выбор не совсем случаен, поскольку из физических соображений ясно, что, во-первых, интенсивность излучения — величина заведомо положительная, и, во-вторых, отраженное излучение должно быть намного меньше падающего. Решение задачи Коши с помощью функции rkfixed приведено в листинге 10.1, причем в его последней строке выведены расчетные значения у0 и уг на правой крайней точке интервала х=1. Разумеется, они не совпадают, т. е. правое краевое условие не выполнено, из чего следует, что полученный результат не является решением поставленной краевой задачи.
О постановке задачи
10.5.1. О постановке задачиНа примере решения нелинейной задачи взаимодействия световых пучков при помощи алгоритма стрельбы разберемся с физикой явления, которую описывает нелинейность. Рассмотрим простейшее усложнение модели (10.1) (см. разд. 10.1), добавив в уравнения малые слагаемые, зависящие квадратично от интенсивности света. Чтобы быть ближе к физической сущности явления, мы введем эту зависимость в формулу коэффициента поглощения среды а(х), сделав а(х) из (10.1) функцией не только координаты, но и суммарной интенсивности:
a(x,Y,y)=a0(x)-?(x)(Y+y), (10.8)
где а0(х) — постоянная (не зависящая от интенсивности света) составляющая коэффициента поглощения, а е(х) — малый коэффициент пропорциональности. В этом случае модель (10.1) перепишется в виде:
Таким образом, физическая подоплека квадратичной нелинейности задачи (10.9) связана с зависимостью коэффициента поглощения среды от Y+y. Иными словами, чем мощнее свет в некоторой точке среды, тем меньше локальный коэффициент поглощения среды в этой точке.
Механизм этого явления будет хорошо понятен, если предположить, что средой является воздух, насыщенный парами жидкости (попросту говоря, туман). Чем сильнее в некоторой точке туман, тем больше локальное поглощение. Если свет, проходящий через туман, имеет малую интенсивность, то он не оказывает на среду никакого (или почти никакого) воздействия. Однако если пучок света будет очень мощным, за счет его поглощения пары воды будут разогреваться и, следовательно, испаряться. По мере испарения капель воды туман станет рассеиваться, и среда просветлеет (т. е. ее коэффициент поглощения уменьшится). Таким образом, чем мощнее свет, тем активнее пойдет процесс испарения, тем существеннее уменьшится коэффициент поглощения. Именно такой эффект, совершенно понятный с физической точки зрения, и описывает формула (10.8) и, соответственно, модель (10.9).
Подчеркнем еще раз и физический смысл коэффициентов а0(х) и ?(х). Первый описывает постоянный (не зависящий от интенсивности света) коэффициент поглощения среды, а второй — его зависимость от эффекта разогрева среды светом. Чем больше коэффициент а0(х), тем более жесткой является краевая задача; а чем больше ?(х), тем сильнее ее нелинейность. В предельном случае ?(х)->0 задача (10.9) переходит в задачу (10.1), становясь линейной.
Двухточечные краевые задачи
10.2.2. Двухточечные краевые задачиРешение краевых задач для систем ОДУ методом стрельбы в Mathcad достигается применением двух встроенных функций. Первая предназначена для двухточечных задач с краевыми условиями, заданными на концах интервала:
ВНИМАНИЕ!
Решение краевых задач в Mathcad, в отличие от большинства остальных операций, реализовано не совсем очевидным образом. В частности, помните, что число элементов векторов о и load равно количеству уравнений N, а векторов z, score и результата действия функции sbval— количеству правых граничных условий L. Соответственно, левых граничных условий в задаче должно быть (N-L) .
Как видно, функция sbval предназначена не для поиска собственно решения, т. е. неизвестных функций yi(x), а для определения недостающих начальных условий в первой точке интервала, т. е. yi(x0). Чтобы вычислить yi(х) на всем интервале, требуется дополнительно решить задачу Коши. Разберем особенности использования функции sbval на конкретном примере (листинг 10.2), описанном выше (см. разд. 10.1.1). Краевая задача состоит из системы двух уравнений (N=2), одного левого (L=1) и одного правого (N-L=2-1=1) граничного условия.
Метод стрельбы
10.5.2. Метод стрельбыВстроенная функция sbval, реализующая в Mathcad алгоритм стрельбы (см. разд. 10.2), может справляться и с нелинейными задачами. Приведем пример решения краевой задачи (10.9) (листинг 10.8 и Рисунок 10.11) с теми же граничными условиями, что были поставлены для модели (10.1). Интерес представляет решение у, описывающее рост интенсивности отраженного пучка по мере его распространения справа-налево (обратите внимание, что на Рисунок 10.11 оно отложено с десятикратным увеличением).
Жесткие краевые задачи
10.4.2. Жесткие краевые задачиОдин из случаев, когда применение разностных схем может быть очень полезным, связан с решением жестких краевых задач (подробнее о жестких ОДУ читайте в главе 11). В частности, рассматриваемая задача о встречных световых пучках становится жесткой при увеличении коэффициента ослабления а(х) в несколько десятков раз. Например, при попытке решить ее с а(х) :=100 с помощью листинга 10.2 вместо ответа выдается сообщение об ошибке "Can't converge to a solution. Encountered too many integration steps" (He сходится к решению. Слишком много шагов интегрирования). Это и неудивительно, поскольку жесткие системы характерны тем, что требуют исключительно малого значения шага в стандартных алгоритмах.
Для жестких задач неприменимы и явные разностные схемы, о которых рассказывалось в предыдущем разделе (см. разд. 10.3.1). Результат расчетов по программе листинга 10.6, например с а(х) :=20 (Рисунок 10.9), дает характерную для неустойчивых разностных схем "разболтку" — колебания нарастающей амплитуды, не имеющие ничего общего с реальным решением.
Краевые задачи с условием во внутренней точке
10.2.3. Краевые задачи с условием во внутренней точкеИногда дифференциальные уравнения определяются с граничными условиями не только на концах интервала, но и с дополнительным условием в некоторой промежуточной точке расчетного интервала. Чаще всего такие задачи содержат данные о негладких в некоторой внутренней точке интервала решениях. Для них имеется встроенная функция bvaifit, также реализующая алгоритм стрельбы.
Обычно функция bvaifit применяется для задач, в которых производная решения имеет разрыв во внутренней точке xf. Некоторые из таких задач не удается решить обычным методом пристрелки, поэтому приходится вести пристрелку сразу из обеих граничных точек. В этом случае дополнительное внутреннее условие в точке xf является просто условием сшивки в ней левого и правого решений. Поэтому для данной постановки задачи оно записывается В форме score (xf, у) :=у.
Рассмотрим действие функции bvaifit на знакомом примере модели взаимодействия пучков света (см. Рисунок 10.1), предположив, что в промежутке между xf=0.5 и x1=1.0 находится другая, оптически более плотная среда с другим коэффициентом ослабления излучения а(х)=3. Соответствующая краевая задача решена в листинге 10.3, причем разрывный показатель ослабления определяется в его второй строке.
Разностные схемы
10.5.3. Разностные схемыРешение краевых задач при помощи разностных схем связано с необходимостью разработки собственного алгоритма для каждой конкретной задачи. Как вы помните (см. разд. 10.4), в случае линейных уравнений в результате построения разностной схемы система алгебраических сеточных уравнений также получалась линейной. Это автоматически означало, что она имеет единственное решение, которое в большинстве случаев могло быть найдено стандартными численными методами. В ситуации, когда исходные дифференциальные уравнения нелинейны, система разностных уравнений тоже является нелинейной, и их решение существенно усложняется, хотя бы потому, что оно не является единственным. Поэтому подход к построению разностных схем нелинейных уравнений должен быть специфическим, но наградой за него станет решение задач, с которыми не справляется алгоритм стрельбы (например, жестких).
Подробная разработка алгоритмов и соответствующих программ Mathcad выходит далеко за пределы данной книги, поэтому мы конспективно представим один из приемов решения нелинейных краевых задач, сводящийся к их линеаризации. В общих чертах подход заключается в следующем. Предположим, что некоторое приближение (обозначим его Y(х) и у(х)) к решению нелинейной задачи (10.9) нам известно, и можно считать, что Y->Y+Z и y->y+z, где z и z — близкие к нулю функции х. Тогда, пользуясь их малостью (по сравнению с Y0 и у0), можно разложить нелинейные слагаемые в уравнениях (10.9) в ряд Тейлора по z и z. Получим:
Y'+Z' = -aY + ry + ?(Y2 + Yy) - aZ + rz + 2?Y(Z + z); (10.10)
y'+z' = ay - rY - ?(y2 + Yy) + az - rz - 2?y(Z + z) .
Теперь, поскольку Y(x) и у(х) являются приближением к решению исходной задачи, то можно считать, что они (приблизительно) удовлетворяют и уравнению, и граничным условиям. Тогда, вычитая (10.9) из (10.10), получим краевую задачу для z (х) и z (х):
Z' = -aZ + rz + 2?Y(Z + z);
z' = az - rZ - 2?y(Z + z) ; (10.11)
Z(0) = z(0) = Z(l) = z(l) = 0.
Это и есть та самая новая задача для "маленьких" функций z (х) и z (х), которую надо решить, и которая, благодаря малости z и z, является линейной. Вся беда в том, что мы не знаем "больших" функций Y(x) и у(х), а они, увы, входят в задачу (10.11). Тем не менее рецепт получения этих функций напрашивается сам собой: если нелинейность исходных ОДУ не слишком сильная, то в качестве Y и у можно взять решение линейной краевой задачи, т. е. задачи (10.9) с ?(х)=0 (см. разд. 10.4.1).
Сказанное иллюстрируют листинги 10.9 и 10.10, первый из которых решает линейную краевую задачу (являясь, фактически, небольшой модификацией листинга 10.7 из разд. 10.4.2), а второй решает линеаризованную задачу (10.11), учитывая результат листинга 10.9. В листинге 10.9 матрица о и вектор правых частей в являются разностной аппроксимацией ОДУ (ее первые N строк аппроксимируют первое уравнение, а оставшиеся N строк — второе). Такой же смысл и точно такую же структуру имеют матрица с и вектор правых частей с для второй (линеаризованной) задачи (10.11). Для решения сеточных уравнений Dу=B и Cz=G используется (конечно, весьма неэкономично) встроенная функция isolve, реализующая алгоритм Гаусса.
Важно привлечь внимание читателя к последним строкам листинга 10.9. В них осуществляется интерполяция полученного решения системы сеточных уравнений для того, чтобы в нелинейной задаче (в листинге 10.10) можно было использовать непрерывные "большие" функции из линейной задачи. В последней строке листинга 10.10 осуществляется сложение "больших" и "маленьких" функций (результатов листингов 10.9 и 10.10) для получения полного решения нелинейной задачи (10.9), которое изображено на Рисунок 10.12.
Не будем давать дополнительных комментариев к Mathcad-программам, надеясь, что читатель, заинтересовавшийся нелинейным примером со световыми пучками и эффектом разогрева светом среды, сам разберется в листингах, тем более что техника разработки разностных схем была нами детально разобрана раньше (см. разд. 10.4).
Обыкновенные дифференциальные уравнения краевые задачи
Глава 10. Обыкновенные дифференциальные уравнения: краевые задачи| Содержание |
Иллюстрация метода стрельбы (продолжение листинга 10 1)
Рисунок 10.2. Иллюстрация метода стрельбы (продолжение листинга 10.1)
Решение пробной задачи Коши для модели (10 1)
Листинг 10.1. Решение пробной задачи Коши для модели (10.1)
График полученных решений показан на Рисунок 10.2 (слева). Из него также видно, что взятое наугад второе начальное условие не обеспечило выполнение граничного условия при х=1. В целях лучшего выполнения этого граничного условия следует взять большее значение y1(0), например, y1(0)=15, и вновь решить задачу Коши. Соответствующий результат показан на том же Рисунок 10.2 (в центре). Граничное условие выполняется с лучшей точностью, но опять-таки оказалось недостаточным. Для еще одного значения y1(0)=20 получается решение, показанное на Рисунок 10.2 (справа). Из сравнения двух правых графиков легко заключить, что недостающее начальное условие больше 15, но меньше 20. Продолжая подобным образом "пристрелку" по недостающему начальному условию, возможно отыскать правильное решение краевой задачи.
В этом и состоит принцип алгоритма стрельбы. Выбирая пробные начальные условия (проводя пристрелку) и решая соответствующую серию задач Коши, можно найти то решение системы ОДУ, которое (с заданной точностью) удовлетворит граничному условию (или, в общем случае, условиям) на другой границе расчетного интервала.
Конечно, описанный алгоритм несложно запрограммировать самому, оформив его как решение системы заданных алгоритмически уравнений, выражающих граничные условия на второй границе, относительно неизвестных пристрелочных начальных условий. Но делать этого нет необходимости, поскольку он оформлен в Mathcad в виде встроенных функций.
Решение линеаризованной задачи (продолжение листинга 10 9)
Листинг 10.10. Решение линеаризованной задачи (продолжение листинга 10.9)="27.gif">
Последний важный момент, который следует обозначить, связан с решением задач, обладающих значительной нелинейностью. Решение, приведенное в листингах 10.9, 10.10 и на Рисунок 10.12, согласно самой постановке, должно не сильно отличаться от решения линейной краевой задачи, поскольку функции z(x) и z(x) малы по сравнению с Y(x) и у(х). Если же нелинейность сильная, то решение может сильно отличаться от Y и у, и линеаризация (10.11) будет просто неправильной. В этом случае следует слегка усложнить алгоритм решения нелинейной краевой задачи.
Решение краевой задачи
Листинг 10.2. Решение краевой задачи
Первые три строки листинга задают необходимые параметры задачи и саму систему ОДУ. В четвертой строке определяется вектор z. Поскольку правое граничное условие всего одно, то недостающее начальное условие тоже одно, соответственно, и вектор z имеет только один элемент z0. Ему необходимо присвоить начальное значение (мы приняли z0=10, как в листинге 10.1), чтобы запустить алгоритм стрельбы (см. разд. 10.1.2).
Примечание 1
Примечание 1
Начальное значение фактически является параметром численного метода и поэтому может решающим образом повлиять на решение краевой задачи.
В следующей строке листинга векторной функции load(x,z) присваиваются левые граничные условия. Эта функция аналогична векторной переменной, определяющей начальные условия для встроенных функций, решающих задачи Коши. Отличие заключается в записи недостающих условий. Вместо конкретных чисел на соответствующих местах пишутся имена искомых элементов вектора z. В нашем случае вместо второго начального условия стоит аргумент z0 функции load. Первый аргумент функции load — это точка, в которой ставится левое граничное условие. Ее конкретное значение определяется непосредственно в списке аргументов функции sbval. Следующая строка листинга определяет правое граничное условие, для введения которого используется функция score(х,у). Оно записывается точно так же, как система уравнений в функции о. Аргумент х функции score аналогичен функции load и нужен для тех случаев, когда граничное условие явно зависит от координаты х. Вектор score должен состоять из такого же числа элементов, что и вектор z.
Реализованный в функции sbval алгоритм стрельбы ищет недостающие начальные условия таким образом, чтобы решение полученной задачи Коши делало функцию score (х, у) как можно ближе к нулю. Как видно из листинга, результат применения sbval для интервала (0,1) присваивается векторной переменной и. Этот вектор похож на вектор z, только в нем содержатся искомые начальные условия вместо приближенных начальных значений, заданных в z. Вектор и содержит, как и z, всего один элемент и0. С его помощью можно определить решение краевой задачи у (х) (последняя строка листинга). Тем самым функция sbval сводит решение краевых задач к задачам Коши. График решения краевой задачи показан на Рисунок 10.3.
На Рисунок 10.4 показано решение той же самой краевой задачи, но с другим правым граничным условием, соответствующим R=0, т. е. без зеркала на правой границе. В этом случае слабый обратный пучок света образуется исключительно за счет обратного рассеяния излучения от лазера. Конечно, многие из читателей уже обратили внимание, что реальная физическая среда не может создавать такого большого рассеяния назад. Иными словами, более реальны значения r(х)<<а (х). Однако когда коэффициенты в системе ОДУ при разных yi очень сильно (на порядки) различаются, система ОДУ становится жесткой, и функция sbval не может найти решения, выдавая вместо него сообщение об ошибке "Could not find a solution" (Невозможно найти решение).
ВНИМАНИЕ!
Метод стрельбы не годится для решения жестких краевых задач. Поэтому алгоритмы решения жестких ОДУ в Mathcad приходится программировать самому (см. разд. 10.4).
Краевая задача с граничным условием в промежуточной точке
Листинг 10.3. Краевая задача с граничным условием в промежуточной точке
Система уравнений и левое краевое условие вводится так же, как и в предыдущем листинге для функции sbval. Обратите внимание, что таким же образом записано и правое краевое условие. Для того чтобы ввести условие отражения на правой границе, пришлось определить еще один неизвестный пристрелочный параметр z20. Строка листинга, в которой определена функция score, задает условие стрельбы — сшивку двух решений в точке xf. В самой последней строке листинга выдан ответ — определенные численным методом значения обоих пристрелочных параметров, которые объединены в вектор а (мы применили в предпоследней строке операцию транспонирования, чтобы результат получился в форме вектора, а не матрицы-строки). Для корректного построения графика решения лучше составить его из двух частей — решения задачи Коши на интервале (x0,xf) и другой задачи Коши для интервала (xf,x1). Реализация этого способа приведена в листинге 10.4, который является продолжением листинга 10.3. В последней строке листинга 10.4 выведено значение второй искомой функции на правой границе интервала. Всегда полезно проконтролировать, что оно совпадает с соответствующим пристрелочным параметром (выведенным в последней строке листинга 10.3).
Решение краевой задачи (продолжение листинга 10 3)
Листинг 10.4. Решение краевой задачи (продолжение листинга 10.3)
Решение краевой задачи приведено на Рисунок 10.5. С физической точки зрения естественно, что интенсивность света уменьшается быстрее по мере распространения в более плотной среде в правой половине расчетного интервала. В средней точке xf=0.5, как и ожидалось, производные обоих решений имеют разрыв.
Примечание 1
Примечание 1
Еще один пример решения краевых задач с разрывными коэффициентами ОДУ приведен в справочной системе Mathcad.
Решение задачи о собственных колебаниях струны
Листинг 10.5. Решение задачи о собственных колебаниях струны
В первых двух строках листинга определяются функции, входящие в задачу, в том числе р' (х) :=0, и границы расчетного интервала (0,1). В третьей строке дается начальное приближение к собственному значению ?0, в четвертой вводится система ОДУ. Обратите внимание, что она состоит не из двух, а из трех уравнений. Первые два из них определяют эквивалентную (10.3) систему ОДУ первого порядка, а третье необходимо для задания собственного значения в виде еще одного компонента у2 искомого вектора у. Поскольку по определению собственное значение постоянно при всех х, то его производная должна быть приравнена нулю, что отражено в последнем уравнении. Важно также, что во втором из уравнений собственное значение записано как у2, поскольку является одним из неизвестных.
В следующих двух строках листинга задается левое граничное условие, включающее и недостающее условие на собственное значение для третьего уравнения, и правое граничное условие у0=0. В предпоследней строке листинга обычным образом применяется функция sbval, а в последней выводится результат ее работы вместе с известным аналитически собственным значением п2-я2. Как легко убедиться, мы нашли первое собственное значение для n=1, а чтобы найти другие собственные значения, необходимо задать другие начальные приближения к ним (в третьей строке листинга 10.5). Например, выбор ?0=50 приводит ко второму собственному значению 22?2, а ?0=100 — к третьему 32?2.
Чтобы построить график соответствующей собственной функции, надо добавить в листинг строку, программирующую решение задачи Коши, например, такую: u:=rkfixed(load(a,A) ,а,b, 100,D) . Полученные кривые показаны на Рисунок 10.6 в виде коллажа трех графиков, рассчитанных для трех собственных значений.
Примечание 1
Примечание 1
Примеры решения нескольких задач на собственные значения можно найти в различных Ресурсах Mathcad.
Реализация явной разностной схемы
Листинг 10.6. Реализация явной разностной схемы
Дадим минимальные комментарии, надеясь, что заинтересовавшийся читатель с карандашом в руках разберется в порядке индексов и соответствии матричных элементов, а возможно, составит и более удачную программу.
В первой строке листинга определяются функции и константы, входящие в модель, во второй задается число точек сетки N=5 и ее равномерный шаг. Следующие две строки определяют матричные коэффициенты, аппроксимирующие уравнения для Yi; а пятая и шестая — для уi. Седьмая и восьмая строки листинга задают соответственно левое и правое граничное условие, а строки с девятой по одиннадцатую — правые части системы (10.6). В следующей строке завершается построение матрицы А вырезанием из нее левого нулевого столбца. В предпоследней строке листинга применена встроенная функция isolve для решения системы (10.6), а в последней выведены рассчитанные ею неизвестные граничные значения. Графики решения приведены на Рисунок 10.8, причем первые N элементов итогового вектора есть вычисленное излучение вперед, а последние N элементов — излучение назад.
Реализация неявной разностной схемы для жесткой краевой задачи
Листинг 10.7. Реализация неявной разностной схемы для жесткой краевой задачи
Не будем специально останавливаться на обсуждении листинга 10.7, поскольку он почти в точности повторяет предыдущий листинг. Отличие заключается лишь в формировании матрицы А другим способом, согласно неявной схеме. Решение, показанное на Рисунок 10.10, демонстрирует, что произошло "небольшое чудо": "разболтка" исчезла, а распределение интенсивностей стало физически предсказуемым. Обратите внимание, что (из-за взятого нами слишком большого коэффициента ослабления излучения) отраженный пучок света имеет очень маленькую интенсивность, и ее пришлось построить на графике с увеличением в тысячу раз.
Решение нелинейной краевой задачи методом стрельбы
Листинг 10.8. Решение нелинейной краевой задачи методом стрельбы
Следует помнить, что чем сильнее нелинейность и чем жестче краевая задача, тем более сильные требования предъявляются к начальному значению алгоритма, ввод которого осуществляется посредством функции load. Попробуйте повторить расчеты листинга 10.8 с иным сочетанием параметров а0(х) и ?(х), и вы увидите, что с немного более жесткими задачами алгоритм стрельбы уже не справляется.
Примечание 1
Примечание 1
На компакт-диске, прилагаемом к книге, вы найдете еще более интересное решение для другой задачи, с нелинейностью типа (y+Y)2.
Решение линейной (приближенной) краевой задачи
Листинг 10.9.Решение линейной (приближенной) краевой задачи
Модель краевой задачи
Рисунок 10.1. Модель краевой задачи
Полученную задачу называют краевой (boundary value problem), поскольку условия поставлены не на одной, а на обеих границах интервала (0,1). И, в связи с этим, их не решить методами предыдущей главы, предназначенными для задач с начальными условиями. Далее для показа возможностей Mathcad будем использовать этот пример с R=1 и конкретным видом a(x)=const=1 и r(x)=const=0.1, описывающим случай изотропного (не зависящего от координаты х) рассеяния.
Примечание 2
Примечание 2
Модель, представленная на Рисунок 10.1, привела к краевой задаче для системы линейных ОДУ. Она имеет аналитическое решение в виде комбинации экспонент. Более сложные, нелинейные, задачи возможно решить только численно. Нетрудно сообразить, что модель станет нелинейной, если сделать коэффициенты ослабления и рассеяния зависящими от интенсивности излучения. Физически это будет соответствовать изменению оптических свойств среды под действием мощного излучения.
Примечание 3
Примечание 3
Модель встречных световых пучков привела нас к системе уравнений (10.1), в которые входят производные только по одной переменной х. Если бы мы стали рассматривать более сложные эффекты рассеяния в стороны (а не только вперед и назад), то в уравнениях появились бы частные производные по другим пространственным переменным у и z. В этом случае получилась бы краевая задача для уравнений в частных производных, решение которой во много раз сложнее ОДУ.
Нелинейные краевые задачи
10.5. Нелинейные краевые задачиВсе примеры, приведенные пока в этом разделе, связаны с краевыми задачами для линейных ОДУ. Между тем на практике часто приходится исследовать именно нелинейные задачи, которые несравненно сложнее с вычислительной точки зрения. Рассмотрим оба подхода к решению нелинейных задач (алгоритм стрельбы и разностный метод), не забывая о том, что нелинейные ОДУ также могут быть в той или иной степени жесткими.
Неверное решение жесткой краевой задачи по неустойчивой явной разностной схеме
Рисунок 10.9. Неверное решение жесткой краевой задачи по неустойчивой явной разностной схеме
Выходом из положения будет использование неявных разностных схем. Применительно к нашей задаче достаточно заменить правые части уравнений (10.1) значениями не на левой, а на правой границе каждого шага:
Граничные условия, конечно, можно оставить в том же виде (10.5). Поскольку мы имеем дело с линейными дифференциальными уравнениями, то и схему (10.7) легко будет записать в виде матричного равенства (10.6), перегруппировывая соответствующим образом выражение (10.7) и приводя подобные слагаемые. Разумеется, полученная матрица А будет иной, нежели матрица А для явной схемы (10.4). Поэтому и решение (реализация неявной схемы) может отличаться от изображенного на Рисунок 10.9 результата расчетов по явной схеме. Программа, составленная для решения системы (10.7), приведена в листинге 10.7.
О постановке задач
10.1.О постановке задачПостановка краевых задач для ОДУ отличается от задач Коши, рассмотренных в главе 9, тем, что граничные условия для них ставятся не в одной начальной точке, а на обеих границах расчетного интервала. Если имеется система N обыкновенных дифференциальных уравнений первого порядка, то часть из N условий может быть поставлена на одной границе интервала, а оставшиеся условия — на противоположной границе.
Примечание 1
Примечание 1
Дифференциальные уравнения высших порядков можно свести к эквивалентной системе ОДУ первого порядка (см. главу 9).
Чтобы лучше понять, что из себя представляют краевые задачи, рассмотрим их постановочную часть на конкретном физическом примере модели взаимодействия встречных световых пучков. Предположим, что надо определить распределение интенсивности оптического излучения в пространстве между источником (лазером) и зеркалом, заполненном некоторой средой (Рисунок 10.1). Будем считать, что от зеркала отражается R-Я часть падающего излучения (т. е. его коэффициент отражения равен R), а среда как поглощает излучение с коэффициентом ослабления а(х), так и рассеивает его. Причем коэффициент рассеяния назад равен г(х). В этом случае закон изменения интенсивности у0(х) излучения, распространяющегося вправо, и интенсивности y1 (х) излучения влево определяется системой двух ОДУ первого порядка:
Для правильной постановки задачи требуется, помимо уравнений, задать такое же количество граничных условий. Одно из них будет выражать известную интенсивность излучения 10, падающего с левой границы х=0, а второе — закон отражения на его правой границе x=1:
y0(0) = 10; (10.2)
y1(l) = Rxy0(l).
Обыкновенные дифференциальные уравнения краевые задачи
Обыкновенные дифференциальные уравнения: краевые задачиВ этой главе рассматриваются краевые задачи для обыкновенных дифференциальных уравнений (ОДУ). Средства Mathcad, реализующие алгоритм стрельбы (см. разд. 10.2), позволяют решать краевые задачи для систем ОДУ, в которых часть граничных условий поставлена в начальной точке интервала, а остальная часть — в его конечной точке. Также возможно решать уравнения с граничными условиями, поставленными в некоторой внутренней точке.
Краевые задачи во множестве практических приложений часто зависят от некоторого числового параметра. При этом решение существует не для всех его значений, а лишь для счетного их числа. Такие задачи называют задачами на собственные значения (см. разд. 10.3).
Несмотря на то, что, в отличие от задач Коши для ОДУ, в Mathcad не предусмотрены встроенные функции для решения жестких краевых задач, с ними все-таки можно справиться, применив программирование разностных схем, подходящих для решения задач этого класса (см. разд. 104). О подходе к решению нелинейных краевых задач написано в конце главы (см. разд. 10.5).
Первые три собственные функции
Рисунок 10.6. Первые три собственные функции задачи колебаний струны (коллаж трех графиков, рассчитанных листингом 10.5 для разных ?0)
Разностные схемы для ОДУ
10.4. Разностные схемы для ОДУМногие краевые задачи не поддаются решению методом стрельбы. Однако в Mathcad 12 других встроенных алгоритмов нет. Тем не менее это не означает, что по-другому решать краевые задачи невозможно, ведь другие численные алгоритмы несложно запрограммировать самому пользователю. Рассмотрим возможную реализацию наглядного метода, называемого разностным, которым можно решать краевые задачи как для ОДУ, так и для дифференциальных уравнений в частных производных.
Решение краевой задачи для R=0 (продолжение листинга 10 2 при R=0)
Рисунок 10.4. Решение краевой задачи для R=0 (продолжение листинга 10.2 при R=0)
Решение краевой задачи для R=l (продолжение листинга 10 2)
Рисунок 10.3. Решение краевой задачи для R=l (продолжение листинга 10.2)
Решение краевой задачи по неявной разностной схеме (продолжение листинга 10 7)
Рисунок 10.10. Решение краевой задачи по неявной разностной схеме (продолжение листинга 10.7)
Решение краевой задачи разностным методом (продолжение листинга 10 6)
Рисунок 10.8. Решение краевой задачи разностным методом (продолжение листинга 10.6)
Как мы увидели, реализация в Mathcad разностных схем вполне возможна и не слишком трудоемка — предложенная программа состоит всего из двух десятков математических выражений. Конечно, для их написания требуется и время, и часто кропотливые расчеты, но, собственно, в этом и состоит работа математика. Кстати говоря, при небольшом числе шагов расчеты по разностным схемам не требуют существенного времени (программа, приведенная в листинге 10.6, работает быстрее, чем метод стрельбы, встроенный в функцию sbval). Существуют, кроме того, весьма очевидные для многих читателей пути ускорения расчетов, связанные с применением более подходящих методов решения систем линейных уравнений с разреженной матрицей.
Решение краевой задачи с разрывом в xf=0 5 (продолжение листингов 10 3—10 4)
Рисунок 10.5. Решение краевой задачи с разрывом в xf=0.5 (продолжение листингов 10.3—10.4)
Ради справедливости необходимо заметить, что разобранную краевую задачу легко решить и с помощью функции sbval, заменив в листинге 10.2 зависимость а (х) на третью строку листинга 10.3. В этом случае листинг 10.2 даст в точности тот же ответ, что показан на Рисунок 10.5. Однако в определенных случаях (в том числе из соображений быстроты расчетов) удобнее использовать функцию bvaifit, т. е. вести пристрелку с обеих границ интервала.
СОВЕТ
Если вы имеете дело с подобными уравнениями, попробуйте сначала решить их как обычную краевую задачу с помощью более надежной и легкой в применении функции sbval.
Решение краевых задач средствами Mathcad
10.2. Решение краевых задач средствами MathcadДля решения краевых задач в Mathcad реализован наиболее популярный алгоритм, называемый методом стрельбы или пристрелки (shooting method). Он, по сути, сводит решение краевой задачи к решению серии задач Коши с различными начальными условиями. Рассмотрим сначала основной принцип алгоритма стрельбы на примере модели световых пучков (см. Рисунок 10.1), а встроенные функции, реализующие этот алгоритм, приведем в следующем разделе.
Решение нелинейной краевой задачи (продолжение листинга 10 8)
Рисунок 10.11. Решение нелинейной краевой задачи (продолжение листинга 10.8)
Решение нелинейной краевой задачи разностным методом (продолжение листингов 10 9 и 10 10)
Рисунок 10.12. Решение нелинейной краевой задачи разностным методом (продолжение листингов 10.9 и 10.10)
Обозначим полученное в результате решение, как и в листинге 10.10, вектором J(?), подчеркивая тем самым его зависимость от параметра нелинейности. Очевидно, что J(0) есть решение линейной задачи. Для того чтобы решить задачу с сильной нелинейностью, т. е. довольно большим ?=?1, можно организовать продолжение по ? как по параметру. Иными словами, используя в качестве начального приближения J (0), можно решить задачу для другого, малого ?=Д?, получив J(A?), затем, взяв это J(?) в качестве приближенного решения, получить J(2A?) и т. д. малыми шагами добраться до желаемого ?1.
Примечание 1
Примечание 1
Упрощенную реализацию этого алгоритма вы найдете на компакт-диске, прилагаемом к книге. Она связана с выводом во внешний файл данных результата задачи из листинга 10.10 и считыванием из него же этих данных в качестве входной информации для следующей итерации. В качестве нулевой итерации используется решение линейной задачи, выводимое предварительно в файл из усовершенствованного листинга 10.9.
Сетка покрывающая расчетный интервал
Рисунок 10.7. Сетка, покрывающая расчетный интервал
Примечание 1
Примечание 1
Существует множество способов аппроксимации дифференциальных уравнений разностными. От выбора конкретного варианта зависит не только простота, быстрота и удобство вычислений, но и сама возможность получения правильного ответа.
Получилась система (по числу шагов) 2(N-1) разностных линейных алгебраических уравнений с 2N неизвестными Yi и уi. Для того чтобы она имела единственное решение, надо дополнить число уравнений до 2N.
Это можно сделать, записав в разностном виде оба граничных условия: \
Y0=10, YN=RYN. (10.5)
Сформированная полная система алгебраических уравнений называется разностной схемой, аппроксимирующей исходную краевую задачу. Обратите внимание, что правые части разностных уравнений системы (10.4) на каждом шаге записаны для левой границы шага. Такие разностные схемы называют явными, т. к. все значения Yi+1 и yi+1 находятся в левой части уравнений. Полученную явную разностную схему легко записать в матричной форме
Аz=В, (10.6)
где z — неизвестный вектор, получающийся объединением векторов Y и у. Решив систему (10.6), мы получим решение краевой задачи.
Примечание 2
Примечание 2
На самом деле все несколько сложнее, поскольку, вообще говоря, необходимо еще доказать, что, во-первых, разностная схема действительно аппроксимирует дифференциальные уравнения и, во-вторых, при N->? разностное решение действительно сходится к дифференциальному.
Процесс решения системы разностных уравнений называют также реализацией разностной схемы. Программа, которая решает рассматриваемую краевую задачу разностным методом, приведена в листинге 10.6.
Задачи на собственные значения для ОДУ
10.3. Задачи на собственные значения для ОДУЗадачи на собственные значения — это краевые задачи для системы ОДУ, в которой правые части зависят от одного или нескольких параметров X. Значения этих параметров неизвестны, а решение краевой задачи существует только при определенных ?k, которые называются собственными значениями (eigenvalues) задачи. Решения, соответствующие этим ?k, называют собственными функциями (eigenflmctions) задачи. Правильная постановка таких задач требует формулировки количества граничных условий, равного сумме числа уравнений и числа собственных значений. Физическими примерами задач на собственные значения являются, например, уравнение колебаний струны, уравнение Шредингера в квантовой механике, уравнения волн в резонаторах и многие другие.
С вычислительной точки зрения, задачи на собственные значения очень похожи на рассмотренные выше краевые задачи. В частности, для многих из них также применим метод стрельбы (см. разд. 10.1.2). Отличие заключается в пристрелке не только по недостающим левым граничным условиям, но еще и по искомым собственным значениям. В Mathcad для решения задач на собственные значения используются те же функции sbval и bvaifit. В их первый аргумент, т. е. вектор, присваивающий начальные значения недостающим начальным условиям, следует включить и начальное приближение для собственного значения.
Рассмотрим методику решения на конкретном примере определения собственных упругих колебаний струны. Профиль колебаний струны у(х) описывается линейным дифференциальным уравнением второго порядка
где р(х) и q(x) — жесткость и плотность, которые, вообще говоря, могут меняться вдоль струны. Если струна закреплена на обоих концах, то граничные условия задаются в виде у(0)=у(1)=0. Сформулированная задача является частным случаем задачи Штурма—Лиувилля. Поскольку решается система двух ОДУ, содержащая одно собственное значение x, то по идее задача требует задания трех (2+1) условий. Однако, как легко убедиться, уравнение колебаний струны — линейное и однородное, поэтому в любом случае решение у(х) будет определено с точностью до множителя. Это означает, что производную решения можно задать произвольно, например, у' (0)=1, что и будет третьим условием. Тогда краевую задачу можно решать как задачу Коши, а пристрелку вести только по одному параметру — собственному значению.
Процедура поиска первого собственного значения представлена в листинге 10.5.
Матричные вычисления в Mathcad
Явная схема Эйлера
11.2.1. Явная схема ЭйлераРассмотрим сначала математические аспекты построения разностной схемы для уравнения диффузии тепла, а затем приведем примеры работы разработанного алгоритма применительно к линейному и нелинейному уравнениям.
Построение разностной схемы
Используем для решения уравнения теплопроводности шаблон, изображенный на Рисунок 11.6. Для дискретизации второй производной по пространственной координате необходимо использовать три последовательных узла, в то время как для разностной записи первой производной по времени достаточно двух узлов. Записывая на основании данного шаблона дискретное представление для (i, k) -го узла, получим разностное уравнение:
Приведем в разностной схеме (11.8) подобные слагаемые, перенеся в правую часть значения сеточной функции с индексом k (как часто говорят, с предыдущего слоя по времени), а в левую — с индексом k+i (т. е. со следующего временного слоя). Кроме этого, введем коэффициент с, который будет характеризовать отношение шагов разностной схемы по времени и пространству
Несколько забегая вперед, заметим, что значение
параметра с, называемого коэффициентом Куранта, имеет большое значение для анализа устойчивости разностной схемы. С учетом этих замечаний, разностная схема (11.8) запишется в виде:
Классификация уравнений в частных производных
11.1.1. Классификация уравнений в частных производныхПостановка задач для уравнений в частных производных включает определение самого уравнения (или системы нескольких уравнений), а также необходимого количества краевых условий (число и характер задания которых определяются спецификой уравнения). По своему названию уравнения должны содержать частные производные неизвестной функции и (или нескольких функций, если уравнений несколько) по различным аргументам, например, пространственной переменной х и времени t. Соответственно, для решения задачи требуется вычислить функцию нескольких переменных, например, u(x,t) в некоторой области определения аргументов 0
Сами уравнения в частных производных (несколько условно) можно разделить на три основных типа:
Некоторые более сложные уравнения нельзя однозначно подогнать под приведенную классификацию, тогда говорят о гибридных типах уравнений.
Параболические и гиперболические уравнения
11.3.1. Параболические и гиперболические уравненияРазработчики впервые применили дополнительные встроенные функции для решения параболических и гиперболических уравнений в частных производных в версии Mathcad 11, отлично осознавая значимость этих задач для современного исследователя и инженера. Предусмотрены два варианта решения: при помощи вычислительного блока Given/pdesolve, а также при помощи встроенной функции numol. Первый путь проще в применении и нагляднее, зато второй позволяет автоматизировать процесс решения уравнений в частных производных, например, если нужно включить его в качестве составного шага в более сложную Mathcad-программу.
Вычислительный блок Given /pdesolve
Встроенная функция pdesolve применяется в рамках вычислительного блока, начинающегося ключевым словом Given, и пригодна для решения различных гиперболических и параболических уравнений. Она предназначена для решения одномерного уравнения (или системы уравнений) в частных производных (того, которое определит пользователь в рамках вычислительного блока Given), зависящего от времени t и пространственной координаты х, имеет целый набор различных аргументов и работает следующим образом:
В качестве примера использования функции pdesolve (листинг 11.4) используем то же самое одномерное уравнение теплопроводности (11.5) с граничными и начальными условиями (11.6) и (11.7).
Эллиптические уравнения
11.3.2. Эллиптические уравненияРешение эллиптических уравнений в частных производных реализовано только для единственного типа задач — двумерного уравнения Пуассона.
Это уравнение содержит вторые производные функции u(х,у) по двум пространственным переменным;
Уравнение Пуассона описывает, например, распределение электростатического поля u(х,у) в двумерной области с плотностью заряда f (x,y), или (см. разд. 11.1.2) стационарное распределение температуры u(х,у) на плоскости, в которой имеются источники (или поглотители) тепла с интенсивностью f (х,у) .
Примечание 1
Примечание 1
Несмотря на то, что применение встроенных функций, описанных в данном разделе, анонсировано разработчиками Mathcad только для уравнения Пуассона, их можно применять и для решения других уравнений, даже необязательно эллиптического типа. О том, как осуществить такие расчеты, написано в конце данного раздела.
Уравнение Пуассона с нулевыми граничными условиями
Корректная постановка краевой задачи для уравнения Пуассона требует задания граничных условий. В Mathcad решение ищется на плоской квадратной области, состоящей из (M+1)х(M+1) точек. Поэтому граничные условия должны быть определены пользователем для всех четырех сторон упомянутого квадрата. Самый простой вариант — это нулевые граничные условия, т. е. постоянная температура по всему периметру расчетной области. В таком случае можно использовать встроенную функцию multigrid:
ВНИМАНИЕ!
Сторона квадрата расчетной области должна включать точно N=2n шагов, т. е. 2n+1 узлов, где n — целое число.
Параметр численного метода ncycle в большинстве случаев достаточно взять равным 2.
Неявная схема Эйлера
11.2.2. Неявная схема ЭйлераВ отличие от явной схемы Эйлера, неявная является безусловно-устойчивой (т. е. не выдающей "разболтки" ни при каких значениях коэффициента Куранта). Однако ценой устойчивости является необходимость решения на каждом шаге по времени системы алгебраических уравнений.
Построение неявной разностной схемы
Чтобы построить неявную разностную схему для уравнения диффузии, используем шаблон, изображенный на Рисунок 11.11, т. е. для дискретизации пространственной производной будем брать значения сеточной функции с верхнего (неизвестного) слоя по времени. Таким образом, разностное уравнение для (i,k)-ro узла будет отличаться от уравнения для явной схемы (11.8) только индексами по временной координате в правой части:
Если привести подобные слагаемые, то получится система уравнений, связывающая для каждого 1-го узла три неизвестных значения сеточной функции (в самом этом узле и в соседних с ним слева и справа узлах). Множители при неизвестных значениях сеточной функции в узах шаблона показаны на Рисунок 11.11 в виде подписей, подобно тому, как это было сделано для явной схемы (см. Рисунок 11.6).
Пример уравнение диффузии тепла
11.1.2. Пример: уравнение диффузии теплаНа протяжении всей главы мы будем использовать в качестве примера очень наглядное и имеющее различные, от очевидных до самых неожиданных, решения уравнение теплопроводности.
Двумерное динамическое уравнение
Рассмотрим следующее параболическое уравнение в частных производных, зависящее от трех переменных — двух пространственных х и у, а также от времени t:
Примечание 1
Примечание 1
Выражение в скобках в правой части уравнения (сумму вторых пространственных производных функции и часто, ради краткости, обозначают при помощи оператора Лапласа: ?u).
Это уравнение называется двумерным уравнением теплопроводности или, по-другому, уравнением диффузии тепла. Оно описывает динамику распределения температуры u(x,y,t) на плоской поверхности (например, на металлической пластине) в зависимости от времени (Рисунок 11.1). Физический смысл коэффициента о, который, вообще говоря, может быть функцией как координат, так и самой температуры заключается в задания скорости перетекания тепла от более нагретых областей в менее нагретые. Функция ?(x,y,t,u) описывает приток тепла извне, т. е. источники тепла, которые также могут зависеть как от пространственных координат (что задает локализацию источников), так и от времени и от температуры и.
О возможности решения многомерных уравнений
11.2.3.О возможности решения многомерных уравненийВсе, что было сказано до сих пор, касалось исключительно способов решения одномерных (в смысле пространственных координат) уравнений. И алгоритмы разностных схем, и встроенные функции, включая появившиеся в 11-й версии (см. следующий разд.), относились к уравнениям, зависящим от одной пространственной координаты.
Можно ли при помощи Mathcad решать двумерные или трехмерные (пространственные) уравнения? С точки зрения программирования пользователем численных алгоритмов типа метода сеток принципиальных ограничений нет. Разумеется, если сначала аккуратно выписать разностную схему соответствующего многомерного дифференциального уравнения, то вполне возможно запрограммировать ее при помощи описанных нами средств. Самым главным противодействием будет существенное увеличение времени расчетов. Простая оценка необходимого количества операций показывает, что если ввод зависимости уравнения от второй пространственной координаты многократно увеличивает число разностных уравнений, которые должны решаться при реализации каждого шага по времени. К примеру, если используется пространственная сетка из 100 узлов по каждой координате, то вместо 102 разностных уравнений на каждом шаге придется решать уже 104 уравнений, т. е. объем вычислений сразу же возрастает в 100 раз. Вообще говоря, пакет Mathcad не является экономичной средой вычислений, и бороться с их сильно возрастающим объемом следует пользователю еще на этапе разработки алгоритма. Хорошим примером такой борьбы может служить применение специфических алгоритмов, типа метода прогонки (см. разд. "Алгоритм прогонки"этой главы).
Приведем некоторые дополнительные замечания, связанные с возможностью осуществить редукцию громоздких (в смысле организации вычислений) двумерных задач к более простым. Рассмотрим ради примера двумерное уравнение теплопроводности (11.1) без источника с нулевыми граничными условиями и некоторым начальным двумерным распределением температуры по расчетной поверхности:
Произведем дискретизацию данного уравнения по временной координате, заменяя первую производную ее разностным аналогом и несколько перегруппировывая слагаемые и множители:
Как вы видите, мы используем неявную разностную схему, заранее заботясь о том, чтобы разностная задача была более устойчивой. Здесь ui(x,y) — известная с предыдущего шага по времени функция двух пространственных переменных, а и1+1(х,у) — функция, подлежащая определению при реализации каждого шага по времени.
Посмотрим на полученную задачу с несколько другой стороны — а именно как на дифференциальное уравнение относительно неизвестной функции двух переменных ui+1(x,y). Подчеркнем, что такое уравнение получается для каждого шага по времени, т. е. для реализации всей разностной схемы требуется решить большое число таких уравнений.
С предложенной точки зрения, на каждом временном шаге необходимо решить некоторое двумерное эллиптическое линейное уравнение, причем его граничные условия определяются граничными условиями исходной задачи. Это уравнение очень похоже на уравнение Пуассона с той лишь разницей, что в его правую часть, описывающую источник, входит неизвестная функция (к счастью, линейно). Таким образом, зависимость от найденного на предыдущем шаге по времени решения определяется зависимостью от него источника, т. е. правой части.
Суммируя сказанное, можно констатировать, что если вы имеете запрограммированный алгоритм решения выписанного эллиптического уравнения, чуть более сложного, чем уравнение Пуассона, то его с легкостью можно использовать в качестве подпрограммы реализации разностной схемы двумерного уравнения теплопроводности. Забегая вперед, приходится отметить, что, к сожалению, встроенные функции Mathcad для решения уравнения Пуассона в данном случае не годятся в качестве такой подпрограммы, поскольку предполагают независимость источника от самой неизвестной функции, и могут справляться лишь с правой частью, зависящей только от пространственных координат (см. разд. 11.3.2).
Не будем далее останавливаться на способах решения многомерных уравнений, ограничившись этими замечаниями относительно путей оптимизации алгоритмов их решения.
Численное решение обратного уравнения
Рисунок 11.5. Численное решение обратного уравнения теплопроводности дает совершенно нефизичную картину динамики температуры (см. листинг 11.2 ниже с параметром D=-1)
Численное решение уравнения теплопроводности
Рисунок 11.10. Численное решение уравнения теплопроводности при помощи явной схемы Эйлера (см. листинг 11.1 ниже с временным шагом t=0.0015)
Вместо ожидаемого решения получаются совершенно неожиданные профили температуры, которые быстро осциллируют вдоль пространственной координаты, причем амплитуда и число пиков этих осцилляции быстро увеличиваются от шага к шагу. Совершенно ясно, что полученное решение не имеет ничего общего с физикой моделируемого явления, а является следствием внутренних свойств самой разностной схемы, которые до этого были для нас скрыты.
Характерная "разболтка" решения как раз и является проявлением неустойчивости явной схемы Эйлера для выбранного соотношения шагов по времени и пространству. В теории численных методов показывается, что явная схема Эйлера для уравнения теплопроводности устойчива при значениях коэффициента Куранта, меньших 1, и неустойчива в противоположном случае. Иными словами, существует ограничение для выбора соотношения шагов, заключающееся в том, что для расчета на более частых пространственных сетках необходимо использовать также и малые шаги по времени.
Примечание 4
Примечание 4
Как несложно убедиться, для t=0.0005 коэффициент Куранта C=0.4, для t=0.0010 он все еще меньше единицы: C=0.8, а для t=0. 0015 решение уже больше единицы: C=1.2, в связи с чем схема становится неустойчивой (см. Рисунок 11.10).
Дифференциальные уравнения в частных производных
Дифференциальные уравнения в частных производныхДифференциальные уравнения в частных производных представляют собой одну из наиболее сложных и одновременно интересных задач вычислительной математики. Эти уравнения характеризуются тем, что для их решения не существует единого универсального алгоритма, и большинство задач требует своего собственного особого подхода. Уравнениями в частных производных описывается множество разнообразных физических явлений, и с их помощью можно с успехом моделировать самые сложные явления и процессы (диффузия, гидродинамика, квантовая механика, экология и т. д.).
Дифференциальные уравнения в частных производных требуют нахождения функции не одной, как для ОДУ, а нескольких переменных, например, f (х,у) или f(x,t). Постановка задач (см. разд. 11.Т) включает в себя само уравнение (или систему уравнений), содержащее производные неизвестной функции по различным переменным (частные производные), а также определенное количество краевых условий на границах расчетной области.
Несмотря на то, что Mathcad обладает довольно ограниченными возможностями по отношению к уравнениям в частных производных, в нем имеется несколько встроенных функций (см. разд. 11.3). Решать уравнения в частных производных можно и путем непосредственного программирования пользовательских алгоритмов (см. разд. 11.2). Автор совершенно сознательно сначала рассматривает численные методы для решения уравнений в частных производных, а уже затем описывает предназначенные для этого встроенные функции, чтобы читатель ясно осознавал, каким образом Mathcad производит расчеты. "Слепое" использование встроенных функций для решения уравнений в частных производных не всегда бывает успешным, и ответственность за верный выбор их параметров часто ложится на пользователя, которому необходимо четко представлять основные принципы функционирования численных алгоритмов, примененных во встроенных функциях.
Физическая модель одномерного уравнения теплопроводности
Рисунок 11.3. Физическая модель одномерного уравнения теплопроводности
Физическая модель описываемая двумерным уравнением теплопроводности
Рисунок 11.1. Физическая модель, описываемая двумерным уравнением теплопроводности
Для того чтобы правильно поставить краевую задачу для двумерного уравнения теплопроводности, следует определить следующие дополнительные условия:
Примечание 2
Примечание 2
Если рассматривается не одно уравнение в частных производных, а система уравнений, то соответствующие начальные и граничные условия должны быть поставлены для каждой из неизвестных функций.
Стационарное двумерное уравнение
Частный случай уравнения теплопроводности определяет стационарную, т. е. не зависящую от времени, задачу. Стационарное уравнение описывает физическую картину распределения температуры по пластине, не изменяющуюся с течением времени. Такая картина может возникнуть при условии, что стационарный источник тепла действует довольно продолжительное время, и переходные процессы, вызванные его включением, прекратились. Пример численного решения такого уравнения показан на Рисунок 11.2 в виде поверхности u(х,у).
Дифференциальные уравнения в частных производных
Глава 11. Дифференциальные уравнения в частных производных| Содержание |
График линий уровня решения уравнения Пуассона (продолжение листинга 11 7)
Рисунок 11.18. График линий уровня решения уравнения Пуассона (продолжение листинга 11.7)
Уравнение Пуассона с произвольными граничными условиями
В более сложных случаях, например, для решения краевой задачи с ненулевыми условиями на границах, следует использовать другую встроенную функцию relax, имеющуюся в Mathcad:
Параметр численного алгоритма характеризует скорость сходимости итераций. Он должен быть числом от о до 1. В матрице граничных условий v необходимо задать только граничные элементы, исходя из значения краевых условий по периметру расчетной области. Прочие (внутренние) элементы этой матрицы служат для задания начального приближения к решению.
Суть алгоритма релаксации сводится к тому, что в ходе итераций происходит проверка уравнений и соответствующая коррекция значений искомой функции в каждой точке. Если начальное приближение выбрано удачно, то можно надеяться, что алгоритм сойдется ("релаксирует") к правильному решению.
ВНИМАНИЕ!
Все матрицы, задающие как коэффициенты разностной схемы а, b, с, d, e, граничные условия v, так и само решение F, должны иметь одинаковый размер (M+1)х(M+1), соответствующий размеру расчетной области. При этом целое число м обязательно должно быть степенью двойки: м=2п.
Решение уравнения Пуассона с тремя источниками разной интенсивности при помощи функции relax приведено в листинге 11.8.
График поверхности решения уравнения Пуассона (продолжение листинга 11 7)
Рисунок 11.17. График поверхности решения уравнения Пуассона (продолжение листинга 11.7)
Явная схема для линейного уравнения теплопроводности
Листинг 11.1. Явная схема для линейного уравнения теплопроводности
Начальное распределение температуры вдоль расчетной области и решение для двух моментов времени показано на Рисунок 11.7 сплошной, пунктирной и штриховой линиями соответственно. Физически такое поведение вполне естественно — с течением времени тепло из более нагретой области перетекает в менее нагретую, а зона изначально высокой температуры остывает и размывается.
Неявная схема для линейного уравнения теплопроводности
Листинг 11.2. Неявная схема для линейного уравнения теплопроводности
Алгоритм прогонки
Приведем в данном разделе описание чрезвычайно популярного алгоритма реализации неявных разностных схем, который называется методом прогонки. Этот алгоритм имел историческое значение для становления технологий расчетов уравнений в частных производных, и мы просто не можем не упомянуть о нем в этой книге.
Примечание 2
Примечание 2
Сразу оговоримся, что его применение для решения уравнений в частных производных в среде Mathcad может быть оправдано, только если вы работаете с очень частыми сетками, которые приводят к системам разностных уравнений большой размерности и, соответственно, очень большому времени вычислений.
Основным вычислительным ядром программы, реализующей на Mathcad неявную разностную схему, было решение (на каждом временном слое) системы линейных алгебраических уравнений, задаваемых матрицей А. Заметим, что эта матрица, как говорят, имеет диагональное преобладание, а точнее, является трехдиагоналъной (Рисунок 11.13). Все ее элементы, кроме элементов на главной диагонали и двух соседних диагоналях, равны нулю. С точки зрения оптимизации быстродействия алгоритма применение встроенной функции isolve является весьма расточительным, поскольку основной объем арифметических операций, выполняемых компьютером (а он составляет, как нетрудно убедиться, величину порядка м2), сводится к непроизводительному перемножению нулей.
Алгоритм прогонки (продолжение листинга 11 2)
Листинг 11.3. Алгоритм прогонки (продолжение листинга 11.2)
Решение одномерного уравнения теплопроводности
Листинг 11.4. Решение одномерного уравнения теплопроводности
Для корректного использования функции pdesolve предварительно, после ключевого слова Given, следует записать само уравнение и граничные условия при помощи логических операторов (для их ввода в Mathcad существует специальная панель). Обратите внимание, что уравнение должно содержать имя неизвестной функции u(x,t) вместе с именами аргументов (а не так, как она записывается в пределах встроенной функции pdesolve). Для идентификации частных производных в пределах вычислительного блока следует использовать нижние индексы, например, uxx(,t), для обозначения второй производной функции и по пространственной координате х.
Как видно из Рисунок 11.14, на котором изображены результаты расчетов по листингу 11.4, встроенная функция с успехом справляется с уравнением диффузии, отыскивая уже хорошо знакомое нам решение. Заметим, что использование встроенной функции pdesolve связано с довольно громоздкими вычислениями, которые могут отнимать существенное время.
Примечание 1
Примечание 1
Как вы можете заметить, выбирать величину шага по пространственной и временной переменным может как сам алгоритм, так и пользователь (неявным образом, через число узлов сетки). Читателю предлагается повторить вычисления листинга 11.4 для различных комбинаций параметров (главным образом, числа узлов сетки), чтобы проверить, в каких случаях алгоритм встроенной функции справляется с задачей, выдавая верное решение, а в каких дает сбой.
Решение волнового уравнения
Листинг 11.5. Решение волнового уравнения
Решение волнового уравнения при помощи функции numol
Листинг 11.6. Решение волнового уравнения при помощи функции numol
Как вы видите, функция numol имеет еще большее число аргументов, нежели pdesolve, и позволяет автоматизировать применение метода сеток. Однако пользоваться ею намного сложнее, чем вычислительным блоком, поскольку и уравнения, и начальные и граничные условия должны быть записаны в специальном формате. Применение функции numol оправданно, когда необходимо включить решение уравнений в частных производных в более сложные вычисления в качестве подпрограммы, организовать серию расчетов с меняющимся параметром, подготовить анимацию графиков решения и т. п.
Решение уравнения Пуассона с нулевыми граничными условиями
Листинг 11.7.Решение уравнения Пуассона с нулевыми граничными условиями
В первой строке листинга задается значение м=32, в двух следующих строках создается матрица правой части уравнения Пуассона, состоящая из всех нулевых элементов, за исключением одного, задающего расположение источника. В последней строке матрице с присваивается результат действия функции multigrid. Обратите внимание, первый ее аргумент сопровождается знаком "минус", что соответствует записи правой части уравнения Пуассона (11.11). Графики решения показаны на Рисунок 11.17 и 11.18 в виде трехмерной поверхности и линий уровня соответственно.
содержит пример использования
Листинг 11.7 содержит пример использования функции multigrid для расчета краевой задачи на области 33х33 точки и точечным источником тепла в месте, задаваемом координатами (15,20) внутри этой области.Решение уравнения Пуассона с помощью функции relax
Листинг 11.8. Решение уравнения Пуассона с помощью функции relax
Первые три строки имеют тот же смысл, что и в предыдущем листинге. Только вместо одного источника тепла взято их другое распределение — один сильный источник, один более слабый и один сток тепла. В следующих шести строках задаются коэффициенты разностной схемы. Отложим их обсуждение до последнего раздела этой главы, ограничившись утверждением, что для решения уравнения Пуассона коэффициенты должны быть взяты именно такими, как показано в листинге 11.8. В предпоследней строке задана матрица нулевых граничных условий и нулевых начальных приближений, а в последней матрице с присваивается результат действия функции relax. График полученного решения в виде линий уровня показан на Рисунок 11.19.
Решение уравнения теплопроводности при помощи функции relax
Листинг 11.9. Решение уравнения теплопроводности при помощи функции relax
Результат действия программы листинга 11.9 показан на Рисунок 11.21 в виде трехмерной поверхности. Если сравнить Рисунок 11.21 с Рисунок 11.4, полученным при расчетах по запрограммированной разностной схеме, то в графиках Рисунок 11.4 нетрудно узнать сечения этой поверхности плоскостями t=const. Еще раз подчеркнем, что использовать встроенную функцию можно только для тех уравнений, которые допускают построение разностной схемы типа "крест" (см. Рисунок 11.17) или составного фрагмента этой схемы.
Матрица системы линейных разностных уравнений для неявной схемы (листинг 11 2 для M=10)
Рисунок 11.13. Матрица системы линейных разностных уравнений для неявной схемы (листинг 11.2 для M=10)
Для отыскания решения линейных систем алгебраических уравнений имеется чрезвычайно эффективный алгоритм, называемый прогонкой, который позволяет уменьшить число арифметических операций на целый порядок, т. е. до значения порядка м. Это означает, что при использовании пространственных сеток с юоо узлами выигрыш во времени вычислений составит величину порядка ю3! Реализация данного алгоритма приведена в листинге 11.3, который является продолжением листинга 11.2, используя определенные в нем коэффициенты матрицы А, а также начальное условие.
Программа листинга 11.3 осуществляет пересчет одного шага по времени, т. е. заменяет содержимое столбца и с предыдущего временного слоя вычисленными значениями неизвестной функции со следующего слоя. Первые пять строк листинга 11.3 представляют так называемый обратный ход прогонки, а последние две строки — ее прямой ход. Заинтересовавшемуся читателю предлагается самому оформить представленный алгоритм прогонки в виде программы решения разностных уравнений для вычисления произвольного временного слоя по примеру листингов 11.1 и 11.2. Заметим, что описание этого знаменитого алгоритма можно отыскать практически в любом современном учебнике по численным методам.
О постановке задач
11.1.О постановке задачОстановимся сначала на общей классификации уравнений в частных производных, принятой в математике, а затем более детально поговорим о постановочной части соответствующих задач на примере различных вариаций уравнения диффузии тепла, которое допускает наглядную физическую интерпретацию.
Разностные схемы
11.2. Разностные схемыРассмотрим одномерное уравнение теплопроводности (11.3) и на его примере разберем наиболее часто использующийся для численного решения уравнений в частных производных метод сеток. Выпишем еще раз само уравнение
а также и начальное
u(x,0) = u0(x), (11.6)
и граничные условия
u(0, t) = f0(t) , u(L, t) = f1(t), (11.7)
которые необходимы для правильной с математической точки зрения постановки задачи.
Основная идея численного решения уравнений в частных производных очень похожа на метод решения краевых задач для ОДУ, рассмотренный нами в предыдущей главе. Основным отличием от ОДУ является необходимость дискретизации уравнения не по одной, а по нескольким переменным (в зависимости от размерности задачи).
Таким образом, сначала следует покрыть расчетную область (x,t) сеткой и использовать затем узлы этой сетки для разностной аппроксимации уравнения. В результате, вместо поиска непрерывных зависимостей u(x,t) достаточно будет отыскать значения функции в узлах сетки (а ее поведение в промежутках между узлами может быть получено при помощи построения какой-либо интерполяции). По этой причине дискретное представление функции и часто называют сеточной функцией.
Поскольку уравнения в частных производных по определению зависят от производных неизвестных функций по нескольким переменным, то способов дискретизации этих уравнений может быть, как правило, несколько. Конфигурацию узлов, используемую для разностной записи уравнений в частных производных на сетке, называй шаблоном, а полученную систему разностных уравнений — разностной схемой. О принципах построения разностных схем и, в частности, о классах явных и неявных схем мы уже подробно говорили на примере краевых задач для ОДУ (см. разд. 10.4.1), поэтому, излишне не повторяясь, перейдем к рассмотрению типичных особенностей уравнений в частных производных, которые возникают при разработке и реализации разностных схем.
Решение линейного уравнения теплопроводности (продолжение листинга 11 1)
Рисунок 11.7. Решение линейного уравнения теплопроводности (продолжение листинга 11.1)
Примечание 2
Примечание 2
Несколько забегая вперед, заметим, что показанное на Рисунок 11.7 решение и соответствует коэффициенту Куранта С=0.4. Попробуйте осуществить расчет с увеличенным временным шагом, чтобы коэффициент с был больше 1, и посмотрите, что из этого получится (такой расчет и его объяснение приведены ниже в разд. "Устойчивость").
Нелинейное уравнение
Намного более интересные решения можно получить для нелинейного уравнения теплопроводности, например, с нелинейным источником тепла ф(u)=103(u-u3). Заметим, что в листинге 11.1 мы предусмотрительно определили коэффициент диффузии и источник тепла в виде пользовательских функций, зависящих от аргумента и, т. е. от температуры. Если бы мы собирались моделировать явную зависимость их от координат, то следовало бы ввести в пользовательскую функцию в качестве аргумента переменную х, как это сделано для источника тепла ф. Поэтому нет ничего проще замены определения этих функций с констант D(U)=1 и ф(х,u)=о на новые функции, которые станут описывать другие модели диффузии тепла. Начнем с того, что поменяем четвертую строку листинга 11.1 на ф(х,и)=ю3-(и-и3), не изменяя пока постоянного значения коэффициента диффузии.
Примечание 3
Примечание 3
С физической точки зрения, зависимость коэффициента диффузии и функции источника тепла от температуры означает, что эти параметры будут меняться от точки к точке среды, определяясь локальными значениями текущей температуры в этих точках. Ввод ненулевого источника тепла означает, что среда получает определенное количество тепла, тем большее, чем больше локальная температура. Можно догадаться, что введение такой зависимости может моделировать, в частности, горение среды.
Если осуществить расчеты с упомянутым источником (имеющим кубическую нелинейность), то получится очень интересное решение уравнения теплопроводности, имеющее профиль тепловых фронтов. С течением времени граница раздела высокой и низкой температуры распространяется в обе стороны от зоны первичного нагрева, оставаясь весьма четко выделенной (Рисунок 11.8).
Еще более неожиданные решения возможны при нелинейности также и коэффициента диффузии. Например, если взять квадратичный коэффициент диффузии D(x,u)=u2 (что с учетом его умножения на неизвестные функции создаст кубическую нелинейность уравнения), а также ф(х,u)=103-u3-5, то вы сможете наблюдать совсем иной режим горения среды. В отличие от рассмотренного эффекта распространения тепловых фронтов, горение оказывается локализованным в области первичного нагрева среды, причем температура в центре нагрева со временем возрастает до бесконечной величины (Рисунок 11.9). Такое решение описывает так называемый режим горения "с обострением".
Решение линейного
Рисунок 11.12. Решение линейного уравнения теплопроводности при помощи неявной схемы на первом слое по времени (листинг 11.2)
Решение одномерного уравнения теплопроводности (см листинг 11 1 ниже)
Рисунок 11.4. Решение одномерного уравнения теплопроводности (см. листинг 11.1 ниже)
Линейное и нелинейное уравнения
Если присмотреться к уравнению диффузии тепла внимательнее, то можно условно разделить практические случаи его использования на два типа.
Решения линейных уравнений в частных производных, как правило, получаются вполне предсказуемыми, и их часто можно получить аналитически (этим проблемам посвящены соответствующие разделы науки, называемой математической физикой). В случае уравнения теплопроводности линейная задача описывает физически ожидаемое решение, выражающее остывание пластины или стержня в форме перетекания тепла от нагретого центра к холодной периферии.
Нелинейные уравнения, напротив, могут демонстрировать самые неожиданные решения, причем в подавляющем большинстве практических задач их можно получить только численно, а никак не аналитически.
Примечание 3
Примечание 3
Различные линейные и нелинейные варианты рассматриваемого уравнения теплопроводности описывают различные модели физических сред, которые характеризуются определенными зависимостями D(u) и ф(и). В частности, для металлов в большинстве случаев можно считать, что D=const, в то время как для плазмы имеется специфическая зависимость коэффициента диффузии от температуры.
Обратное уравнение теплопроводности
Замечательными свойствами обладает так называемое обратное уравнение диффузии тепла, которое получается путем замены в исходном (прямом) уравнении переменной t на -t. Согласно постановке задачи, обратное уравнение теплопроводности описывает реконструкцию динамики профиля температуры остывающего стержня, если известно начальное условие в виде профиля температуры в некоторый момент времени после начала остывания. Таким образом, требуется определить, как происходило остывание стержня. Мы ограничимся самым простым линейным уравнением с D=const без источников тепла:
Это уравнение гиперболического типа и оно, несмотря на кажущуюся близость к рассмотренным вариантам уравнения теплопроводности, обладает замечательными свойствами.
Если попробовать осуществить расчет обратного уравнения диффузии тепла по тем же самым алгоритмам, что и для обычных уравнений (для этого достаточно в листинге 11.1 или 11.2 заменить значение коэффициента диффузии на отрицательное число, например, D=-1), то мы получим заведомо нефизичное решение. Оно показано на Рисунок 11.5 в виде профилей распределения температуры для нескольких последовательных моментов времени. Как видно, решение выражается в появлении все более быстрых пространственных осцилляции профиля температуры для каждого нового момента времени. Очень существенно, что такое решение является не проявлением неустойчивости численного алгоритма (как для ситуации, рассмотренной в разд. "Устойчивость"этой главы), а определяется спецификой самой задачи.
Оказывается, что обратное уравнение теплопроводности принадлежит к довольно широкому классу задач, называемых некорректными. Некорректные задачи нельзя решать стандартными методами, а для того, чтобы с ними справиться (т. е., чтобы получить осмысленное физическое решение), приходится несколько менять саму их постановку, вводя в нее дополнительную априорную информацию о строении решения.
Решение стационарного двумерного уравнения теплопроводности (см листинг 11 7 ниже)
Рисунок 11.2. Решение стационарного двумерного уравнения теплопроводности (см. листинг 11.7 ниже)
Как несложно сообразить, если искомая функция не зависит от времени, то частная производная по времени в левой части уравнения равна нулю, и само уравнение можно переписать (переобозначив ради упрощения ?<-?/D) следующим образом:
Полученное уравнение, согласно классификации предыдущего раздела, является эллиптическим. Его называют уравнением Пуассона, а для его решения в Matcad предусмотрены две встроенные функции. Если, к тому же, источники равны нулю, то уравнение (11.2), принимающее вид ?u=0, называют уравнением Лапласа.
Одномерное динамическое уравнение
Предположим, что мы рассматриваем задачу распределения тепла не по плоской поверхности, а по удлиненному телу типа металлического стержня (Рисунок 11.3). В этом случае зависимость от координаты у в общем уравнении теплопроводности пропадает, и получается одномерное уравнение:
Одномерное уравнение намного проще двумерного, поскольку объем вычислений для реализации алгоритма его численного решения не так велик. Типичное решение одномерного уравнения диффузии тепла с коэффициентом диффузии о=2, нулевым источником ф=о и начальным распределением температуры в форме нагретой центральной области стержня показано (в виде графика поверхности) на Рисунок 11.4.
Начиная с версии Mathcad 11, для решения одномерных параболических и гиперболических уравнений можно применять встроенную функцию pdesolve.
Решение уравнения диффузии тепла при помощи встроенной функции pdesoдve (листинг 11 4)
Рисунок 11.14. Решение уравнения диффузии тепла при помощи встроенной функции pdesoдve (листинг 11.4)
Пример: волновое уравнение
Приведем еще один пример применения функции pdesoive для решения уравнений в частных производных. Рассмотрим одномерное волновое уравнение, которое описывает, например, свободные колебания струны музыкального инструмента:
Здесь неизвестная функция u(x,t) описывает динамику смещения профиля струны относительно невозмущенного (прямолинейного) положения, а параметр с характеризует материал, из которого изготовлена струна.
Как вы видите, уравнение (11.11) содержит производные второго порядка, как по пространственной координате, так и по времени. Для того чтобы можно было использовать встроенную функцию pdesolve, необходимо переписать волновое уравнение в виде системы двух уравнений в частных производных, введя вторую неизвестную функцию v=ut. Программа для решения волнового уравнения приведена в листинге 11.5, а результат— на Рисунок 11.15.
Решение уравнения Пуассона с помощью функции relax (продолжение листинга 11 8)
Рисунок 11.19. Решение уравнения Пуассона с помощью функции relax (продолжение листинга 11.8)
Разностная схема для решения уравнения Пуассона
Несмотря на отсутствие сведений в справочной системе Mathcad о решении других линейных дифференциальных уравнений в частных производных, кроме уравнения Пуассона, сделать это возможно с помощью той же функции relax (см. предыдущий разд.). Для этого нужно правильным образом задать коэффициенты разностной схемы.
Начнем с пояснения выбора этих коэффициентов (см. листинг 11.8) для уравнения Пуассона. Согласно основным идеям метода сеток (см. разд. 11.2), для дискретизации обеих пространственных производных в уравнении (11.12) следует использовать по три соседних узла вдоль каждой из координат. Поэтому уравнение Пуассона (11.12) может быть записано в разностной форме при помощи шаблона типа "крест" (Рисунок 11.20). В этом случае, после приведения подобных слагаемых в разностных уравнениях коэффициенты разностной схемы будут такими, как показано возле узлов шаблона на этом рисунке (аналогичные коэффициенты для явной и неявных схем решения уравнения теплопроводности см. на Рисунок 11.6 и 11.11 соответственно).
Теперь, если вы сравните полученные числа с константами, которые присвоены элементам матриц-аргументов функции relax (см. листинг 11.8), то увидите, что они как раз и описывают вычисленные нами только что коэффициенты разностной схемы "крест". Таким образом, нетрудно сообразить, что с помощью встроенной функции relax можно решать и другие линейные дифференциальные уравнения в частных производных, которые можно аппроксимировать схемой типа "крест" или схемой, являющейся ее составной частью. Конечно, для того чтобы использовать эту встроенную функцию для другого уравнения, необходимо будет составить соответствующую разностную схему.
Решение уравнения теплопроводности с нелинейным источником (тепловой фронт)
Рисунок 11.8. Решение уравнения теплопроводности с нелинейным источником (тепловой фронт)
Решение уравнения теплопроводности с помощью функции relax (продолжение листинга 11 9)
Рисунок 11.21. Решение уравнения теплопроводности с помощью функции relax (продолжение листинга 11.9)
Решение уравнения теплопроводности
Рисунок 11.9. Решение уравнения теплопроводности с нелинейным источником и коэффициентом диффузии (режим локализации горения)
Читателю предлагается поэкспериментировать с этим и другими нелинейными вариантами уравнения теплопроводности. Существенно, что такие интересные результаты удается получить лишь численно, а в Mathcad только с применением элементов программирования.
Устойчивость
Как мы убедились, явная разностная схема Эйлера дает вполне разумные результаты и вполне может использоваться для практического моделирования задач, связанных с решением уравнений в частных производных. Однако теперь пришло время сказать об очень важной характеристике разностных схем, которая называется их устойчивостью. Не вдаваясь в детали, заметим, что производить расчеты можно только при помощи устойчивых разностных схем, а чтобы пояснить это понятие, обратимся вновь к листингу 11.1, реализующему явную схему для линейного уравнения диффузии.
Слегка изменим соотношение шагов по времени и пространственной координате, произведя расчеты сначала с t=0.0005 (эти результаты показаны на Рисунок 11.7 выше), а также с t=0.0010 и t=0.0015. Результат применения разностной схемы Эйлера для t=0.0010ю примерно тот же, что и для меньшего значения т, приведенного на Рисунок 11.7. А вот следующее (казалось бы, незначительное) увеличение шага по времени приводит к катастрофе (Рисунок 11.10).
Решение волнового уравнения (продолжение листинга 11 5)
Рисунок 11.15. Решение волнового уравнения (продолжение листинга 11.5)
Встроенная функция numol
Альтернативный вариант решения дифференциальных уравнений в частных производных заключается в применении еще одной встроенной функции numo|, которая реализует тот же самый алгоритм сеток, позволяя вручную задать большинство его параметров:
Пример решения волнового уравнения при помощи функции numol приведен в листинге 11.6, особое внимание в котором мы призываем уделить формату представления векторов rhs, init и be, а также принципу извлечения отдельных сеточных решений из матрицы-результата. График решения, показанный на Рисунок 11.16, полезно сравнить с результатом применения вычислительного блока из предыдущего раздела (см. листинг 11.5 и Рисунок 11.15).
Решение волнового уравнения (продолжение листинга 11 6)
Рисунок 11.16. Решение волнового уравнения (продолжение листинга 11.6)
Именно в целях визуализации решения параболических и гиперболических уравнений в частных производных использование функции numol наиболее полезно. График решения динамических уравнений (зависящих от времени t) выглядит намного эффектнее и воспринимается несравненно лучше, если он оформлен в виде анимации. Для создания анимационных роликов расчетное время следует выразить через константу FRAME и затем применить команду View / Animate (Вид / Анимация) (см. разд. 13.3.2).
Шаблон аппроксимации явной схемы для уравнения теплопроводности
Рисунок 11.6. Шаблон аппроксимации явной схемы для уравнения теплопроводности
Множители для каждого из значений сеточной функции в узлах шаблона, соответствующие разностному уравнению (11.9), приведены рядом с каждой точкой шаблона на Рисунок 11.6. Фактически геометрия шаблона и эти множители задают построенную нами разностную схему.
Несложно убедиться в том, что для получения замкнутой системы разностных алгебраических уравнений систему (11.9) необходимо дополнить дискретным представлением начального и граничных условий (11.6) и (11.7). Тогда число неизвестных будет в точности равно числу уравнений, и процесс формирования разностной схемы будет окончательно завершен.
Примечание 1
Примечание 1
Важно подчеркнуть, что возможная нелинейность полученной системы алгебраических уравнений определяется зависимостями от температуры функций D(u) и ф (u), т. е. как коэффициент диффузии, так и источник тепла могут быть функциями сеточной функции ui,k.
Если присмотреться к разностным уравнениям (11.9) повнимательнее, то можно сразу предложить несложный алгоритм реализации этой разностной схемы. Действительно, каждое неизвестное значение сеточной функции со следующего временного слоя, т. е. левая часть соотношения (11.9) явно выражается через три ее значения с предыдущего слоя (правая часть), которые уже известны. Таким образом, в случае уравнения теплопроводности нам очень повезло. Для расчета первого слоя по времени следует попросту подставить в (11.9) начальное условие (известные значения и с нулевого слоя в узлах сетки), для расчета второго слоя достаточно использовать вычисленный таким образом набор и с первого слоя и т. д. Из-за того что разностная схема сводится к такой явной подстановке, ее и называют явной, а благодаря пересчету значений с текущего слоя через ранее вычисленные слои — схемой бегущего счета.
Линейное уравнение
Сделанные замечания относительно реализации явной схемы для уравнения диффузии тепла сразу определяют алгоритм ее программирования в Mathcad. Для решения задачи нужно аккуратно ввести в листинг соответствующие формулы при помощи элементов программирования.
Решение системы разностных уравнений (11.9) для модели без источников тепла, т.е. ф(x,T,t)=0 и постоянного коэффициента диффузии D=const приведено в листинге 11.1. В его первых трех строках заданы шаги по временной и пространственной переменным t и А, а также коэффициент диффузии о, равный единице. В следующих двух строках заданы начальные (нагретый центр области) и граничные (постоянная температура на краях) условия соответственно. Затем приводится возможное программное решение разностной схемы, причем функция пользователя v(t) задает вектор распределения искомой температуры в каждый момент времени (иными словами, на каждом слое), задаваемый целым числом t.
Шаблон аппроксимации уравнения Пуассона "крест"
Рисунок 11.20. Шаблон аппроксимации уравнения Пуассона "крест"
Решение уравнения диффузии тепла при помощи функции relax
Приведем пример применения встроенной функции relax для решения другого уравнения в частных производных (т. е. не уравнения Пуассона, для которого она изначально предназначена). Вычислим при помощи этой функции решение уже хорошо нам знакомого однородного линейного уравнения теплопроводности (см. разд. 11.1.2). Будем использовать явную разностную схему, шаблон которой изображен на Рисунок 11.6. Для того чтобы "приспособить" для явной схемы функцию relax, требуется только задать ее аргументы в соответствии с коэффициентами, показанными на шаблоне (см. Рисунок 11.6). Программа, реализующая таким способом явную схему, представлена на листинге 11.9. Число Куранта в этом листинге обозначено переменной с, как и положено явной разностной схеме, она выдает устойчивое решение только для C<1.
Шаблон неявной схемы для уравнения теплопроводности
Рисунок 11.11. Шаблон неявной схемы для уравнения теплопроводности
Очень важно, что если само уравнение теплопроводности линейно, то с в левой части разностного уравнения является константой, а ф в его правой части может зависеть только от первой степени и. Поэтому система уравнений (11.10) для всех пространственных узлов 1=1.. .м-l является линейной системой, что существенно упрощает ее решение (поскольку известно, что для линейных систем с ненулевым определителем решение существует и является единственным). Напомним, что для получения замкнутой системы линейных уравнений необходимо дополнить данный набор разностных уравнений граничными условиями, т. е. известными значениями сеточной функции для i=0 и i=M.
Примечание 1
Примечание 1
Если рассматривать нелинейный случай, то на каждом шаге по времени пришлось бы решать систему нелинейных уравнений, число решений которых могло бы быть большим, и среди них требовалось бы отыскать нужное, а не паразитное решение.
Для реализации неявной схемы, таким образом, можно использовать комбинацию средств программирования Mathcad и встроенной функции решения системы линейных уравнений isolve. Один из возможных способов решения предложен в листинге 11.2. Большая часть этого листинга является вводом параметров задачи (шагов, начальных и граничных условий), и только в последней его строке определяется функция пользователя, вычисляющая сеточную функцию на каждом временном слое (при помощи встроенной функции решения системы линейных уравнений isolve). В нескольких предыдущих строках листинга (после расчета коэффициента Куранта сои) формируется матрица системы уравнений А, которая записывается в подходящем для Mathcad виде, как это сделано в листинге 11.2. Как несложно убедиться, столбец правых частей разностных уравнений выражается вычисленными значениями сеточной функции с предыдущего слоя.
Результаты расчетов по неявной схеме показаны на Рисунок 11.12 и, как видно, они дают примерно те же результаты, что и в случае применения явной схемы (см. Рисунок 11.7). Обратите внимание, что решение устойчиво при любых значениях коэффициента Куранта (в том числе и больших 1), поскольку, как следует из соответствующих положений теории численных методов, неявная схема является безусловно-устойчивой.
Встроенные функции для решения уравнений в частных производных
11.3. Встроенные функции для решения уравнений в частных производныхКак видно из предыдущего раздела, с уравнениями в частных производных вполне можно справиться и не прибегая к специфическим средствам Mathcad. Между тем в версиях Mathcad 11 (и выше) имеется несколько встроенных функций, при помощи которых можно автоматизировать процесс решения дифференциальных уравнений в частных производных. Рассмотрим в данном разделе основные аспекты их применения, отмечая не только инструкции по их применению, описанные разработчиками Mathcad, но также и некоторые "скрытые" возможности этих функций.
Матричные вычисления в Mathcad
Генерация псевдослучайных чисел
12.3.1. Генерация псевдослучайных чиселКак уже отмечалось (см. разд. 12.1.1), для генерации M-компонентного вектора независимых псевдослучайных чисел имеется ряд встроенных функций, реализующих различные типы статистических распределений и имеющих вид r*(M,раr), где * — идентификатор, a par— список параметров конкретного распределения. В частности, генератор нормальных псевдослучайных чисел был рассмотрен ранее (см. разд. 12.1.2), а для равномерного распределения предусмотрено две встроенных функции:
Чаще всего в несложных программах применяется последняя функция, которая приводит к генерации одного псевдослучайного числа. Наличие такой встроенной функции в Mathcad — дань традиции, применяемой в большинстве сред программирования. Возможно, именно скалярный тип этой функции определяет простоту и привычность ее использования в расчетах моделей типа Монте-Карло.
Гистограммы
12.2.1. ГистограммыГистограммой называется график, аппроксимирующий по случайным данным плотность их распределения. При построении гистограммы область значений случайной величины (а,b) разбивается на некоторое количество bin сегментов, а затем подсчитывается процент попадания данных в каждый сегмент. Для построения гистограмм в Mathcad имеется несколько встроенных функций. Рассмотрим их, начиная с самой сложной по применению, чтобы лучше разобраться в возможностях каждой из функций.
Гистограммы с произвольными интервалами
Если вектор intvls имеет bin элементов, то и результат hist имеет столько же элементов. Построение гистограммы иллюстрируется листингом 12.8 и Рисунок 12.6.
Статистические функции
12.1.1. Статистические функцииВ Mathcad имеется ряд встроенных функций, задающих используемые в математической статистике законы распределения. Они вычисляют как значение плотности вероятности различных распределений по значению случайной величины х, так и некоторые сопутствующие функции. Все они, по сути, являются либо встроенными аналитическими зависимостями, либо специальными функциями. Большой интерес представляет наличие генераторов случайных чисел, создающих выборку псевдослучайных данных с соответствующим законом распределения, что является основой методов Монте-Карло (см. разд. 12.2).
В Mathcad заложена информация о большом количестве разнообразных статистических распределений, включающая, с одной стороны, табулированные функции вероятности, и, с другой, возможность генерации последовательности случайных чисел с соответствующим законом распределения. Для реализации этих возможностей имеются четыре основных категории встроенных функций. Их названия являются составными и устроены одинаковым образом: первая литера идентифицирует определенный закон распределения, а оставшаяся часть (ниже в списке функций она условно обозначена звездочкой) задает смысловую часть встроенной функции:
Чтобы получить функции, относящиеся, например, к равномерному распределению, вместо * надо поставить unif и ввести соответствующий список параметров par. Он будет состоять в данном случае из двух чисел: а,b — границ интервала распределения случайной величины.
Перечислим все типы распределения, реализованные в Mathcad, вместе с их параметрами, на этот раз обозначив звездочкой * недостающую первую букву встроенных функций. Некоторые из плотностей вероятности показаны на Рисунок 12.1.
Генерация коррелированных выборок
12.3.2. Генерация коррелированных выборокДо сих пор мы рассматривали наиболее простой случай применения генераторов независимых случайных чисел. В методах Монте-Карло часто требуется создавать случайные числа с определенной корреляцией. Приведем пример программы, создающей два вектора x1 и х2 одинакового размера и одним и тем же распределением, случайные элементы которых попарно коррелированны с коэффициентом корреляции R (листинг 12.21).
Пример нормальное (Гауссово) распределение
12.1.2. Пример: нормальное (Гауссово) распределениеВ теории вероятности доказано, что сумма различных независимых случайных слагаемых (независимо от их закона распределения) оказывается случайной величиной, распределенной согласно нормальному закону (так называемая центральная предельная теорема). Поэтому нормальное распределение хорошо моделирует самый широкий круг явлений, для которых известно, что на них влияют несколько независимых случайных факторов.
Перечислим еще раз встроенные функции, имеющиеся в Mathcad для описания нормального распределения вероятностей:
Математическое ожидание и дисперсия являются, по сути, параметрами распределения. Плотность распределения для трех пар значений параметров показана на Рисунок 12.3. Напомним, что плотность распределения dnorm задает вероятность попадания случайной величины х в малый интервал от х до х+?х. Таким образом, например, для первого графика (сплошная линия) вероятность того, что случайная величина х примет значение в окрестности нуля, приблизительно в три раза больше, чем вероятность того, что она примет значение в окрестности х=2. А значения случайной величины большие 5 и меньшие -5 и вовсе очень маловероятны.
Среднее и дисперсия
12.2.2. Среднее и дисперсияВ Mathcad имеется ряд встроенных функций для расчетов числовых статистических характеристик рядов случайных данных:
Пример использования первых четырех функций приведен в листинге 12.10.
Моделирование случайного процесса
12.3.3. Моделирование случайного процессаВстроенные функции для генерации случайных чисел создают выборку из случайных данных Ai. Часто требуется создать непрерывную или дискретную случайную функцию A(t) одной или нескольких переменных (случайный процесс или случайное поле), значения которой будут упорядочены относительно своих переменных. Создать псевдослучайный процесс можно способом, представленным в листинге 12.22.
Примеры Выборочная оценка дисперсии и среднего нормальной случайной величины
12.2.3. Примеры: Выборочная оценка дисперсии и среднего нормальной случайной величиныТиповые задачи математической статистики связаны с получением тех или иных интервальных и точечных оценок различных параметров случайной выборки. Приведем пример двух задач, иллюстрирующих назначение и принципы применения введенных в предыдущих разделах статистических функций.
Интервальная оценка дисперсии
Требуется определить числовой интервал (L,U), внутри которого будет лежать с вероятностью 1-а=75% дисперсия нормальной случайной величины, исходя из объема выборки в N чисел. Эта задача решается в статистике с помощью ?2-распределения (листинг 12.13).
Корреляция
12.2.4. КорреляцияФункции, устанавливающие связь между парами двух случайных векторов, называются ковариацией и корреляцией (или, по-другому, коэффициентом корреляции). Они различаются нормировкой, как следует из их определения (листинг 12.16):
Пример огибающая и фаза нормального случайного процесса
12.3.4. Пример: огибающая и фаза нормального случайного процессаЗавершим разговор о моделировании случайных процессов примером, часто встречающимся в задачах статистической радиофизики. Рассмотрим модель, представляющую собой сумму гармонической функции и нормально распределенной шумовой компоненты, которая хорошо описывает передачу сигнала в электронных устройствах в условиях помех (листинги 12.24—12.25) и называется узкополосным нормальным процессом. Как известно, узкополосный процесс представим в виде E(t)exp(iф(t)), где случайные функции E(t) и ф(t) называются, соответственно, его огибающей и фазой. Мы приведем пример нулевого сигнала (т. е. расчет огибающей и фазы чистого случайного процесса), хотя минимальное изменение листинга 12.25 даст вам возможность промоделировать и ненулевые значения сигнала.
Первая половина листинга 12.24 представляет собой подготовительный этап, заключающийся в генерации двух векторов с нормальным распределением вероятности. Из курса математической статистики известно, что узкополосный гауссов процесс можно представить в виде, приведенном в последней строке листинга 12.24, причем случайные функции A(t) и C(t) называются квадратурными составляющими нормального случайного процесса. Графики A(t) и C(t) показаны на Рисунок 12.13.
Новые функции корреляционного анализа сигналов
12.2.5. Новые функции корреляционного анализа сигналовВ Mathcad 12 появились две функции, связанные с корреляционной обработкой сигналов и изображений:
Смысл действия одномерной функции correl заключается в последовательном сдвиге одного вектора относительно другого, перемножении и суммировании их элементов, стоящих в таком положении друг напротив друга. Это не совсем ковариация в терминах математической статистики, поскольку не осуществляется нормировки (деления полученной суммы на число перемноженных элементов). Работа двумерной функции correl2d состоит в последовательном позиционировании "маленькой" матрицы (окна) на фоне "крупной" (прототипа), перемножении их элементов, находящихся друг над другом, и суммировании их. В результате получается элемент матрицы корт реляции, соответствующий центрам наложения матрицы-прототипа и матрицы-окна. Обе эти функции играют определенную роль в задачах обработки сигналов.
Коэффициенты асимметрии и эксцесса
12.2.6. Коэффициенты асимметрии и эксцессаКоэффициент асимметрии задает степень асимметричности плотности вероятности относительно оси, проходящей через ее центр тяжести. Коэффициент асимметрии определяется третьим центральным моментом распределения. В любом симметричном распределении с нулевым математическим ожиданием, например, нормальным, все нечетные моменты, в том числе и третий, равны нулю, поэтому коэффициент асимметрии тоже равен нулю.
Степень сглаженности плотности вероятности в окрестности главного максимума задается еще одной величиной — коэффициентом эксцесса. Он показывает, насколько острую вершину имеет плотность вероятности по сравнению с нормальным распределением. Если коэффициент эксцесса больше нуля, то распределение имеет более острую вершину, чем распределение Гаусса, если меньше нуля, то более плоскую.
Для расчета коэффициентов асимметрии и эксцесса в Mathcad имеются две встроенные функции:
Примеры расчета коэффициентов асимметрии и эксцесса для распределения Вейбулла приведены в листинге 12.17.
Статистические функции матричного аргумента
12.2.7. Статистические функции матричного аргументаВсе рассмотренные примеры работы статистических функций относились к векторам, элементы которых были случайными числами. Но точно так же все эти функции применяются и по отношению к выборкам случайных данных, сгруппированных в матрицы. При этом статистические характеристики рассчитываются для совокупности всех элементов матрицы, без разделения ее на строки и столбцы. Например, если матрица имеет размерность MxN, то и объем выборки будет равен M-N.
Соответствующий пример вычисления среднего значения приведен в листинге 12.18. В его первой строке определяется матрица данных х размера 4x2. Действие встроенной функции mean матричного аргумента (последняя строка листинга) иллюстрируется явным суммированием элементов матрицы х (предпоследняя строка). Действие прочих встроенных функций на матрицы совершенно аналогично действию их на векторы (листинг 12.19).
Демонстрирует широко распространенный
Рисунок 12.9 демонстрирует широко распространенный прием моделирования, основанный на смеси некоторого полезного сигнала f(х) и шумовой компоненты а, что характерно для типичного физического эксперимента. В качестве шума используется серия равномерно распределенных псевдослучайных чисел. Различные соотношения интенсивностей сигнала и шума определяют различные условия модельной задачи и позволяют эффективно протестировать алгоритмы, которые разработаны для ее решения.Примечание 1
Более сложные статистические расчеты для модели сигнал/шум приведены в листинге 12.24 (см. разд. 12.3.4).
Параметры генераторов псевдослучайных чисел
В Mathcad применяются типичные алгоритмы генерации последовательностей псевдослучайных чисел, которые используют в качестве "отправной точки" некоторое начальное (или отправное) значение (seed value). Это начальное значение используется для того, чтобы совершить над ним определенные математические действия (к примеру, взять остаток от деления на некоторое другое число) и получить в итоге первое псевдослучайное число последовательности. Затем те же математические операции совершаются с первым числом для получения второго, и т. д.
Несложно догадаться, что если использовать все время одно и то же начальное значение генератора псевдослучайных чисел, то, открывая всякий раз новый документ со встроенной функцией получения тех или иных псевдослучайных чисел, будет выдаваться в точности одна и та же их последовательность. Сами числа внутри последовательности будут "почти" случайными (значимость этого "почти" будет зависеть только от качества алгоритма генерации), но вот сама последовательность при каждом открытии документа будет одной и той же.
Примечание 2
Примечание 2
Вы можете сами убедиться в идентичности выдаваемых последовательностей псевдослучайных чисел, открывая и закрывая несколько раз документ с Рисунок 12.9 или любой документ с иным датчиком случайных чисел, описанных в данном разделе. Любопытным также будет эксперимент по последовательному применению нескольких встроенных функций генерации псевдослучайных чисел в пределах одного документа.
Для того чтобы иметь возможность поменять саму последовательность сгенерированных случайным датчиком чисел, в Mathcad 11 (и выше) предусмотрена новая возможность явного определения используемого им начального значения. Для этого вызовите командой меню Tools/ Worksheet Options (Сервис / Опции документа) диалог Worksheet Options (Опции документа) и на вкладке Built-in Variables (Встроенные переменные) введите желаемое (целое) число в поле ввода Seed value for random Bombers (Начальное значение для генератора псевдослучайных чисел).
Альтернативный способ заключается в использовании соответствующей встроенной функции непосредственно в документе:
Диалоговое окно Insert Function
Рисунок 12.2. Диалоговое окно Insert Function
Статистика
Глава 12. Статистика| Содержание |
График и матрица гистограммы (продолжение листинга 12 9)
Рисунок 12.7. График и матрица гистограммы (продолжение листинга 12.9)
Примечание 2
Примечание 2
Для того чтобы назначить двумерному графику тип гистограммы, в диалоговом окне Formatting Currently Selected Graph (Форматирование) установите на вкладке Traces (Графики) тип списка bar (Столбцы) или solidbar (Гистограмма). На Рисунок 12.6 и 12.7 применены установки второго типа: закрашенными столбиками (solidbar).
К задаче проверки статистических гипотез (продолжение листинга 12 14)
Рисунок 12.8. К задаче проверки статистических гипотез (продолжение листинга 12.14)
В листинге 12.15 приводится альтернативный способ проверки той же самой гипотезы, связанный с вычислением значения не квантиля, а самого распределения Стьюдента.
Корреляционная функция (продолжение листингов 12 22—12 23)
Рисунок 12.12. Корреляционная функция (продолжение листингов 12.22—12.23)
Примечание 2
Примечание 2
Внимательному читателю предлагается самостоятельно ответить на вопрос: почему при таком расчете корреляционной функции её значение R (0) не равно 1, как должно быть по определению?
Квадратурные составляющие случайного процесса (продолжение листинга 12 24)
Рисунок 12.13. Квадратурные составляющие случайного процесса (продолжение листинга 12.24)
Вероятность того, что
Листинг 12.1. Вероятность того, что х будет меньше 1.881
Расчет числовых характеристик случайного вектора
Листинг 12.10. Расчет числовых характеристик случайного вектора
Определение статистических характеристик случайных величин приведено в листинге 12.11 на еще одном примере обработки выборки малого объема (по пяти данным). В том же листинге иллюстрируется применение еще двух функций, которые имеют смысл дисперсии и стандартного отклонения в несколько другой нормировке. Сравнивая различные выражения, вы без труда освоите связь между встроенными функциями.
ВНИМАНИЕ!
Осторожно относитесь к написанию первой литеры в этих функциях, особенно при обработке малых выборок (листинг 12.11).
К определению статических характеристик
Листинг 12.11. К определению статических характеристик
Иногда в статистике встречаются и иные функции, например, помимо арифметического среднего, применяются другие средние значения:
Математическое определение этих функций и пример их использования в Mathcad приведены в листинге 12.12.
Вычисление различных средних значений
Листинг 12.12. Вычисление различных средних значений
Интервальное оценивание дисперсии
Листинг 12.13. Интервальное оценивание дисперсии
Указанный интервал называется (1-a) доверительным интервалом. Обратите внимание на использование при решении данной задачи функции stdev (с прописной буквы) для расчета выборочного стандартного отклонения. В статистике часто встречаются выражения, которые более удобно записывать через функции в такой нормировке, именно для этого они и появились в Mathcad.
Проверка статистических гипотез
В статистике рассматривается огромное число задач, связанных с проверкой тех или иных гипотез н. Разберем пример простой гипотезы. Пусть имеется выборка N чисел с нормальным законом распределения и неизвестными дисперсией и математическим ожиданием. Требуется принять или отвергнуть гипотезу H о том, что математическое ожидание закона распределения равно некоторому числу µ0=0. 2.
Задачи проверки гипотез требуют задания уровня критерия проверки гипотезы а, который описывает вероятность ошибочного отклонения истинной н. Если взять ее очень малым, то гипотеза, даже если она ложная, будет почти всегда приниматься; если, напротив, взять а близкими 1, то критерий будет очень строгим, и гипотеза, даже верная, скорее всего, будет отклонена.
В нашем случае гипотеза состоит в том, что µ0=0.2, а альтернатива— что µ0#0.2. Оценка математического ожидания, как следует из курса классической статистики, решается с помощью распределения Стьюдента с параметром N-1 (этот параметр называется степенью свободы распределения).
Для проверки гипотезы (листинг 12.14) рассчитывается (а/2)-квантиль распределения Стьюдента т, который служит критическим значением для принятия или отклонения гипотезы. Если соответствующее выборочное значение T по модулю меньше т, то гипотеза принимается (считается верной). В противном случае гипотезу следует отвергнуть.
Квадратурные составляющие нормального случайного процесса
Листинг 12.14. Квадратурные составляющие нормального случайного процесса
Проверка гипотезы о математическом ожидании при неизвестной дисперсии
Листинг 12.14. Проверка гипотезы о математическом ожидании при неизвестной дисперсии
В последней строке листинга вычисляется истинность или ложность условия, выражающего решение задачи. Поскольку условие оказалось ложным (равным не 1, а 0), то гипотезу необходимо отвергнуть.
На Рисунок 12.8 показано распределение Стьюдента с N-l степенью свободы, вместе с критическими значениями, определяющими соответствующий интервал. Если t (оно тоже показано на графике) попадает в него, то гипотеза принимается; если не попадает (как произошло в данном случае) — отвергается. Если увеличить а, ужесточив критерий, то границы интервала будут сужаться по сравнению с показанными на рисунке.
Другой вариант проверки гипотезы (продолжение листинга 12 14)
Листинг 12.15. Другой вариант проверки гипотезы (продолжение листинга 12.14)
Мы разобрали только два характерных примера статистических вычислений. Однако с помощью Mathcad легко решаются самые разнообразные задачи математической статистики.
Примечание 1
Примечание 1
Большое количество задач разобрано в Ресурсах в рубрике Statistics (Статистика) Быстрых шпаргалок.
Расчет ковариации и корреляции
Листинг 12.16. Расчет ковариации и корреляции
Примечание 1
Примечание 1
Как видно, полученное значение корреляции близко к нулю, т. е. псевдослучайные числа генерируются практически независимо. О том, как создать массив коррелированных чисел, написано ниже (см. разд. 12.3.2).
Расчет выборочных коэффициентов асимметрии и эксцесса
Листинг 12.17. Расчет выборочных коэффициентов асимметрии и эксцесса
Вычисление среднего значения элементов матрицы
Листинг 12.18. Вычисление среднего значения элементов матрицы
Действие различных статистических функций на матрицу
Листинг 12.19. Действие различных статистических функций на матрицу="31.gif" >
Примечание 1
Примечание 1
Некоторые статистические функции (например, вычисления ковариации) имеют два аргумента. Они также могут быть матрицами, но, в соответствии со смыслом функции, должны иметь одинаковую размерность.
Большинству статистических функций позволяется иметь в качестве аргументов даже не одну матрицу, а любое количество матриц, векторов и скаляров. Числовые характеристики будут рассчитаны для всей совокупности значений аргументов функции. Соответствующий пример приведен в листинге 12.20.
квантиль нормального распределения
Листинг 12.2. 97-% квантиль нормального распределения
Статистические функции нескольких аргументов
Листинг 12.20. Статистические функции нескольких аргументов
Генерация попарно коррелированных случайных чисел
Листинг 12.21. Генерация попарно коррелированных случайных чисел
Результат действия программы для R=0.4 показан на Рисунок 12.10 (слева). Сравните полученную выборку с правым графиком, полученным для высокой корреляции (R=0.9) и с Рисунок 12.5 (см. разд. 12.1.2) для независимых данных, т. е. R=0.
Генерация псевдослучайного процесса
Листинг 12.22. Генерация псевдослучайного процесса
В первых строках листинга 12.22 определено количество N независимых случайных чисел, которые будут впоследствии сгенерированы, и радиус временной корреляции т. В следующих трех строках определяются моменты времени т,, которым будут отвечать случайные значения A(tj). Создание нормального случайного процесса сводится к генерации обычным способом вектора независимых случайных чисел х и построению интерполяционной зависимости в промежутках между ними. В листинге 12.22 используется сплайн-интерполяция (см. главу 13).
В результате получается случайный процесс A(t), радиус корреляции которого определяется расстоянием т между точками, для которых строится интерполяция. График случайного процесса A(t) вместе с исходными случайными числами показан на Рисунок 12.11. Случайное поле можно создать несколько более сложным способом с помощью многомерной интерполяции.
Примечание 1
Примечание 1
Простой пример генерации двумерного случайного поля вы найдете на компакт-диске, прилагаемом к книге.
Дискретизация случайного
Листинг 12.23. Дискретизация случайного процесса и вычисление корреляционной функции (продолжение листинга 12.19)
Дискретизация интервала (0,Tmах) для случайного процесса A(t) произведена с различным элементарным интервалом А (первая строка листинга). В зависимости от значения А получается различный объем n выборки случайных чисел Yi, являющихся значениями случайной функции A(t) в точках дискретизации. В последних четырех строках определяются различные характеристики случайной величины Y, являющиеся, по сути, характеристиками случайного процесса A(t). График рассчитанной в 2-М+1 точках корреляционной функции R(j) показан на Рисунок 12.12.
Огибающая и фаза нормального случайного процесса (продолжение листинга 12 24)
Листинг 12.25. Огибающая и фаза нормального случайного процесса (продолжение листинга 12.24)
содержит суммирование
Листинг 12.25 содержит суммирование полученных шумовых компонент с составляющими гармонического сигнала и выдает в качестве результата функции E(t) и ф(t) (они показаны на Рисунок 12.14).Вероятность того, что
Листинг 12.3. Вероятность того, что х будет больше 2
Вероятность того что х будет находиться в интервале (2 3)
Листинг 12.4. Вероятность того, что х будет находиться в интервале (2,3)
Обратите внимание, что задачи двух
Листинг 12.5. Вероятность того, что |x|<2
Обратите внимание, что задачи двух последних листингов решаются двумя разными способами. Второй из них связан с еще одной встроенной функцией erf, называемой функцией ошибок (или интегралом вероятности, или функцией Крампа).
Математический смысл функции ошибок ясен из листинга 12.5. Интеграл вероятности имеет всего один аргумент, в отличие от функции нормального распределения. Исторически последняя пересчитывалась через табулированный интеграл вероятности по формулам, приведенным в листинге 12.6 для произвольных значений параметров µ и ?(листинг 12.6).
Вероятность того что х будет в интервале (2 3)
Листинг 12.6. Вероятность того, что х будет в интервале (2,3)
Если вы имеете дело с моделированием методами Монте-Карло, то в качестве генератора случайных чисел с нормальным законом распределения применяйте встроенную функцию rnorm. В листинге 12.7 ее действие показано на примере создания двух векторов по M=500 элементов в каждом с независимыми псевдослучайными числами x1i и х2i распределенными согласно нормальному закону. О характере распределения случайных элементов векторов можно судить по Рисунок 12.5. В дальнейшем мы будем часто сталкиваться с генерацией случайных чисел и расчетом их различных средних характеристик.
Генерация двух векторов с нормальным законом распределения
Листинг 12.7. Генерация двух векторов с нормальным законом распределения
Построение гистограммы
Листинг 12.8. Построение гистограммы
Для анализа взято N=1000 данных с нормальным законом распределения, созданных генератором случайных чисел (третья строка листинга). Далее определяются границы интервала (upper, lower), содержащего внутри себя все случайные значения, и осуществляется его разбиение на количество (bin) одинаковых сегментов, начальные точки которых записываются в вектор int (предпоследняя строка листинга).
Примечание 1
Примечание 1
В векторе int можно задать произвольные границы сегментов разбиения так, чтобы они имели разную ширину.
Упрощенный вариант построения гистограммы
Листинг 12.9. Упрощенный вариант построения гистограммы
Методы МонтеКарло
12.3. Методы Монте-КарлоДля моделирования различных физических, экономических и прочих эффектов широко распространены методы, называемые методами Монте-Карло, обязанные своим названием европейскому центру азартных игр, основанных на случайных событиях. Основная идея этих методов состоит в создании определенной последовательности случайных чисел, моделирующей тот или иной эффект, например, шум в физическом эксперименте, случайную динамику биржевых индексов и т. п. Фактически все содержание предыдущего раздела (см. разд. 12.2) являлось примером реализации методов Монте-Карло, поскольку содержало расчет тех или иных статистических характеристик имитационных данных, полученных при помощи генератора случайных чисел.
Рассмотрим еще несколько типовых задач, относящихся к методам Монте-Карло, попутно изучая соответствующие возможности, заложенные в Mathcad.
Модель сигнал/шум с равномерным законом распределения
Рисунок 12.9. Модель сигнал/шум с равномерным законом распределения
Нормальные функции распределения
Рисунок 12.4. Нормальные функции распределения
Функция распределения F(x) (cumulative probability) — это вероятность того, что случайная величина примет значение, меньшее или равное х. Как следует из математического смысла, она является интегралом от плотности вероятности в пределах от -? до х. Функции распределения для упомянутых нормальных законов изображены на Рисунок 12.4. Функция, обратная F(x) (inverse cumulative probability), называемая еще квантилем распределения, позволяет по заданному аргументу р определить значение х, причем случайная величина будет меньше или равна х с вероятностью р.
Примечание 1
Примечание 1
Здесь и далее графики различных статистических функций, показанные на рисунках, получены с помощью Mathcad без каких-либо дополнительных выражений в рабочей области.
Приведем несколько примеров, позволяющих почувствовать математический смысл рассмотренных функций на примере случайной величины х, распределенной по нормальному закону с µ=0 и ?=1 (листинги 12.1—12.5).
Огибающая и фаза нормального случайного процесса (продолжение листингов 12 24—25)
Рисунок 12.14. Огибающая и фаза нормального случайного процесса (продолжение листингов 12.24—25)
Плотность вероятности некоторых распределений
Рисунок 12.1. Плотность вероятности некоторых распределений
Примечание 1
Примечание 1
Математический смысл каждой из четырех типов функций будет объяснен в следующем разделе на примере распределения Гаусса.
Вставку рассмотренных статистических функций в программы удобно осуществлять с помощью диалогового окна Insert Function (Вставка функции). Для этого необходимо выполнить следующие действия:
1. Установите курсор на место вставки функции в документе.
2. Вызовите диалоговое окно Insert Function нажатием кнопки f(x) на стандартной панели инструментов, или командой меню Insert / Function (Вставка / Функция), или нажатием клавиш
3. В списке Function Category (Категория функции) (Рисунок 12.2) выберите одну из категорий статистических функций. Категория Probability Density (Плотность вероятности) содержит встроенные функции для плотности вероятности, Probability Distribution (Функция распределения) — для вставки функций или квантилей распределения, Random Numbers (Случайные числа) — для вставки функции генерации случайных чисел.
4. В списке Function Name (Имя функции) выберите функцию в зависимости от требующегося закона распределения. При выборе того или иного элемента списка в текстовых полях в нижней части окна будет появляться информация о назначении выбранной функции.
5. Нажмите кнопку ОК для вставки функции в документ.
Плотность вероятности нормальных распределений
Рисунок 12.3. Плотность вероятности нормальных распределений
Построение гистограммы (продолжение листинга 12 8)
Рисунок 12.6. Построение гистограммы (продолжение листинга 12.8)
Обратите внимание, что в последней строке листинга осуществлена нормировка значений гистограммы, с тем чтобы она правильно аппроксимировала плотность вероятности, также показанную на графике.
Гистограммы с равными интервалами
Если нет необходимости задавать сегменты гистограммы разной ширины, то удобнее воспользоваться упрощенным вариантом функции hist:
Для того чтобы использовать этот вариант функции hist вместо предыдущего, достаточно заменить первый из ее аргументов в листинге 12.8 следующим образом:
Недостаток упрощенной формы функции hist в том, что по-прежнему необходимо дополнительно определять вектор сегментов построения гистограммы.
От этого недостатка свободна функция histogram:
Примеры использования функции histogram приведены в листинге 12.9 и на Рисунок 12.7. Сравнение с предыдущим листингом подчеркивает простоту построения гистограммы этим способом (стоит отметить, что в листинге 12.9, в отличие от предыдущего, мы не нормировали гистограмму).
Псевдослучайные числа
Рисунок 12.10. Псевдослучайные числа с корреляцией R=0.4 (продолжение листинга"12.21) и R=0.9
Псевдослучайные числа с нормальным законом распределения (продолжение листинга 12 7)
Рисунок 12.5. Псевдослучайные числа с нормальным законом распределения (продолжение листинга 12.7)
Псевдослучайный процесс (продолжение листинга 12 22)
Рисунок 12.11. Псевдослучайный процесс (продолжение листинга 12.22)
К случайным процессам, сгенерированным таким способом, как и к данным эксперимента, применяются любые статистические методы обработки, например, корреляционный или спектральный анализ. Приведем в качестве примера листинг 12.23, показывающий, как организовать расчет корреляционной функции случайного процесса.
Статистические распределения
12.1. Статистические распределенияФундаментальным понятием математической статистики является понятие случайного числа (случайной величины). Согласно определению, случайная величина принимает то или иное значение, но какое конкретно — зависит от принципиально непредсказуемых обстоятельств опыта и заранее точно предсказано быть не может. Можно лишь говорить о вероятности Р(хк) принятия случайной дискретной величиной того или иного значения хk или о вероятности попадания непрерывной случайной величины в тот или иной числовой интервал (х,х+?х). Вероятность Р(хk) или P(х)(?х) соответственно может принимать значения от о (такое значение случайной величины совершенно невероятно) до i (случайная величина заведомо примет значение от х до (х+?х). Соотношение Р(хk) называют законом распределения случайной величины, а зависимость P(х) между возможными значениями непрерывной случайной величины и вероятностями попадания в их окрестность называется ее плотностью вероятности (probability density).
Статистика
СтатистикаMathcad имеет развитый аппарат работы с задачами математической статистики и обработки эксперимента.
Во-первых, имеется большое количество встроенных специальных функций, позволяющих рассчитывать плотности вероятности и другие основные характеристики основных законов распределения случайных величин (см. разд. 12.1). Наряду с этим, в Mathcad запрограммировано соответствующее количество генераторов псевдослучайных чисел для каждого закона распределения (см. разд. 12. Г), что позволяет эффективно проводить моделирование методами Монте-Карло.
Во-вторых, предусмотрена возможность построения гистограмм и расчета статистических характеристик выборок случайных чисел и случайных процессов, таких как средние, дисперсии, корреляции и т. п. (см. разд. 12.2). При этом случайные последовательности могут как создаваться генераторами случайных чисел (методы Монте-Карло, см. разд. 12.3), так и вводиться пользователем из файлов.
В-третьих, имеется целый арсенал средств, направленных на интерполяцию-экстраполяцию данных, построение регрессии по методу наименьших квадратов, фильтрацию сигналов. Наконец, реализован ряд численных алгоритмов, осуществляющих расчет различных интегральных преобразований, что позволяет организовать спектральный анализ различного типа.
Выборочные статистические характеристики
12.2. Выборочные статистические характеристикиВ большинстве статистических расчетов вы имеете дело с выборками: либо со случайными данными, полученными в ходе какого-либо эксперимента (которые выводятся из файла или печатаются непосредственно в документе), либо с результатами генерации случайных чисел, рассмотренными в предыдущих разделах встроенными функциями, моделирующими то или иное явление методом Монте-Карло (см. разд. 12.3). Рассмотрим возможности Mathcad по оценке функций распределения и расчету числовых характеристик случайных данных.
Матричные вычисления в Mathcad
Линейная интерполяция
13.1.1. Линейная интерполяцияСамый простой вид интерполяции — линейная, которая представляет искомую зависимость А(Х) в виде ломаной линии. Интерполирующая функция А(х) состоит из отрезков прямых, соединяющих точки (Рисунок 13.2).
Линейная регрессия
13.2.1. Линейная регрессияСамый простой и наиболее часто используемый вид регрессии — линейная. Приближение данных (xi,yi) осуществляется линейной функцией у(х) = =b+ах. На координатной плоскости (х,у) линейная функция, как известно, представляется прямой линией (Рисунок 13.14). Еще линейную регрессию часто называют методом наименьших квадратов, поскольку коэффициенты а и b вычисляются из условия минимизации суммы квадратов ошибок |b+axi-yi|.
Примечание 1
Примечание 1
Чаще всего такое же условие ставится и в других задачах регрессии, т. е. приближения массива данных (xi,yi) другими зависимостями у(х). Исключение рассмотрено в листинге 13.9.
Примечание 2
Примечание 2
Различным расчетным аспектам реализации метода наименьших квадратов, в большинстве случаев сводящимся к решению систем алгебраических линейных уравнений, была посвящена значительная часть главы 8.
Для расчета линейной регрессии в Mathcad имеются два дублирующих друг друга способа. Правила их применения представлены в листингах 13.7 и 13.8. Результат обоих листингов получается одинаковым (Рисунок 13.14):
Ввод/вывод в текстовые файлы
13.3.1. Ввод/вывод в текстовые файлыПеречислим встроенные функции для работы с текстовыми файлами, которые имеются в Mathcad 2001—12.
Примечание 1
Примечание 1
В Mathcad 12 имеется дополнительная универсальная встроенная функция READFILE, значительно облегчающая процесс импорта данных.
Примечание 2
Примечание 2
Можно задавать как полный путь к файлу, например, С:\Мои документы, так и относительный, имея в виду, что он будет отсчитываться от папки, в которой находится файл с документом Mathcad. Если вы задаете в качестве аргумента просто имя файла (как в листингах 13.16—13.17), то файл будет записан или прочитан из той папки, в которой находится сам документ Mathcad.
Примеры использования встроенных функций иллюстрируются, листингами 13.16—13.18. Результат действия листингов 13.16 и 13.18 можно понять, просмотрев получающиеся текстовые файлы, например, с помощью Блокнота Windows (Рисунок 13.20 и 13.21 соответственно).
Кубическая сплайнинтерполяция
13.1.2. Кубическая сплайн-интерполяцияВ большинстве практических приложений желательно соединить экспериментальные точки не ломаной линией, а гладкой кривой. Лучше всего для этих целей подходит интерполяция кубическими сплайнами, т. е. отрезками кубических парабол (Рисунок 13.4):
Сплайн-интерполяция в Mathcad реализована чуть сложнее линейной. Перед применением функции interp необходимо предварительно определить первый из ее аргументов — векторную переменную s. Делается это при помощи одной из трех встроенных функций тех же аргументов (х, у):
Выбор конкретной функции сплайновых коэффициентов влияет на интерполяцию вблизи конечных точек интервала. Пример сплайн-интерполяции приведен в листинге 13.2.
Полиномиальная регрессия
13.2.2. Полиномиальная регрессияВ Mathcad реализована регрессия одним полиномом, отрезками нескольких полиномов, а также двумерная регрессия массива данных.
Полиномиальная регрессия
Полиномиальная регрессия означает приближение данных (xi,yi) полиномом k-й степени A(x)=a+bx+cx2+dx3+.. .+hxk (Рисунок 13.16). При k=i полином является прямой линией, при k=2 — параболой, при k=3 — кубической параболой и т. д. Как правило, на практике применяются k<5.
Примечание 1
Примечание 1
Для построения регрессии полиномом k-й степени необходимо наличие, по крайней мере, (k+1) точек данных.
В Mathcad полиномиальная регрессия осуществляется комбинацией встроенной функции regress и полиномиальной интерполяции (см. разд. 13.1.2):
ВНИМАНИЕ!
Для построения полиномиальной регрессии после функции regress вы обязаны использовать функцию interp
.
Ввод/вывод в файлы других типов
13.3.2. Ввод/вывод в файлы других типовПодобно вводу/выводу в текстовые файлы можно организовать чтение и запись данных в графические звуковые и файлы.
Графические файлы
При записи и чтении числовой информации в файлы различных графических форматов данные отождествляются с интенсивностью того или иного цвета пиксела изображения, находящегося в файле. Перечислим основные встроенные функции, предназначенные для графического ввода/вывода:
Примечание 1
Примечание 1
Имеется также большое количество функций специального доступа к графическим файлам, например, чтение интенсивности цветов в других цветовых моделях (яркость-насыщенмость-оттенок), а также чтение только одного из основных цветов и т. п. Вы без труда найдете информацию об этих функциях в справочной системе Mathcad, а их применение полностью эквивалентно описанным встроенным функциям.
Действие функций доступа к графическим файлам иллюстрируется листингами 13.19—13.21. Заметим, что для создания изображения используется встроенная функция identity, создающая единичную матрицу. Изображение, созданное листингом 13.19, приведено на Рисунок 13.22.
Другие типы регрессии
13.2.3. Другие типы регрессииКроме рассмотренных, в Mathcad встроено еще несколько видов трехпараметрической регрессии. Их реализация несколько отличается от приведенных выше вариантов регрессии тем, что для них, помимо массива данных, требуется задать некоторые начальные значения коэффициентов а, b, с. Используйте соответствующий вид регрессии, если хорошо представляете себе, какой зависимостью описывается ваш массив данных. Когда тип регрессии плохо отражает последовательность данных, то ее результат часто бывает неудовлетворительным и даже сильно различающимся в зависимости от выбора начальных значений. Каждая из функций выдает вектор уточненных параметров а, b, с.
Примечание 1
Примечание 1
Правильность выбора начальных значений можно оценить по результату регрессии — если функция, выданная Mathcad, хорошо приближает зависимость у (х), значит, они были подобраны удачно.
Пример расчета одного из видов трехпараметрической регрессии (экспоненциальной) приведен в листинге 13.13 и на Рисунок 13.19. В предпоследней строке листинга выведены в виде вектора вычисленные коэффициенты а, ь, с, а в последней строке через эти коэффициенты определена искомая функция f (х) .
Мастер импорта данных и функция READFILE
13.3.3. Мастер импорта данных и функция READFILEВ Mathcad 12 появились две новых, более универсальных, возможности для импорта данных из внешнего файла. Они связаны с появлением Мастера импорта данных, позволяющего осуществить импорт в нужном формате в диалоговом режиме с подсказками, а также новой встроенной функции READFILE, призванной унифицировать процесс импорта. Первый путь позволяет импортировать данные "вручную", проследив процесс считывания данных последовательно, шага за шагом, а второй — автоматизировать его, не путаясь в других многочисленных функциях импорта.
Примечание 1
Примечание 1
Оба способа подразумевают возможность импорта файлов данных самых разных форматов: текстовых с разнообразными символами-разделителями, а также файлы формата xls (Microsoft Excel).
Рассмотрим реализацию первой возможности на примере считывания данных из файла, представленного ранее на Рисунок 13.20 и 13.21 и в листингах 13.16 и 13.17 (см. разд. 13.3.1):
1. Введите команду меню Insert / Component (Вставка / Компонент), а затем выберите в списке тип компонента Data Import Wizard (Мастер импорта данных). В результате появится окно Мастера, которое в пошаговом диалоговом режиме позволит осуществить считывание нужной информации (Рисунок 13.25).
2. Выберите в раскрывающемся списке File format (Формат файла) желаемый формат файла, из которого вы осуществляете импорт (Рисунок 13.26). Если вы затрудняетесь с его точной идентификацией, лучшим решением будет задание типа Delimited Text (Текст с разделителями), что позволит возложить распознавание типа данных и формат их записи на Mathcad.
3. Нажмите кнопку Browse (Пролистать) и отыщите в открывшемся диалоговом окне местоположение нужного вам файла.
4. Если вы уверены (например, основываясь на накопленном опыте работы с Мастером импорта), что данные будут считаны правильно, то можете сразу нажать кнопку Finish (Завершить). Если вы хотите просмотреть и, при необходимости, изменить те или иные опции импорта, нажмите кнопку Next (Вперед). В последнем случае в диалоговом режиме еще на двух страницах Мастера (подобных Рисунок 13.26, но с новыми параметрами) можно будет отредактировать многочисленные параметры импорта (такие как тип разделителя между данными, интервалы импорта и т. д.).
Полиномиальная сплайнинтерполяция
13.1.3. Полиномиальная сплайн-интерполяцияБолее сложный тип интерполяции — так называемая интерполяция В-сплайнами. В отличие от обычной сплайн-интерполяции (см. разд. 13.1.2), сшивка элементарных В-сплайнов производится не в точках хi а в других точках ui, координаты которых предлагается ввести пользователю. Сплайны могут быть полиномами 1, 2 или 3 степени (линейные, квадратичные или кубические). Применяется интерполяция В-сплайнами точно так же, как и обычная сплайн-интерполяция, различие состоит только в определении вспомогательной функции коэффициентов сплайна.
Примечание 1
Примечание 1
Размерность вектора и должна быть на 1, 2 или 3 меньше размерности векторов х и у. Первый элемент вектора и должен быть меньше или равен первому элементу вектора х, а последний элемент и — больше или равен последнему элементу х.
Интерполяция В-сплайнами иллюстрируется листингом 13.3 и Рисунок 13.7.
Регрессия общего вида
13.2.4. Регрессия общего видаВ Mathcad можно осуществить регрессию в виде линейной комбинации C1f1(x)+C2f2(x) + ..., где fi(x) — любые функции пользователя, a Ci — подлежащие определению коэффициенты. Кроме того, имеется путь проведения регрессии более общего вида, когда комбинацию функций и искомых коэффициентов задает сам пользователь.
Приведем встроенные функции для регрессии общего вида и примеры их использования (листинги 13.14 и 13.15), надеясь, что читатель при необходимости найдет более подробную информацию об этих специальных возможностях в справочной системе и Ресурсах Mathcad.
Примечание 1
Примечание 1
Число данных (количество элементов в векторах х и у) должно быть не меньше, чем N. Это менее жесткое требование появилось в Mathcad 12, а до этого для функции регрессии общего вида genfit было необходимо задание не менее N+1 данных.
Сплайнэкстраполяция
13.1.4. Сплайн-экстраполяцияВсе описанные в предыдущих разделах типы интерполяции работают также и как функции экстраполяции данных. Для вычисления экстраполяции достаточно просто указать соответствующее значение аргумента, которое лежит за границами рассматриваемого интервала. С этой точки зрения разницы в применении в Mathcad между интерполяцией и экстраполяцией нет.
На практике при построении экстраполяции следует соблюдать известную осторожность, не забывая о том, что ее успех определяется значимостью ближайших к границе интервала точек. Чем дальше от них вы будете пытаться экстраполировать зависимость, заданную экспериментальными точками, тем сомнительнее будет результат. Сказанное иллюстрируется Рисунок 13.8 и 13.9, на которых изображена линейная (пунктир на обоих графиках) и сплайн- (сплошные кривые) экстраполяция. На Рисунок 13.8 используется линейная сплайн-экстраполяция при помощи функции ispline (см. Рисунок 13.5 в качестве примера интерполяции), а на Рисунок 13.9 — функции кубического сплайна cspline (что соответствует листингу 13.2 и Рисунок 13.4). Видно, что вдали от рассматриваемого интервала результаты экстраполяции совершенно различны, что, конечно, объясняется тем, что она является ни чем иным как параболической зависимостью.
Экстраполяция функцией предсказания
13.1.5. Экстраполяция функцией предсказанияКак мы увидели (см. разд. 13.1.4), стандартные функции интерполяции-экстраполяции стоит применять только в непосредственной близости границ интервала данных. В Mathcad имеется более развитый инструмент экстраполяции, который учитывает распределение данных вдоль всего интервала. В функцию predict встроен линейный алгоритм предсказания поведения функции, основанный на анализе, в том числе осцилляции:
Пример использования функции предсказания на примере экстраполяции осциллирующих данных уj с меняющейся амплитудой приведен в листинге 13.4. Полученный график экстраполяции, наряду с самой функцией, показан на Рисунок 13.10. Аргументы и принцип действия функции predict отличаются от рассмотренных выше встроенных функций интерполяции-экстраполяции. Значений аргумента для данных не требуется, поскольку по определению функция действует на данные, идущие друг за другом с равномерным шагом. Обратите внимание, что результат функции predict вставляется "в хвост" исходных данных.
Многомерная интерполяция
13.1.6. Многомерная интерполяцияДвумерная сплайн-интерполяция приводит к построению поверхности z (х,у), проходящей через массив точек, описывающий сетку на координатной плоскости (х,у). Поверхность создается участками двумерных кубических сплайнов, являющихся функциями (х,у) и имеющих непрерывные первые и вторые производные по обеим координатам.
Многомерная интерполяция строится с помощью тех же встроенных функций, что и одномерная (см. разд. 13.1.2), но имеет в качестве аргументов не векторы, а соответствующие матрицы. Существует одно важное ограничение, связанное с возможностью интерполяции только квадратных NxN массивов данных:
Примечание 1
Примечание 1
Вспомогательные функции построения вторых производных имеют те же матричные аргументы, что и interp: lspline (X, Y), pspline (X, Y), cspline (X, Y).
Пример исходных данных приведен на Рисунок 13.12 в виде графика линий уровня, программная реализация двумерной интерполяции показана в листинге 13.6, а ее результат — на Рисунок 13.13.
Двумерная полиномиальная регрессия (продолжение листинга 13 12)
Рисунок 13.18. Двумерная полиномиальная регрессия (продолжение листинга 13.12)
Обратите внимание, что, если вы
Рисунок 13.21. Файл, созданный листингами 13.16 и 13.18
Обратите внимание, что, если вы выводите данные в файл, пользуясь встроенной функцией WRITEPRN, то в любом случае создается новый текстовый файл. Если даже до записи данных файл с таким именем существовал, то его содержимое будет уничтожено, заменившись новыми данными. Если вы хотите сохранить прежнее содержимое текстового файла с данными, пользуйтесь функцией APPENDPRN. Эта встроенная функция может применяться и для создания нового файла. Иными словами, если файла с заданным именем не существовало, то он будет создан и наполнен теми данными, которые вами определены в документе.
Примечание 3
Примечание 3
Создание нового файла путем использования функции APPENDPRN добавлено разработчиками только в версии Mathcad 11. В прежних версиях программы попытка добавить данные к несуществующему файлу при помощи этой функции вызовет сообщение об ошибке.
созданный листингом
Рисунок 13.20. Файл, созданный листингом 13.16
созданный листингом
Рисунок 13.22. Файл, созданный листингом 13.19
Интерполяция и регрессия
Глава 13. Интерполяция и регрессия| Содержание |
Интерполяция и регрессия
Интерполяция и регрессияПосвятим данную главу самым простым методам обработки данных — интерполяции-экстраполяции и регрессии. Будем считать, что основным объектом исследования будет выборка экспериментальных данных, которые, чаще всего, представляются в виде массива, состоящего из пар чисел (xi,yi) (проблеме ввода/вывода числовых данных во внешние файлы посвящен заключительный раздел этой главы). В связи с этим возникает задача аппроксимации дискретной зависимости y(xi) непрерывной функцией f(x). Функция f (х), в зависимости от специфики задачи, может отвечать различным требованиям:
Различные виды построения аппроксимирующей зависимости f (х) иллюстрирует Рисунок 13.1. На нем исходные данные обозначены кружками, интерполяция отрезками прямых линий — пунктиром, линейная регрессия — наклонной прямой линией, а фильтрация — жирной гладкой кривой. Эти зависимости приведены в качестве примера и отражают лишь малую часть возможностей Mathcad по обработке данных. Вообще говоря, в Mathcad имеется целый арсенал встроенных функций, позволяющий осуществлять самую различную регрессию, интерполяцию и экстраполяцию.
Интерполяция
13.1. ИнтерполяцияДля построения интерполяции-экстраполяции в Mathcad имеется несколько встроенных функций, позволяющих "соединить" точки выборки данных (xi,yi) кривой разной степени гладкости. По определению интерполяция означает построение функции А(Х), аппроксимирующей зависимость у(х) в промежуточных точках (между xi). Поэтому интерполяцию еще по-другому называют аппроксимацией. В точках xi значения интерполяционной функции должны совпадать с исходными данными, т. е. A(xi)=y (xi).
Примечание 1
Примечание 1
Везде в этом разделе, при рассказе о различных типах интерполяции будем использовать вместо обозначения А(х) другое имя ее аргумента A(t), чтобы не путать вектор данных х и скалярную переменную t.
Исходное двумерное поле данных (продолжение листинга 13 6)
Рисунок 13.12. Исходное двумерное поле данных (продолжение листинга 13.6)
Экспоненциальная регрессия (продолжение листинга 13 13)
Рисунок 13.19. Экспоненциальная регрессия (продолжение листинга 13.13)
Экстраполяция при помощи функции предсказания (продолжение листинга 13 4)
Рисунок 13.10. Экстраполяция при помощи функции предсказания (продолжение листинга 13.4)
Квадратичная сплайнэкстраполяция (продолжение листинга 13 2)
Рисунок 13.9. Квадратичная сплайн-экстраполяция (продолжение листинга 13.2)
Линейная интерполяция (продолжение листинга 13 1)
Рисунок 13.2. Линейная интерполяция (продолжение листинга 13.1)
Для построения линейной интерполяции служит встроенная функция linterp (листинг 13.1):
ВНИМАНИЕ!
Элементы вектора х должны быть определены в порядке возрастания, т. е. Xl
Линейная регрессия по методу наименьших
Рисунок 13.15. Линейная регрессия по методу наименьших квадратов и методу медиан (продолжение листингов 13.7 и 13.9)
Линейная регрессия (продолжение
Рисунок 13.14. Линейная регрессия (продолжение листинга 13.7 или 13.8)В Mathcad имеется альтернативный алгоритм, реализующий не минимизацию суммы квадратов ошибок, а медиан-медианную линейную регрессию для расчета коэффициентов а и b (листинг 13.9):
Линейная сплайнэкстраполяция
Рисунок 13.8. Линейная сплайн-экстраполяция
Линейная интерполяция
Листинг 13.1. Линейная интерполяция
Как видно из листинга, чтобы осуществить линейную интерполяцию, надо выполнить следующие действия:
1. Ввести векторы данных х и у (первые две строки листинга).
2. Определить функцию linterp (х, у, t).
3. Вычислить значения этой функции в требуемых точках, например, linterp(х,у,2.4)=3.52, или linterp(х,у,6)=5.9, или постройте ее График, как показано на Рисунок 13.2.
Примечание 1
Примечание 1
Обратите внимание, что функция A(t) на графике имеет аргумент t, а не х. Это означает, что функция A(t) вычисляется не только при значениях аргумента (т. е. в семи точках), а при гораздо большем числе аргументов в интервале (0,6), что автоматически обеспечивает Mathcad. Просто в данном случае эти различия незаметны, т. к. при обычном построении графика функции А(х) от векторного аргумента х (Рисунок 13.3) Mathcad по умолчанию соединяет точки графика прямыми линиями (т. е. скрытым образом осуществляет их линейную интерполяцию).
Полиномиальная регрессия
Листинг 13.10. Полиномиальная регрессия
Регрессия отрезками полиномов
Помимо приближения массива данных одним полиномом имеется возможность осуществить регрессию сшивкой отрезков (точнее говоря, участков, т. к. они имеют криволинейную форму) нескольких полиномов. Для этого имеется встроенная функция loess, применение которой аналогично функции regress (листинг 13.11 и Рисунок 13.17):
Параметр span задает степень сглаженности данных. При больших значениях span регрессия практически не отличается от регрессии одним полиномом (например, span=2 дает почти тот же результат, что и приближение точек параболой).
Регрессия отрезками полиномов
Листинг 13.11. Регрессия отрезками полиномов
СОВЕТ
Регрессия одним полиномом эффективна, когда множество точек выглядит как полином, а регрессия отрезками полиномов оказывается полезной в противоположном случае.
Двумерная полиномиальная регрессия
Листинг 13.12. Двумерная полиномиальная регрессия
Примечание 2
Примечание 2
Обратите внимание на знаки транспонирования в листинге. Они применены для корректного представления аргументов (например, z, в качестве вектора, а не строки).
Экспоненциальная регрессия
Листинг 13.13. Экспоненциальная регрессия
Примечание 2
Примечание 2
Многие задачи регрессии данных различными двухпараметрическими зависимостями у(х) можно свести к более надежной, с вычислительной точки зрения, линейной регрессии. Делается это с помощью соответствующей замены переменных.
Регрессия линейной комбинацией функций пользователя
Листинг 13.14. Регрессия линейной комбинацией функций пользователя
Регрессия общего вида
Листинг 13.15. Регрессия общего вида
Запись матрицы в текстовый файл
Листинг 13.16. Запись матрицы в текстовый файл
Чтение данных из текстового файла в матрицу
Листинг 13.17. Чтение данных из текстового файла в матрицу
Дозапись вектора k в существующий текстовый файл
Листинг 13.18. Дозапись вектора k в существующий текстовый файл
Запись матрицы I в графический файл
Листинг 13.19. Запись матрицы I в графический файл
Кубическая сплайнинтерполяция
Листинг 13.2. Кубическая сплайн-интерполяция
Смысл сплайн-интерполяции заключается в том, что в промежутках между точками осуществляется аппроксимация в виде зависимости A(t)=at3+bt2+ ct+d. Коэффициенты а, b, с, d рассчитываются независимо для каждого промежутка, исходя из значений у* в соседних точках. Этот процесс скрыт от пользователя, поскольку смысл задачи интерполяции состоит в выдаче значения A(t) в любой точке t (Рисунок 13.4).
Чтение из графического файла
Листинг 13.20. Чтение из графического файла
Запись в цветной графический файл
Листинг 13.21. Запись в цветной графический файл="46.gif" >
Звуковые файлы
Начиная с версии Mathcad 2001, появилась возможность записывать и считывать амплитуду акустических сигналов в звуковые файлы с расширением wav.
Использование этих встроенных функций позволяет организовать обработку звука.
Анимация
Во многих случаях самый зрелищный способ представления результатов математических расчетов — это анимация. В Mathcad очень просто создавать анимационные ролики и сохранять их в видеофайлах.
Основной принцип анимации в Mathcad — покадровая анимация. Ролик анимации — это просто последовательность кадров, представляющих собой некоторый участок документа, который выделяется пользователем. Расчеты производятся обособленно для каждого кадра, причем формулы и графики, которые в нем содержатся, должны быть функцией от номера кадра. Номер кадра задается системной переменной FRAME, которая может принимать лишь натуральные значения. По умолчанию, если не включен режим подготовки анимации, FRAME=0.
Рассмотрим последовательность действий для создания ролика анимации, например, демонстрирующего перемещение гармонической бегущей волны. При этом каждый момент времени будет задаваться переменной FRAME.
1. Введите в документ необходимые выражения и графики, в которых участвует переменная номера кадра FRAME. Подготовьте часть документа, которую вы желаете сделать анимацией, таким образом, чтобы она находилась в поле вашего зрения на экране. В нашем примере подготовка сводится к определению функции f (x,t) :=sin(x-t) и созданию ее декартова графика у (х, FRAME) .
2. Выполните команду Tools / Animation / Record (Сервис / Анимация / Запись).
3. В диалоговом окне Record Animation (Анимация) задайте номер первого кадра в поле From (От), номер последнего кадра в поле То (До) и скорость анимации в поле At (Скорость) в кадрах в секунду (Рисунок 13.23).
Импорт данных при помощи универсальной функции READFILE
Листинг 13.22. Импорт данных при помощи универсальной функции READFILE
Интерполяция Всплайнами
Листинг 13.3. Интерполяция В-сплайнами
Экстраполяция при помощи функции предсказания
Листинг 13.4. Экстраполяция при помощи функции предсказания
Как видно из Рисунок 13.11, функция предсказания может быть полезна при экстраполяции данных на небольшие расстояния. Вдали от исходных данных результат часто бывает неудовлетворительным. Кроме того, функция predict хорошо работает в задачах анализа подробных данных с четко прослеживающейся закономерностью (типа Рисунок 13.10), в основном осциллирующего характера.
Экстраполяция при помощи функции предсказания
Листинг 13.5. Экстраполяция при помощи функции предсказания
Двумерная интерполяция
Листинг 13.6. Двумерная интерполяция
Линейная регрессия
Листинг 13.7. Линейная регрессия
Другая форма записи линейной регрессии
Листинг 13.8. Другая форма записи линейной регрессии
Построение линейной регрессии двумя разными методами (продолжение листинга 13 7)

Различие результатов среднеквадратичной и медиан-медианной регрессии иллюстрируется на Рисунок 13.15.
Начало создания анимации
Рисунок 13.23. Начало создания анимации
4. Выделите протаскиванием указателя мыши при нажатой левой кнопке мыши область в документе, которая станет роликом анимации.
5. В диалоговом окне Record Animation (Анимация) нажмите кнопку Animate (Анимация). После этого в окошке диалогового окна Record Animation (Анимация) будут появляться результаты расчетов выделенной области, сопровождающиеся выводом текущего значения переменной FRAME. По окончании этого процесса на экране появится окно проигрывателя анимации (Рисунок 13.24).
6. Запустите просмотр анимации в проигрывателе нажатием кнопки воспроизведения в левом нижнем углу окна проигрывателя.
7. В случае если вид анимации вас устраивает, сохраните ее в виде видеофайла, нажав кнопку Save As (Сохранить как) в диалоговом окне Record Animation (Анимация). В появившемся диалоговом окне Save Animation (Сохранить анимацию) обычным для Windows способом укажите имя файла и его расположение на диске.
8. Закройте диалог Record Animation (Анимация) нажатием кнопки Cancel (Отмена) или кнопки управления его окном.
Обычное построение графика функции от векторной переменной х (продолжение листинга 13 1)
Рисунок 13.3. Обычное построение графика функции от векторной переменной х (продолжение листинга 13.1)
Ошибочное построение графика сплайнинтерполяции (продолжение листинга 13 2)
Рисунок 13.6. Ошибочное построение графика сплайн-интерполяции (продолжение листинга 13.2)
Просмотр созданного ролика анимации
Рисунок 13.24. Просмотр созданного ролика анимации
После того как вы сохранили видеофайл, его можно использовать за пределами Mathcad. Скорее всего, если вы, находясь в обозревателе Windows, дважды щелкнете на имени этого файла, он будет загружен в проигрыватель видеофайлов Windows, и вы увидите его на экране компьютера. Таким образом, запуская видеофайлы в обычном проигрывателе, вы имеете возможность устроить красочную презентацию результатов вашей работы как на своем, так и на другом компьютере.
Примечание 2
Примечание 2
При создании файлов анимации допускается выбирать программу видеосжатия (кодек) и качество компрессии. Делается это с помощью кнопки Options (Опции) в диалоговом окне Record Animation (Анимация).
Работа функции предсказания в случае малого количества данных (продолжение листинга 13 5)
Рисунок 13.11. Работа функции предсказания в случае малого количества данных (продолжение листинга 13.5)
Если данных мало, то предсказание может оказаться бесполезным. В листинге 13.5 приведена экстраполяция небольшой выборки данных (из примеров, рассмотренных в предыдущих разделах). Соответствующий результат показан на Рисунок 13.11 для различных крайних точек массива исходных данных, для которых строится экстраполяция.
Разные задачи аппроксимации данных
Рисунок 13.1. Разные задачи аппроксимации данных
Регрессия отрезками полиномов (продолжение листинга 13 11)
Рисунок 13.17. Регрессия отрезками полиномов (продолжение листинга 13.11)
Двумерная полиномиальная регрессия
По аналогии с одномерной полиномиальной регрессией и двумерной интерполяцией (см. разд. 13.1.5), Mathcad позволяет приблизить множество точек Zi,j(xi,yj) поверхностью, которая определяется многомерной полиномиальной зависимостью. В качестве аргументов встроенных функций для построения полиномиальной регрессии должны стоять в этом случае не векторы, а соответствующие матрицы.
ВНИМАНИЕ!
Для построения регрессии не предполагается никакого предварительного упорядочивания данных (как, например, для двумерной интерполяции, которая требует их представления в виде матрицы NxN). В связи с этим данные представляются как вектор.
Двумерная полиномиальная регрессия иллюстрируется листингом 13.12 и Рисунок 13.18. Сравните стиль представления данных для двумерной регрессии с представлением тех же данных для двумерной сплайн-интерполяции (см. листинг 13.6) и ее результаты с исходными данными (см. Рисунок 13.12) и их сплайн-интерполяцией (см. Рисунок 13.13).
Регрессия полиномами разной степени (коллаж результатов листинга 13 10 для разных k)
Рисунок 13.16. Регрессия полиномами разной степени (коллаж результатов листинга 13.10 для разных k)
Пример полиномиальной регрессии квадратичной параболой приведен в листинге 13.10.
Регрессия
13.2. РегрессияЗадачи математической регрессии имеют смысл приближения выборки данных (xi,yi) некоторой функцией f(x), определенным образом минимизирующей совокупность ошибок |f(xi)-yil. Регрессия сводится к подбору неизвестных коэффициентов, определяющих аналитическую зависимость f (х). В силу производимого действия большинство задач регрессии являются частным случаем более общей проблемы сглаживания данных.
Как правило, регрессия очень эффективна, когда заранее известен (или, по крайней мере, хорошо угадывается) закон распределения данных (xi,yi).
Результат двумерной интерполяции (продолжение листинга 13 6)
Рисунок 13.13. Результат двумерной интерполяции (продолжение листинга 13.6)
Результат импорта данных из файла
Рисунок 13.27. Результат импорта данных из файла
Новая функция READFILE облегчает процесс "программного" считывания данных из файла (листинг 13.22):
Следующая страница окна Data Import Wizard
Рисунок 13.26. Следующая страница окна Data Import Wizard
5. После нажатия кнопки Finish (Завершить) в диалоге Data Import Wizard (Мастер импорта данных) и возвращения на рабочую область документа Mathcad введите в местозаполнитель, появившийся слева от таблицы импортированных данных, желаемое имя переменной. В дальнейших расчетах ее можно будет использовать как обычную матрицу.
Итог работы Мастера показан на Рисунок 13.27. Его первая строка является результатом описанных шагов по считыванию данных в матрицу, а вторая строка показывает вывод этой матрицы в стандартной для Mathcad форме.
Сплайнинтерполяция (продолжение листинга 13 2)
Рисунок 13.4. Сплайн-интерполяция (продолжение листинга 13.2)
Чтобы подчеркнуть различия, соответствующие разным вспомогательным функциям cspline, pspline, ispline, покажем результат действия листинга 13.2 при замене функции cspline в предпоследней строке на линейную ispline (Рисунок 13.5). Как видно, выбор вспомогательных функций существенно влияет на поведение A(t) вблизи граничных точек рассматриваемого интервала (0,6) и особенно разительно меняет результат экстраполяции данных за его пределами.
В заключение остановимся на уже упоминавшейся в предыдущем разделе распространенной ошибке при построении графиков интерполирующей функции (см. Рисунок 13.3). Если на графике, например, являющемся продолжением листинга 13.2, задать построение функции А<Х) вместо A(t), то будет получено просто соединение исходных точек ломаной (Рисунок 13.6). Так происходит потому, что в промежутках между точками вычисления интерполирующей функции не производятся.
Сплайнинтерполяция с выбором коэффициентов линейного сплайна lspline
Рисунок 13.5. Сплайн-интерполяция с выбором коэффициентов линейного сплайна lspline
Стартовая страница окна Data Import Wizard
Рисунок 13.25. Стартовая страница окна Data Import Wizard
Всплайнинтерполяция (продолжение листинга 13 3)
Рисунок 13.7. В-сплайн-интерполяция (продолжение листинга 13.3)
Ввод/вывод данных
13.3. Ввод/вывод данныхЗавершим главу, посвященную интерполяции и регрессии, реализации в Mathcad функции ввода/вывода во внешние файлы, поскольку, как правило, анализ данных чаще всего связан с их импортом из внешних источников (например, результатов эксперимента из текстовых или графических файлов) и экспортом на внешние носители. В большинстве случаев ввод внешних данных в документы Mathcad применяется чаще вывода, поскольку Mathcad имеет гораздо лучшие возможности представления результатов расчетов, чем многие пользовательские программы. Для общения с внешними файлами данных в Mathcad имеет несколько разных способов, причем в но-( вой 12-й версии добавлены новые, намного более удобные опции импорта.
Матричные вычисления в Mathcad
Фурьеспектр действительных данных
14.1.1. Фурье-спектр действительных данныхПреобразование Фурье имеет огромное значение для различных математических приложений, и для него разработан очень эффективный алгоритм, называемый БПФ (быстрое преобразование Фурье). Рассмотрим сначала наиболее типичную для физического эксперимента ситуацию расчета Фурье-спектра действительного сигнала, для которой алгоритм БПФ реализован в нескольких встроенных функциях Mathcad, различающихся нормировками:
ВНИМАНИЕ!
Аргумент прямого Фурье-преобразования, т. е. вектор у, должен иметь ровно 2n элементов (n— целое число). Результатом является вектор с 1+2n-1 элементами. Если число данных не совпадает со степенью 2, то необходимо дополнить недостающие элементы нулями, иначе вместо решения появится сообщение об ошибке.
Встроенная функция вейвлетпреобразования
14.2.1. Встроенная функция вейвлет-преобразованияMathcad имеет одну встроенную функцию для расчета вейвлет-преобразования на основе вейвлетобразующей функции Добеши:
Аргумент функции вейвлет-преобразования, т. е. вектор у, должен так же, как и в преобразовании Фурье, иметь ровно 2n элементов (n — целое число). Результатом функции wave является вектор, скомпонованный из нескольких коэффициентов с двухпараметрического вейвлет-спектра. Использование функции wave объясняется на примере анализа суммы двух синусоид в листинге 14.5, а три семейства коэффициентов вычисленного вейвлет-спектра показаны на Рисунок 14.16.
Встроенные функции для сглаживания ВЧфильтр
14.3.1. Встроенные функции для сглаживания: ВЧ-фильтрВ Mathcad имеется несколько встроенных функций, реализующих различные алгоритмы сглаживания данных:
Все функции имеют в качестве аргумента векторы, составленные из массива данных, и выдают в качестве результата вектор сглаженных данных того же размера. Функция medsmooth предполагает, что данные расположены равномерно.
Примечание 1
Примечание 1
Подробную информацию об алгоритмах, заложенных в функции сглаживания, вы найдете в справочной системе Mathcad в статье Smoothing (Сглаживание), находящейся в разделе Statistics (Статистика). Очень полезные сведения о разных типах фильтрации можно отыскать в Быстрых шпаргалках.
Часто бывает полезным совместить сглаживание с последующей интерполяцией или регрессией. Соответствующий пример приведен в листинге 14.7 для функции supsmooth. Результат работы листинга показан на Рисунок 14.18 (кружки обозначают исходные данные, крестики — сглаженные, пунктирная кривая — результат сплайн-интерполяции). Сглаживание тех же данных при помощи "бегущих медиан" и функции Гаусса с разным значением ширины окна пропускания показаны на Рисунок 14.19 и 14.20 соответственно.
Обратное преобразование Фурье
14.1.2. Обратное преобразование ФурьеДля расчета обратного Фурье-преобразования (восстановления сигнала по имеющемуся действительному спектру) следует использовать следующие встроенные функции (они также реализуют алгоритм БПФ):
Примечание 1
Примечание 1
Аргумент (вектор v) функций, реализующих обратное преобразование Фурье, может быть как действительным, так и комплексным. А вот результат их работы является вектором, составленным из действительных чисел. Если аргумент является N-компонентным вектором, где N=l+2n, то в результате получается в два раза больший вектор из 2 (N-1) =2n+1 компонент.
Результат обратного преобразования Фурье спектра, представленного на Рисунок 14.2 и 14.3, показан в виде кружков на Рисунок 14.5 вместе с исходными данными.
Программирование вейвлетпреобразований
14.2.2. Программирование вейвлет-преобразованийНаряду со встроенной функцией wave Mathcad снабжен пакетом расширения для осуществления вейвлет-анализа. Пакет расширения содержит большое число дополнительных встроенных функций, имеющих отношение к вейвлет-преобразованиям. Обзор пакетов расширения выходит за рамки данной книги, поэтому ограничимся простым упоминанием об этой возможности. Напомним, что дополнительную информацию об использовании данных встроенных функций можно найти в соответствующей электронной книге, которую можно открыть при помощи меню Help / E-Books / Wavelet extension pack (Справка / Электронные книги / Вейвлет-анализ данных).
Помимо встроенной функции вейвлет-спектра Добеши и возможностей пакета расширения Mathcad, возможно непосредственное программирование алгоритмов пользователя для расчета вейвлет-спектров. Оно сводится к аккуратному расчету соответствующих семейств интегралов. Один из примеров такой программы приведен в листинге 14.6, а ее результат на Рисунок 14.17. Анализу подвергается та же функция, составленная из суммы двух гармонических функций, сильно различающихся по частоте. Сам график двухпара-метрического вейвлет-спектра с(а,b) на плоскости (а,b) выведен в виде привычных для вейвлет-анализа линий уровня.
Скользящее усреднение ВЧфильтр
14.3.2. Скользящее усреднение: ВЧ-фильтрПомимо встроенных в Mathcad, существует несколько популярных алгоритмов сглаживания, на одном из которых хочется остановиться особо. Самый простой и очень эффективный метод — это скользящее усреднение. Его суть состоит в расчете для каждого значения аргумента среднего значения по соседним w данным. Число w называют окном скользящего усреднения; чем оно больше, тем больше данных участвуют в расчете среднего, тем более сглаженная кривая получается. На Рисунок 14.21 показан результат скользящего усреднения одних и тех же данных (кружки) с разным окном w=3 (пунктир), w=5 (штрихованная кривая) и w=l5 (сплошная кривая). Видно, что при малых w сглаженные кривые практически повторяют ход изменения данных, а при больших w — отражают лишь закономерность их медленных вариаций.
Преобразование Фурье комплексных данных
14.1.3. Преобразование Фурье комплексных данныхАлгоритм быстрого преобразования Фурье для комплексных данных встроен в соответствующие функции, в имя которых входит литера "с":
Функции действительного преобразования Фурье используют тот факт, что в случае действительных данных спектр получается симметричным относительно нуля, и выводят только его половину (см. разд. 14.1. Г). Поэтому, в частности, по 128 действительным данным получалось всего 65 точек спектра Фурье. Если к тем же данным применить функцию комплексного преобразования Фурье (Рисунок 14.6), то получится вектор из 128 элементов. Сравнивая Рисунок 14.3 и 14.6, можно уяснить соответствие между результатами действительного и комплексного Фурье-преобразования.
Устранение тренда НЧфильтр
14.3.3. Устранение тренда: НЧ-фильтрЕще одна типичная задача возникает, когда интерес исследований заключается не в анализе медленных (или низкочастотных) вариаций сигнала у(х) (для чего применяется сглаживание данных), а в анализе быстрых его изменений. Часто бывает, что быстрые (или высокочастотные) вариации накладываются определенным образом на медленные, которые обычно называют трендом. Часто тренд имеет заранее предсказуемый вид, например, линейный. Чтобы устранить тренд, можно предложить последовательность действий, реализованную в листинге 14.9.
1. Вычислить регрессию f (х), например, линейную, исходя из априорной информации о тренде (предпоследняя строка листинга).
2. Вычесть из данных у(х) тренд f (х) (последняя строка листинга).
Полосовая фильтрация
14.3.4. Полосовая фильтрацияВ предыдущих разделах была рассмотрена фильтрация быстрых вариаций сигнала (сглаживание) и его медленных вариаций (снятие тренда). Иногда требуется выделить среднемасштабную составляющую сигнала, уменьшив как более быстрые, так и более медленные его компоненты. Одна из возможностей решения этой задачи связана с применением полосовой фильтрации на основе последовательного скользящего усреднения.
Алгоритм Полосовой фильтрации приведен в листинге 14.10, а результат его применения показан на Рисунок 14.23 сплошной кривой. Алгоритм реализует такую последовательность операций:
1. Выставление ноль-линии, т. е. приведение массива данных у к нулевому среднему значению путем его вычитания из каждого элемента у (третья и четвертая строки листинга).
2. Устранение из сигнала у высокочастотной составляющей, имеющее целью получить сглаженный сигнал middle, например, с помощью скользящего усреднения с малым окном w (в листинге 14.10 w=3).
3. Выделение из сигнала middle низкочастотной составляющей slow, например, путем скользящего усреднения с большим окном w (в листинге 14.9 w=7), либо с помощью снятия тренда (см. разд. 14.3.3).
4. Вычесть из сигнала middle тренд slow (последняя строка листинга), тем самым выделяя среднемасштабную составляющую исходного сигнала у.
Пример артефакты дискретного Фурьепреобразования
14.1.4. Пример: артефакты дискретного Фурье-преобразованияПри численном нахождении преобразования Фурье следует очень внимательно относиться к таким важнейшим параметрам, как объем выборки (в терминах листинга 14.1, xМАХ) и интервал дискретизации (?). Соотношение этих двух величин определяет диапазон частот (?0,?N), для которых возможно вычисление значений Фурье-спектра. В этой связи хотелось бы обратить внимание на три типичные опасности, которые могут подстерегать неподготовленного исследователя при расчете дискретного Фурье-преобразования и быть для него весьма неожиданными.
Влияние конечности интервала выборки
Во-первых, следует обратить внимание на само определение преобразования Фурье как интеграла с бесконечными пределами. Его численное отыскание подразумевает принципиальную ограниченность интервала интегрирования (просто в силу конечности объема выборки). Поэтому самым очевидным несоответствием будет поиск, вообще говоря, другого интеграла, отличного от интеграла Фурье. Влияние конечности интервала выборки проявляется главным образом на искажении его низкочастотной области. В качестве примера приведем Фурье-спектр гармонической функции с частотой 0.015. Для его расчета достаточно заменить в листинге 14.1 четвертую строку на равенство yi:=sin(2?0.915xi). Соответствующий Фурье-спектр изображен на Рисунок 14.7 (сверху — в обычном, а снизу — в более крупном масштабе) и демонстрирует не совсем правильное поведение в низкочастотной области. Как видно, спектр содержит довольно широкий максимум вместо одиночного пика, как было в случае средних частот сигнала на Рисунок 14.3.
Примечание 1
Примечание 1
Если быть точным, вместо спектра некоторой функции f (х) дискретное преобразование Фурье подразумевает вычисление спектра другой функции f (х)Ф(х), где Ф(х) — это функция-ступенька, равная единице в пределах расчетного интервала и нулю за его пределами. В частотной области это соответствует операции свертки означенных двух функций, что, конечно, искажает (неизвестный) точный спектр f (х). Для борьбы с проявлением конечности интервала выборки используются специальные методы, основанные на применении техники спектральных окон.
Примечание 2
Примечание 2
Из сказанного ясно, почему исследователя не должна смущать необходимость дополнения массива исходных данных нулями до размера 2" (чтобы можно было использовать алгоритм БПФ). По самому определению дискретного Фурье-преобразования, исходная функция и так предполагается равной нулю за пределами расчетного интервала, что приводит к неминуемым искажениям.
Пример спектр модели сигнал/шум
14.1.5. Пример: спектр модели сигнал/шумПока мы использовали в качестве примера детерминированный сигнал, представляющий собой сумму трех синусоид. Несмотря на единство термина "дискретное преобразование Фурье", прикладное применение спектрального анализа можно довольно четко разделить на две категории.
Фурье-спектр
Внесем минимальное добавление в расчеты листинга 14.1, а именно добавим к его четвертой строке (в которой определяется yi) еще одно (четвертое) слагаемое: псевдослучайную величину ?rnd(1), где значение 1/? характеризует отношение сигнал/шум. Явный вид изменений, которые следует внести в листинг 14.1, приведен на Рисунок 14.10, наряду с графиком сигнала у(х). Расчет Фурье-спектра данного сигнала (в соответствии с алгоритмом, представленным выше, см. листинг 14.1) показан на Рисунок 14.11. Как видно, присутствие шумовой компоненты может значительно искажать спектр сигнала и затруднять его интерпретацию.
Примечание 1
Примечание 1
Максимальное значение спектра на левом крае частотного интервала является ни чем иным, как проявлением искажающего влияния конечности выборки и сдвига ноль-линии (см. разд. 14.1.3), произошедшим из-за внесения шума с математическим ожиданием, равным 0.5.
Спектральная фильтрация
14.3.5. Спектральная фильтрацияАльтернативой всем представленным до сих пор алгоритмам, в частности, методу полосовой фильтрации (см. предыдущий разд.) является фильтрация на основе интеграла Фурье. Пока мы использовали для подавления в сигнале тех или иных частотных диапазонов определенные процедуры, основанные на арифметических преобразованиях. Между тем, для той же цели (правда, с большими вычислительными затратами) можно применять методы Фурье-анализа. Действительно, если вычислить спектр сигнала, удалить из него (или существенно уменьшить) определенные частоты, а затем выполнить обратное преобразование Фурье, то результатом будет фильтрованный сигнал.
Пример фильтрации на основе преобразования Фурье приведен в листинге 14.11. В качестве модельного сигнала использовалась смесь двух гармонических сигналов и равномерно распределенного шума (Рисунок 14.24). Фурье-спектр данных z, вычисленный при помощи встроенной функции fft, показан на Рисунок 14.25.
Двумерный спектр Фурье
14.1.6. Двумерный спектр ФурьеВ Mathcad имеется возможность применять встроенные функции комплексного преобразования Фурье не только к одномерным, но и к двумерным массивам, т. е. матрицам. Соответствующий пример приведен в листинге 14.4 и на Рисунок 14.14 в виде графика линий уровня исходных данных и рассчитанного Фурье-спектра.
Пример вычисление спектра мощности
14.3.6. Пример: вычисление спектра мощностиЗавершим главу, посвященную спектральному анализу, еще одним примером вычисления спектра мощности (его определение приведено в разд. 14.1.5) модельного сигнала, связанного с использованием его Фурье-спектра. В листинге 14.12 сначала выполняется преобразование Фурье суммы гармонического сигнала и шума (распределенного равномерно), а затем (в трех последних строках листинга) производится его сглаживание путем скользящего усреднения с окном, равным 3. Из курса математической статистики известно, что квадратная степень сглаженного преобразования Фурье может считаться оценкой спектра мощности, и для вычисленного таким образом спектра уже можно применять вероятностные оценки погрешностей его отсчетов. Результаты работы листинга 14.12 (Фурье- и энергетический спектр) показаны на Рисунок 14.26 (кружками и сплошной кривой соответственно).
Адаптивное сглаживание (продолжение листинга 14 7)
Рисунок 14.18. Адаптивное сглаживание (продолжение листинга 14.7)
Автокорреляционная функция модельной зависимости сигнал / шум (продолжение листинга 14 3)
Рисунок 14.12. Автокорреляционная функция модельной зависимости сигнал / шум (продолжение листинга 14.3)
Данные (слева) и их Фурьеспектр (справа) (продолжение листинга 14 4)
Рисунок 14.14. Данные (слева) и их Фурье-спектр (справа) (продолжение листинга 14.4)
Фурьепреобразование модельного сигнала и его спектр мощности (продолжение листинга 14 12)
Рисунок 14.26. Фурье-преобразование модельного сигнала и его спектр мощности (продолжение листинга 14.12)
Фурьепреобразование модельного сигнала и спектральное окнофильтр (продолжение листинга 14 11)
Рисунок 14.25. Фурье-преобразование модельного сигнала и спектральное окно-фильтр (продолжение листинга 14.11)
Фурьеспектр суммы гармонического сигнала и константы (влияние конечности выборки)
Рисунок 14.8. Фурье-спектр суммы гармонического сигнала и константы (влияние конечности выборки)
Фурьеспектр
14.1. Фурье-спектрИнтегральные преобразования массива сигнала у(х) ставят в соответствие всей совокупности данных у(х) некоторую функцию другой координаты F(v). Рассмотрим встроенные функции для расчета интегральных преобразований, реализованных в Mathcad.
Математический смысл преобразования Фурье состоит в представлении сигнала у(х) в виде бесконечной суммы синусоид вида F(v)sin(v-x). Функция F(V) называется преобразованием Фурье, или интегралом Фурье, или Фурье-спектром сигнала. Ее аргумент v имеет смысл частоты соответствующей составляющей сигнала. Обратное преобразование Фурье переводит спектр F(V) в исходный сигнал у (х). Согласно определению,
Как видно, преобразование Фурье является комплексной величиной, даже если сигнал действительный.
Спектральный анализ
Глава 14. Спектральный анализ| Содержание |
График Фурьеспектра данных (продолжение листинга 14 1)
Рисунок 14.3. График Фурье-спектра данных (продолжение листинга 14.1)
График Фурьеспектра данных
Рисунок 14.11. График Фурье-спектра данных
Спектр мощности
В силу стохастичности исходных данных, представляющих сумму полезного сигнала и шума, сами вычисленные значения спектра Фурье носят также случайный характер. В этой связи необходимо знать, с какой погрешностью они рассчитываются. Однако из курса математической статистики известно, что для обычного Фурье-преобразования случайного сигнала (в частности, нормального) не найдено оценок для погрешности. Это слабое место Фурье-спектров делает их практически неприменимыми для анализа случайных сигналов, а вместо них надо применять так называемые спектры мощности (или, по-другому, энергетические спектры), для которых указанные оценки существуют.
Не углубляясь в теорию математической статистики, приведем пример вычисления спектра мощности сигнала (Рисунок 14.10), основанный его определении. Как известно, спектром мощности сигнала называют Фурье-преобразование его корреляционной функции. Таким образом, алгоритм расчета спектра мощности сводится к следующему: во-первых, вычислению автокорреляционной функции (Рисунок 14.12); во-вторых, ее прореживанию и (или) сглаживанию (в целях уменьшения влияния конечности выборки); и, наконец, в-третьих, расчету ее Фурье-преобразования. Результат вычисления спектра мощности (листинг 14.3) в соответствии с приведенным сценарием показан на Рисунок 14.13.
Примечание 2
Примечание 2
Аналогичным образом, через Фурье-преобразование взаимной корреляционной функции, определяются взаимные спектры мощности двух выборок.
Примечание 3
Примечание 3
Еще один способ вычисления спектров мощности, не требующий расчета функции корреляции, приведен ниже (см. разд. 14.3.6).
Примечание 4
Примечание 4
Методика расчета в Mathcad корреляционной функции случайного процесса обсуждалась в главе 12 (см. разд. 12.3.3).
График спектра мощности данных модельной зависимости сигнал / шум (продолжение листинга 14 3)
Рисунок 14.13. График спектра мощности данных модельной зависимости сигнал / шум (продолжение листинга 14.3)
Иллюстрация влияния конечности выборки на расчет низкочастотной части Фурьеспектра
Рисунок 14.7. Иллюстрация влияния конечности выборки на расчет низкочастотной части Фурье-спектра
Сдвиг ноль-линии
Еще одним, наиболее ярким, проявлением вредного влияния конечности интервала выборки может служить расчет Фурье-преобразования суммы гармонического сигнала и константы (Рисунок 14.8). Для того чтобы получить данный рисунок, достаточно еще слегка (по сравнению с Рисунок 14.7) модифицировать строку листинга, касающуюся определения компонент вектора у, добавив к нему i: yi:=sin(2?0.915xi)+1.
Сравнивая Рисунок 14.7 и 14.8, несложно догадаться, почему так разительно изменился вид спектра в низкочастотной области. Пугающий рост спектра на левом крае частотного интервала объясняется совокупностью двух факторов: конечности выборки и добавлением к сигналу ненулевой постоянной составляющей (так называемым сдвигом ноль-линии). Сумма сигнала и константы определяет соответствующее влияние на вычисленный спектр, который также оказывается (просто в силу линейности операции интегрирования) суммой спектров сигнала и ступенчатой функцией (равной той самой константе внутри расчетного интервала и нулю за его пределами).
Исходные модельные данные (продолжение листинга 14 1)
Рисунок 14.1. Исходные модельные данные (продолжение листинга 14.1)
Чтобы смысл преобразования Фурье был более понятен, используем в качестве модельных данных дискретизацию детерминированного сигнала,, равного сумме трех синусоид (Рисунок 14.1).
Исходный модельный сигнал (кружки)
Рисунок 14.24. Исходный модельный сигнал (кружки) и результат его фильтрации на основе Фурье-преобразования (продолжение листинга 14.11)
Комплексное преобразование Фурье (продолжение листинга 14 2)
Рисунок 14.6. Комплексное преобразование Фурье (продолжение листинга 14.2)
Быстрое преобразование Фурье
Листинг 14.1. Быстрое преобразование Фурье
демонстрирует расчет
Листинг 14.1 демонстрирует расчет Фурье-спектра по N=128 точкам, причем предполагается, что интервал дискретизации данных yi равен ?. В середине листинга применяется встроенная функция fft, а его оставшаяся часть предназначена для корректного пересчета соответствующих значений частот ?i (они вычисляются в последней строке листинга). Обратите внимание, что результаты расчета представляются в виде модуля Фурье-спектра (Рисунок 14.2), поскольку сам спектр является комплексным. Очень полезно сравнить полученные амплитуды и местоположение пиков спектра (Рисунок 14.3) с определением синусоид в листинге 14.1.Полосовая фильтрация
Листинг 14.10. Полосовая фильтрация
Фильтрация на основе преобразования Фурье
Листинг 14.11. Фильтрация на основе преобразования Фурье
отличается от листинга
Листинг 14.11 отличается от листинга 14.1 (см. разд. 14.1. Г), главным образом, последними двумя строками, в которых, собственно, и определяется явный вид спектрального фильтра w(?), или, по-другому, спектральное окно. Обратное Фурье-преобразование спектра произведения спектрального окна w(?) и Фурье-спектра сигнала z, представляющее собой результат фильтрации, показано на Рисунок 14.24 сплошной кривой.Вычисление спектра мощности
Листинг 14.12. Вычисление спектра мощности
Комплексное быстрое преобразование Фурье (продолжение листинга 14 1)
Листинг 14.2. Комплексное быстрое преобразование Фурье (продолжение листинга 14.1)
Расчет спектра мощности для модели сигнал/шум
Листинг 14.3. Расчет спектра мощности для модели сигнал/шум
Двумерное преобразование Фурье
Листинг 14.4. Двумерное преобразование Фурье
Поиск вейвлетспектра Добеши
Листинг 14.5. Поиск вейвлет-спектра Добеши
Поиск вейвлетспектра на основе "мексиканской шляпы"
Листинг 14.6. Поиск вейвлет-спектра на основе "мексиканской шляпы"
Примечание 1
Примечание 1
Программа листинга очень проста, но исключительно далека от хорошей в смысле быстродействия. Каждый интеграл вычисляется независимо, без использования методов ускорения, типа применяемых в алгоритме БПФ. Однако простые приемы программирования вполне доступно раскрывают математический смысл вейвлет-преобразования.
Сглаживание с последующей сплайнинтерполяцией
Листинг 14.7. Сглаживание с последующей сплайн-интерполяцией
Сглаживание скользящим усреднением
Листинг 14.8. Сглаживание скользящим усреднением
Примечание 1
Примечание 1
Приведенная программная реализация скользящего усреднения самая простая, но не самая лучшая. Возможно, вы обратили внимание, что все кривые скользящего среднего на Рисунок 14.21 слегка "обгоняют" исходные данные. Почему так происходит, понятно: согласно алгоритму, заложенному в последнюю строку листинга 14.8, скользящее среднее для каждой точки вычисляется путем усреднения значений предыдущих w точек. Чтобы результат скользящего усреднения был более адекватным, лучше применить центрированный алгоритм расчета по w/2 предыдущим и w/2 последующим значениям. Он будет немного сложнее, поскольку придется учитывать недостаток точек не только в начале (как это сделано в программе с помощью функции условия if), но и в конце массива исходных данных.
Устранение тренда
Листинг 14.9. Устранение тренда
На Рисунок 14.22 показаны исходные данные (кружками), выделенный с помощью регрессии линейный тренд (сплошной прямой линией) и результат устранения тренда (пунктир, соединяющий крестики).
Матрицарезультат вычисления Фурьеспектра данных (продолжение листинга 14 1)
Рисунок 14.2. Матрица-результат вычисления Фурье-спектра данных (продолжение листинга 14.1)
Исключительно важными представляются два параметра, заданные в предпоследней строке листинга 14.1, называемые соответственно граничной частотой и частотой Найквиста. Граничная частота ?0 определяет нижнюю, а частота Найквиста ?N — верхнюю границу аргумента вычисленного спектра, как показано маркерами на Рисунок 14.3. Кроме того, важно, что интервал дискретизации Фурье-спектра также равен ?0, а общее число вычисляемых точек спектра составляет N/2 (в нашем примере N/2=64). Последние утверждения иллюстрируются маркерами на Рисунок 14.4, изображающем график Фурье-спектра вблизи нижней границы частот.
Модель сигнал / шум
Рисунок 14.10. Модель сигнал / шум
Низкочастотная область Фурьеспектра (продолжение листинга 14 1)
Рисунок 14.4. Низкочастотная область Фурье-спектра (продолжение листинга 14.1)
Обратное преобразование Фурье (продолжение листинга 14 1)
Рисунок 14.5. Обратное преобразование Фурье (продолжение листинга 14.1)
Видно, что в рассматриваемом случае сигнал у(х) восстановлен с большой точностью, что характерно для плавного изменения сигнала. Если же в качестве аргумента функции ifft использовать модуль Фурье-спектра (мы присвоили этому вектору в листинге 14.1 имя а), то профиль исходного сигнала будет реконструирован правильно, но окажется сдвинутым на определенное расстояние вдоль оси х. Так происходит из-за того, что взятие абсолютной величины комплексного спектра уничтожает информацию об относительной фазе отсчетов данных.
Расчеты Фурьеспектров гармонических сигналов с разной частотой ("маскировка частот")
Рисунок 14.9. Расчеты Фурье-спектров гармонических сигналов с разной частотой ("маскировка частот")
Избавиться от искажений, вызванных сдвигом ноль-линии, довольно просто. Достаточно (до Фурье-преобразования) вычислить среднее значение выборки и затем вычесть его из каждого элемента выборки. Если после этой операции вычислить Фурье-спектр, то он окажется примерно таким, как показано на Рисунок 14.7.
Маскировка частот
Еще один классический пример ошибочного расчета Фурье-спектра связан с возможным присутствием в сигнале гармоник с частотой, превышающей частоту Найквиста, в данном примере ?N=0.б4 (см. разд. 14.1.Г). Иллюстрация эффекта, называемого "маскировкой частот", приведена на Рисунок 14.9, который содержит расчет спектров трех различных синусоидальных сигналов с разной частотой f0. Первый спектр сигнала с частотой, меньшей частоты Найквиста, вычислен верно, а вот два остальных спектра показывают, что в случае превышения частоты Найквиста в спектре начинают присутствовать "лишние" пики. Появление артефактов спектра связано с тем, что дискретных отсчетов начинает не хватать для того, чтобы прописать высокочастотные гармоники с достаточной информативностью.
Примечание 3
Примечание 3
Напоминаем, что все листинги, имеющиеся в книге, а также документы Mathcad с расчетами всех рисунков, вы найдете на прилагаемом компакт-диске.
Результат полосовой фильтрации (продолжение листинга 14 10)
Рисунок 14.23. Результат полосовой фильтрации (продолжение листинга 14.10)
Сглаживание "бегущими медианами"
Рисунок 14.19. Сглаживание "бегущими медианами"
Сглаживание и фильтрация
14.3. Сглаживание и фильтрацияПри анализе данных часто возникает задача их фильтрации, заключающаяся в устранении одной из составляющих зависимости y(xi). Наиболее часто целью фильтрации является подавление быстрых вариаций y(xi), которые чаще всего обусловлены шумом. В результате из быстроосциллирующей зависимости y(xi) получается другая, сглаженная зависимость, в которой доминирует более низкочастотная составляющая. В связи с этим считают, что сглаживание является частным случаем более общей задачи фильтрации сигнала, которая связана с подавлением определенных диапазонов частот спектра.
Наиболее простым и эффективным рецептом сглаживания (smoothing) можно считать регрессию различного вида, разобранную в предыдущей главе (см. разд. 13.2). Однако регрессия часто уничтожает информативную составляющую данных, оставляя лишь наперед заданную пользователем зависимость.
Часто рассматривают противоположную задачу фильтрации — устранение медленно меняющихся вариаций в целях исследования высокочастотной составляющей. В этом случае говорят о задаче устранения тренда. Иногда интерес представляют смешанные задачи выделения среднемасштабных вариаций путем подавления как более быстрых, так и более медленных вариаций. Одна из возможностей решения связана с применением полосовой фильтрации.
Несколько примеров программной реализации различных вариантов фильтрации приведены в данном разделе.
Сглаживание при помощи функции ksmooth
Рисунок 14.20. Сглаживание при помощи функции ksmooth
Скользящее усреднение с разными w=3 5 15 (листинг 14 8 коллаж трех графиков)
Рисунок 14.21. Скользящее усреднение с разными w=3, 5,15 (листинг 14.8, коллаж трех графиков)
Чтобы реализовать в Mathcad скользящее усреднение, достаточно очень простой программы, приведенной в листинге 14.8. Она использует только значения у, оформленные в виде вектора, неявно предполагая, что они соответствуют значениям аргумента х, расположенным через одинаковые промежутки. Вектор х применялся лишь для построения графика результата (Рисунок 14.21).
Спектральный анализ
Спектральный анализМощным инструментом обработки данных, определенных дискретной зависимостью y(xi) или непрерывной функцией f(x) (полученной, например посредством интерполяции или регрессии, как об этом рассказано в главе 13), является спектральный анализ, имеющий в своей основе различные интегральные преобразования. Спектральный анализ используется как в целях подавления шума, так и для решения других проблем обработки данных. Спектром совокупности данных у(х) называют некоторую функцию другой координаты (или координат) F(w), полученную в соответствии с определенным алгоритмом. Примерами спектров являются преобразование Фурье (см. разд. 14. Т) и вейвлет-преобразование (си. разд. 14.2). Напомним, что некоторые преобразования, например, Фурье и Лапласа, можно осуществить в режиме символьных вычислений (см. главу 5). Каждое из интегральных преобразований эффективно для решения своего круга задач анализа данных.
Задачами, непосредственно связанными со спектральным анализом, являются проблемы сглаживания и фильтрации данных (см. разд. 14.3). Они заключаются в построении для исходной экспериментальной зависимости y(xi) некоторой (непрерывной или дискретной) зависимости f (х), которая должна приближать ее, учитывая к тому же, что данные (xi,yi) получены с некоторой погрешностью, выражающей шумовую компоненту измерений. При этом функция f (х) с помощью того или иного алгоритма уменьшает погрешность, присутствующую в данных (xi,yi). Такого типа задачи называют задачами фильтрации. Сглаживание путем построения регрессии данных (см. разд. 13.2) — это частный случай фильтрации.
Примечание 1
Примечание 1
Здесь нельзя не отметить, что в качестве дополнений к Mathcad поставляются три пакета расширения, включающие большое количество дополнительных встроенных функций для обработки данных. Названия пакетов расширения говорят сами за себя: Wavelet extension pack (Вейвлет-анализ данных), Signal processing (Анализ сигналов) и Image processing (Анализ изображений). Работа с пакетами расширения не слишком отличается от обычных приемов работы с Mathcad — следует только установить их определенным образом, как описано в руководстве пользователя Mathcad. Рассмотрение их возможностей выходит далеко за пределы данной книги, поэтому лишь упомянем о том, что описание встроенных функций и примеры их применения рассмотрены в трех электронных книгах, представляющих, соответственно, три упомянутых пакета расширения. Получить доступ к нужной книге можно, наводя указатель мыши на пункт E-Books (Электронные книги) в меню Help (Справка) и выбирая в открывающемся подменю имя нужного пакета расширения.
Сравнение синусоиды и вейвлетобразующей функции
Рисунок 14.15. Сравнение синусоиды и вейвлетобразующей функции
Из-за своего математического смысла вейвлет-спектр имеет не один аргумент, а два. Помимо частоты, вторым аргументом ь является место локализации вейвлетобразующей функции. Поэтому b имеет ту же размерность, что и х.
Устранение тренда (продолжение листинга 14 9)
Рисунок 14.22. Устранение тренда (продолжение листинга 14.9)
Вейвлетспектр на основе функции Добеши (продолжение листинга 14 5)
Рисунок 14.16. Вейвлет-спектр на основе функции Добеши (продолжение листинга 14.5)
Вейвлетспектр на основе "мексиканской шляпы" (продолжение листинга 14 6)
Рисунок 14.17. Вейвлет-спектр на основе "мексиканской шляпы" (продолжение листинга 14.6)
Вейвлетспектры
14.2. Вейвлет-спектрыВ последнее время возрос интерес к другим интегральным преобразованиям, в частности, вейвлет-преобразованию (или дискретному волновому). Оно применяется, главным образом, для анализа нестационарных сигналов и для многих задач подобного рода оказывается более эффективным, чем преобразование Фурье. Основным отличием вейвлет-преобразования является разложение данных не по синусоидам (как для преобразования Фурье), а по другим функциям, называемым вейвлетобразующими. Вейвлетобразующие функции, в противоположность бесконечно осциллирующим синусоидам, локализованы в некоторой ограниченной области своего аргумента, а вдали от нее равны нулю или ничтожно малы. Пример такой функции, называемой "мексиканской шляпой", показан на Рисунок 14.15.
Биржевой анализ: Технический анализ - Инструменты - Софт
- Биржевой анализ - Технический анализ
- Обучение техническому анализу
- Индексы технического анализа
- Индикаторы технического анализа
- Методы технического анализа
- Графика в техническом анализе
- Технический анализ рынков
- Российский технический анализ
- Инструменты технического анализа
- Математика в биржевом анализе
- Разновидности биржевого анализа
- Mathematica в биржевом анализе
- MathCAD в биржевом анализе
- Maple в биржевом анализе
- Matlab в биржевом анализе