Система Mathcad

- n - количество частей, на которые разбивается отрезок [a, b] при решении системы. Чем больше это число, тем точнее получается решение, но расчет занимает больше времени.

- F(x,y) - векторная функция, элементы которой содержат правые части уравнений системы в нормальной форме (когда левые части - первые производные от соответствующих функций, а в правых частях производные отсутствуют). Аргументами этой функции являются вектор y, элементы которого соответствуют различным неизвестным функциям системы, и скалярный аргумент x , соответствующий независимой переменной в системе. В случае одного уравнения функция F может быть скалярной функцией, зависящей от двух скалярных переменных x и y.

Возвращаемым значением всех вышеперечисленных встроенных функций является матрица. Первый столбец этой матрицы - это точки, на которые разбивается отрезок [a, b], а остальные столбцы - это значения функций системы в этих точках. Если в аргументе функции rkfixed было указано количество частей n = 100, то матрица будет содержать 101 строку вместе с начальной.

Решение систем обыкновенных дифференциальных уравнений.

Для численного интегрирования систем ОДУ в MathCAD также имеется выбор - либо использовать вычислительный блок Given/Odesolve, либо встроенные функции rkfixed, Rkadapt и Bulstoer.

При решении систем ОДУ MathCAD требует, чтобы система ОДУ была представлена в нормальной форме (когда левые части - первые производные от соответствующих функций, а в правых частях производные отсутствуют):

где Y и Y' - соответствующие неизвестные векторные функции переменной t, а F - вектор правых частей системы уравнений первого порядка. Именно векторное представление используется для ввода системы ОДУ в среде MathCAD.

Если в систему ОДУ входят и уравнения высших порядков, то оно тоже сводится к системе уравнений первого порядка, как было показано выше. При этом количество нулевых условий для вычислительного блока Given/Odesolve, а также размер вектора начальных условий y и размер вектора правых частей F(x,y) для встроенных функций rkfixed, Rkadapt и Bulstoer должны быть равны сумме порядков всех уравнений.

Вначале покажем решение систем ОДУ первого порядка с использованием вычислительного блока Given/Odesolve

Функция Odesolve для системы ОДУ имеет несколько иной, по сравнению с одним уравнением, синтаксис. Теперь она возвращает вектор функций, составляющих решение системы. Поэтому в качестве первого аргумента функции нужно ввести вектор, состоящий из имен функций, использованных при вводе системы. Второй и третий аргументы то же самое, что и в задаче с одним ОДУ.

Решение системы ОДУ показано на графике слева. Как известно, решения ОДУ часто удобнее изображать не в таком виде, а в фазовом пространстве, по каждой из осей которого откладываются значения каждой из найденных функций (как показано на рисунке справа). При этом аргумент входит в них лишь параметрически. В рассматриваемом случае двух ОДУ такой график - фазовый портрет системы - является кривой на фазовой плоскости. В общем случае, если система состоит из N ОДУ, то фазовое пространство является N - мерным. При N > 3 наглядность теряется, и для визуализации фазового портрета приходится строить его различные проекции.

Рассмотрим решение этой же системы ОДУ первого порядка с использованием встроенной функции rkfixed.

Полученное решение полностью соответствует вышеприведенному решению с использованием вычислительного блока Given/Odesolve. Следует отметить, что начальные условия здесь задаются в виде вектора y, а функциям x(t) и y(t) соответствуют элементы этого вектора y1 и y2. Вектор начальных условий y и вектор правых частей F имеют размер равный двум, т.к. система состоит из двух уравнений первого порядка. Для системы ОДУ, состоящей из двух уравнений второго порядка, размер этих векторов будет равен четырем

Решение жестких ОДУ и систем ОДУ .

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

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

- Radau - метод Radaus для жестких систем. Полностью аналогичен использованию функции odesolve с выбранным в контекстном меню методом Stiff.

- Stiffb - метод Булирша-Штера, адаптированный для жестких систем.

- Stiffr - метод Розенброка.

Краевые задачи для ОДУ

Постановка краевых задач для ОДУ отличается от задач Коши, рассмотренных выше, тем, что граничные условия для них ставятся не в одной начальной точке, а на обеих границах расчетного интервала. Если имеется система N ОДУ первого порядка, то часть из N условий может быть поставлена на одной границе интервала, а оставшиеся условия - на противоположной границе. В связи с тем, что условия поставлены не на одной, а на обеих границах интервала, краевые задачи нельзя решить изложенными выше методами, предназначенными для задач Коши. Для решения краевой задачи в MathCAD нет отдельной функции. Однако есть функции, позволяющие превратить краевую задачу в задачу Коши. Эти функции «угадывают» недостающие начальные условия, исходя из того, что решение должно удовлетворять заданным условиям в конечной точке интервала интегрирования. Простейшей из функций, предназначенных для приведения краевой задачи к задаче Коши, является функция sbval. Для того, чтобы решить двухточечную краевую задачу с помощью этой функции, следует выполнить следующие действия:

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

2. Задайте функцию F(x,y). Эта функция уже описывалась выше. Она представляет собой вектор, каждый элемент которого - это правая часть одного из уравнений системы.

3. Задайте еще одну векторную функцию load(x,v). Это функция от скалярного аргумента x и вектора v, который имеет столько же компонент, сколько недостающих начальных условий в системе. Сам вектор load должен содержать такое же количество элементов, как и вектор F, т.е. столько, сколько должно быть начальных условий в задаче. Если начальное значение какой-либо из функций известно, то соответствующий элемент вектора load должен содержать это значение. Для функций, начальное значение которых неизвестно, соответствующий элемент вектора load должен содержать один из элементов вектора v.

4. Следует задать еще одну некоторую функцию score(x,y). Аргументы этой функции - скаляр x и вектор y, который имеет столько элементов, сколько уравнений в системе. Количество компонент вектора score должно равняться количеству граничных условий, заданных в конечной точке отрезка интегрирования. На самом деле каждая компонента этого вектора задает одно из граничных условий в конечной точке. Например, если в задаче есть граничной условие yi(b)=c, то один из элементов вектора score должен быть функцией, которая обращается в нуль при значениях x=b и y(b)=c. Конкретный вид этой функции не играет особой роли, поэтому проще всего задавать ее в таком виде: scorek(x,y):=yi-c. Таким же образом должны быть заданы все элементы вектора score для всех конечных условий задачи.

5. Теперь все введенные величины нужно использовать как аргументы в функции sbval. Использование этой функции выглядит следующим образом: Y:=sbval(v,a,b,F,load,score). Аргументы a и b - это начало и конец отрезка интегрирования.

6. Результатом функции sbval будет вектор, содержащий недостающие начальные значения. Их последовательность задается той последовательностью, в которой были использованы компоненты вектора v в функции load.

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

8. Теперь можно решить полученную задачу как задачу Коши, с помощью, например, функции Rkadapt.

Ниже приведено решение дифференциального уравнения на отрезке [0,3] с начальным условием y(0)=1 и граничным условием y(3)=4. Недостающее начальное условие - y'(0).

Решение уравнений в частных производных.

Постановка задач для уравнений в частных производных включает определение самого уравнения (или системы нескольких уравнений), а также необходимого количества краевых условий (число и характер которых зависит от специфики уравнения). Согласно своему названию, уравнения должны содержать частные производные неизвестной функции u (или нескольких функций, если уравнений несколько) по различным аргументам, например, пространственной переменной x и времени t. Соответственно, для решения задачи требуется вычислить функцию нескольких переменных, например, u(x,t) в некоторой области определения аргументов 0 x L и 0 t T. Граничные условия определяются как заданные временные зависимости функции u, или производных этой функции, на границах расчетной области 0 и L, а начальные - как заданная функция u(x,0).

Сами уравнения в частных производных (несколько условно) можно разделить на три основных типа:

- параболические - содержащие первую производную по одной переменной и вторую - по другой, причем все эти производные входят в уравнение с одинаковым знаком;

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

- эллиптические - содержащие только вторые производные, причем одного знака.

Средства MathCAD позволяют решать одномерные параболические и гиперболические уравнения (с одной пространственной и одной временной переменной), а также двумерное уравнение Пуассона. Такой, казалось бы, узкий круг решаемых задач на самом деле охватывает подавляющее большинство задач, возникающих в физике и технике.

Дифференциальные уравнения параболического типа.

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

1. Ключевое слово given

2. Уравнение, которое нужно решить. Уравнение должно иметь такой вид: . Для ввода производных в данном случае нельзя пользоваться обычным оператором производной, а нужно пользоваться нижним индексом, как это обычно делается в литературе для записи уравнений в частных производных. При этом нижний индекс набирается не как числовой, а как буквенный (точка).

3. Граничные условия для функции u(x,t). Если уравнение второго порядка по x, то граничных условий должно быть два. Можно использовать как граничные условия Дирихле (u(x0,t) = w(t)), так и граничные условия Неймана (ux(x0,t) = w(t)) или их комбинацию, как в вышеприведенном примере.

4. Начальное значение для неизвестной функции - u(x,0).

5. функция pdsolve(u,x,xrange,t,trange,xpts,tpts). Ее аргументы имеют следующее назначение

- u - имя функции, относительно которой решается уравнение. Для системы уравнений здесь должен быть вектор имен функций (как в odesolve).

- X- имя пространственной переменной.

- Xrange - двухкомпонентный вектор, задающий начало и конец интервала изменения пространственной переменной.

- T - имя временной переменной. Основная разница между пространственной и временной переменными в данном случае - это то, что все уравнения могут содержать только первые производные по временной переменной.

- Trange - еще один двухкомпонентный вектор. Этот вектор задает начало и конец временного интервала, на котором решается задача.

- Xpts, tpts - количество точек, разбивающих для интегрирования пространственный и временной интервалы соответственно. Эти два параметра можно не указывать, тогда количество точек будет выбрано автоматически из соображений достаточной точности. Рекомендуется задавать эти параметры во всех задачах, кроме простых, поскольку во многих случаях высокая точность вычислений теряет смысл из-за погрешности, вносимой самим методом.

Выше приведен пример решения одномерной задачи теплопроводности для одного бруска, один конец которого теплоизолирован, а ругой поддерживается при определенной температуре

Дифференциальные уравнения гиперболического типа.

Функция pdsolve также позволяет решать системы ДУ в частных производных первого порядка по времени. Такая возможность может быть использована для решения задач с ДУ гиперболического типа. Ведь поскольку уравнения гиперболического типа содержат вторую производную по времени, то они не могут быть напрямую введены для решения функцией pdsolve. ДУ гиперболического типа должно быть приведено к системе из двух уравнений первого порядка по времени (как это делалось ранее для обычных ДУ высоких порядков). А далее полученная задача может быть решена с помощью функции pdsolve как система уравнений.

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

Эллиптические уравнения.

Решение эллиптических уравнений в частных производных реализовано только для единственного типа задач - двумерного уравнения Пуассона. Это уравнение содержит вторые производные функции u(x,y) по двум пространственным переменным. Для его решения в MathCAD есть отдельные встроенные функции multigrid и relax.

Вопросы

1. Какие дифференциальные уравнения называются ОДУ первого порядка? Высшего порядка? Что такое нормальная форма ОДУ первого и высшего порядка? К чему сводятся ОДУ высшего порядка при решении?

2. Можно ли решить дифференциальные уравнения в MathCADе символьно?

3. Как решаются ОДУ с помощью вычислительного блока Given/Odesolve? Какой метод решения реализует функция Odesolve? Как можно изменить метод решения для этой функции?

4. Как решаются ОДУ с помощью встроенной функции rkfixed? Чем функция rkfixed отличается от функции Rkadapt?

5. Как осуществляется решение системы ОДУ с помощью вычислительного блока Given/Odesolve? Приведите примеры.

6. Как осуществляется решение системы ОДУ с помощью функции rkfixed? Приведите примеры.

Лекция № 16

Программирование в MathCAD

MathCAD - это система, ориентированная на пользователя, который не обязан знать абсолютно ничего о программировании. Создатели MathCAD изначально поставили перед собой такую задачу, чтобы дать возможность профессионалам - математикам, физикам и инженерам самостоятельно проводить сложные расчеты, не обращаясь за помощью к программистам. Несмотря на блестящее воплощение этих замыслов, выяснилось, что совсем без программирования MathCAD серьезно теряет в своей силе, в основном, из-за недовольства пользователей, знакомых с техникой создания программ.

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

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

Для вставки программного кода в документы MathCAD имеется специальная панель инструментов

.

Большинство кнопок этой панели выполнено в виде текстового представления операторов программирования, поэтому их смысл легко понятен. Операторы программирования могут быть введены только с этой панели, но никак не с клавиатуры (можно еще их вводить с помощью сочетаний клавиш, которые приведены в тексте всплывающей подсказки).

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

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

Для определения функции это будет выглядеть следующим образом:

Внутри программы можно использовать глобальные переменные документа, но изменить их значение внутри программы никак нельзя. Можно создать в программе другие переменные, доступ к которым может осуществляться только из самой программы. Эти переменные называются локальными переменными. Локальные переменные «не видны» извне. Локальная переменная создается с помощью знака локального присвоения с панели Programming. Для оператора локального присваивания, так же как и для операторов обычного и глобального присваивания, можно изменить внешний вид так, чтобы он выглядел как обычный знак равенства. Для этого достаточно вызвать контекстное меню этого оператора и в нем выбрать команду View Definition As/Equal.

Последняя строка любой программы не должна содержать никаких управляющих операторов. Эта строка задает значение, возвращаемое программой. В качестве этого значения может быть скаляр (число или переменная), вектор или матрица. Таким образом, последняя строка программы может содержать имя локальной переменной либо некоторое математическое выражение, куда входят как локальные, так и глобальные переменные, либо вектор или матрицу.

Вставить строку программного кода в уже созданную программу можно в любой момент с помощью той же самой кнопки Add Line. Для этого следует поместить на нужное место внутри программы линии ввода.

Если вертикальная линия ввода находится в начале текущей строки, то нажатие кнопки приведет к появлению новой строки перед текущей строкой, если вертикальная линия ввода находится в конце строки, то нажатие кнопки приведет к появлению новой строки после текущей. Можно выделить только часть текущей строки, что повлияет на положение новой строки в программе (например, если необходимо для некоторого условия выполнить не одно действие, а несколько). Таким образом, основной принцип создания программных модулей заключается в правильном расположении строк кода. Ориентироваться в их действии довольно легко, т.к. фрагменты кода одного уровня сгруппированы в программе с помощью вертикальных черт.

Проверка условий в программах.

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

Для проверки условий в программах MathCAD служит оператор if . В поле ввода справа нужно ввести условие. Для ввода условий служит панель Boolean. На этой панели есть кнопки, предназначенные для проверки условий , а также кнопки, предназначенные для вставки логических операций (логическое отрицание, операция «и», операция «или», операция «исключающее или»). Вторые позволяют создавать сложные условия. В поле ввода слева нужно ввести строку программы, которая должна выполняться, если введенное условие истинно.

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

.

Если невыполнение условия должно привести к выполнению какого-либо иного программного кода, можно в строке, следующей за оператором if, вставить оператор . В поле ввода слева от этого оператора необходимо ввести строку программы, которая будет выполняться только в том случае, если не выполнилось условие, заданное в операторе if.

Следует помнить, что если в программе введено подряд несколько строк с оператором if

,

то выражение слева от оператора otherwise будет выполнено только в том случае, если не выполняются условия, заданные во всех операторах if.

Примеры:

Создание циклов.

Другим очень важным элементом при создании нелинейных программ являются циклы. Циклы позволяют повторять несколько раз выполнение одного и того же программного блока. Для создания циклов в MathCAD предусмотрено два оператора for и while.

Цикл for - цикл со счетчиком.

В таких циклах создается некоторая переменная-счетчик, значение которой изменяется после каждого выполнения тела цикла. Выход из цикла происходит при достижении этой переменной заданного значения. Этот оператор вводится с панели Programming клавишей for . В поле ввода после слова for следует указать имя переменной - счетчика. Это может быть любое имя, которое не использовалось ранее в программе. Внутри цикла можно использовать эту переменную в любых выражениях, нельзя только присваивать ей никакого значения. В поле ввода после знака следует указать диапазон значений переменной-счетчика. Вводить диапазон в данном случае следует так же, как и при создании ранжированной переменной. Вместо диапазона в данном поле ввода можно указать имя некоторого массива (вектора или матрицы). В этом случае переменная-счетчик будет последовательно принимать значения всех элементов этого массива. Возможность перебора элементов массива не может быть реализована с помощью цикла while, поэтому именно в таких случаях цикл for и является незаменимым. В поле ввода под словом for следует ввести тело цикла.

Пример: Заполнить вектор числами от xнач до xкон с шагом h. Затем определить сумму элементов этого вектора и найти их среднее арифметическое значение.

Цикл while - цикл, который выполняется до тех пор, пока выполняется определенное условие.

В поле ввода справа от слова while следует ввести условие. Это условие строится по тем же правилам, что и в операторе if. Оно будет проверяться после каждого выполнения тела цикла и в тот момент, когда условие перестанет выполняться, повторение тела цикла прекратится. В поле ввода ниже слова while следует ввести тело цикла (напомним, что для ввода нескольких строк в теле цикла надо воспользоваться кнопкой Add Line).

Пример1: Вычислить сумму с точностью е.

Пример2: Найти первый элемент, превышающий определенный порог.

Использование операторов break и continue.

Иногда возникает необходимость повлиять на выполнение цикла некоторым образом, например, прервать его выполнение по какому-либо условию или выполнять некоторые итерации не так, как другие. Для этого и служат операторы break и continue.

Оператор break, если он расположен внутри цикла, означает немедленное прекращение выполнения текущей итерации и выход из цикла. Если есть необходимость прекратить выполнение цикла по какому-либо условию, то следует использовать конструкцию следующего вида break if (условие).

Пример: Выделить из массива все элементы от начала и до первого вхождения в него заданного числа.

Оператор continue используется для того, чтобы немедленно перейти в начало цикла и начать следующую итерацию. Этот оператор также обычно используется в составе конструкции вида continue if (условие). Оператор continue используется в случаях, когда необходимо чтобы некоторые вычисления производились для одних итераций и не производились для других.

Пример: Требуется заполнить элементы квадратной матрицы в шахматном порядке.

Оператор return (возврат значения).

Как мы уже указывали, результат выполнения программного модуля помещается, как правило, в последней его строке. Но можно прервать выполнение программы в любой ее точке (например, с помощью условного оператора) и выдать некоторое значение, применив оператор return. В этом случае при выполнении указанного условия значение, введенное в поле ввода после return, возвращается в качестве результата, и никакой другой код больше не выполняется. Вставляется в программу оператор return с помощью одноименной кнопки панели Programming . Пример:

Обработка ошибок.

Система MathCAD предоставляет пользователю некоторый контроль над ошибками, которые могут возникнуть при вычислении выражений или при выполнении программ. Для этой цели служит оператор on error, который можно вставить с помощью одноименной кнопки панели Programming . В поле ввода справа следует ввести выражение или программу, которые необходимо вычислить (известно, что это выражение может содержать ошибку при определенных значениях входных параметров). В поле ввода слева следует ввести выражение, которое будет выполнено вместо правого выражения, если при выполнении последнего возникнет ошибка. Пример: Если аргументу функции присвоено нулевое значение, то в программе возникает ошибка - деление на нуль. Но за счет оператора on error сообщение не выводится, а функции в этой точке присваивается значение, указанное слева от оператора on error - значение машинной бесконечности.

В поле ввода слева может быть введено текстовое выражение, сообщающее об ошибке

Конечно, если программа содержит только простейшие математические операции, то никаких ошибок, кроме деления на нуль или превышения наибольшего допустимого числа в ней возникнуть не может. В таких программах оператор on error используется редко. Но для более сложных программ, которые содержат функции решения дифференциальных уравнений, аппроксимации или другие сложные функции MathCAD, использование on error может избавить от многих трудностей.

Иногда может возникнуть ситуация обратная той, которая была описана выше, т.е. необходимо, чтобы при определенных условиях результатом выражения было сообщение об ошибке, хотя в действительности при этом не возникает ни одной стандартной ошибки MathCAD. Для таких случаев в MathCAD предусмотрена встроенная функция error. В качестве аргумента этой функции нужно в кавычках указать текст сообщения об ошибке, который должен быть выведен. Таким образом, если необходимо, чтобы программа возвращала ошибку при определенном условии, то следует использовать конструкцию вида: error («текст ошибки») if (условие).

Пример:

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

Примеры программирования.

1. Вычислить функцию sin(x) с точностью е.

2. Даны массивы А(5) и В(5). Получить массив С, в который записаны сначала элементы

массива А в порядке возрастания, а затем элементы массива В порядке убывания

3. По введенным значениям коэффициентов А, В, С определить корни квадратного уравнения

4. Дан массив натуральных чисел В(10). Определить, есть ли в нем 4 последовательных числа (например, 1, 2, 3, 4, и т.п.). Напечатать такие группы чисел.

Вопросы

1. Какая панель служит для вставки программного кода в документ MathCAD? Можно ли операторы программирования набрать с клавиатуры?

2. С какой команды начинается создание программного блока? Как с ее помощью можно создавать разветвленный программный блок?

3. Что такое определение программного блока? Обращение к программному блоку?

4. Что такое глобальные и локальные переменные для программного блока? Что может содержать последняя строка программного блока?

5. Как работает оператор if в программном блоке? Приведите пример.

6. Создание цикла с параметром в программном блоке. Приведите пример.

7. Создание цикла while в программном блоке. Приведите пример.

8. Для чего служат операторы break, continue в программном блоке? Приведите примеры.

9. Как работает оператор return в программном блоке? Приведите пример.

10. Как осуществляется обработка ошибок в программном блоке? Приведите пример.

Лекция №3

(установочная лекция)

Численные методы решения задач.

Обработка экспериментальных данных средствами MathCAD

Развитие электронной вычислительной техники, создание алгоритмических языков программирования и обширного математического обеспечения ЭВМ позволяет широко использовать численные методы вычислительной математики при решении различного рода прикладных задач в науке, технике, производстве.

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

Численными методами часто приходится решать следующие математические задачи:

1. решение нелинейных (алгебраических и трансцендентных) уравнений;

2. вычисление определенных интегралов;

3. решение обыкновенных дифференциальных уравнений;

4. решение дифференциальных уравнений в частных производных;

5. решение задач оптимизации;

6. обработка массивов числовых данных.

Каждая из этих задач может представлять собой самостоятельную прикладную задачу или являться составной частью более сложных прикладных задач.

Решение нелинейных уравнений

Обычно нелинейные уравнения делят на трансцендентные и алгебраические. Нелинейные уравнения, содержащие тригонометрические функции или другие специальные функции, например, lg(x) или ex, называются трансцендентными. Методы решения нелинейных уравнений такого типа делятся на аналитические и численные.

Первые позволяют найти решение непосредственно с помощью формул и всегда обеспечивают получение точного решения. Известным примером такого рода является формула корней квадратного уравнения. В численных методах задается процедура решения в виде многократного применения некоторого алгоритма. Задача отыскания корней нелинейного уравнения f(x) = 0 считается решенной, если мы сумеем определить корни с нужной степенью точности.

Для решения нелинейных уравнений известны следующие численные методы: метод половинного деления (метод дихотомии), метод хорд, метод касательных (Ньютона), метод секущих, метод простой итерации. Рассмотрим метод половинного деления.

Графическая интерпретация метода показана на рис.1.

Рис.1 Графическая интерпретация метода половинного деления

В этом методе отыскание корня уравнения f(x) = 0 проходит в два этапа. На первом этапе необходимо отделить корень, т.е.выделить интервал на оси абсцисс, на котором функция f(x) меняет свой знак. Для отделения корня следует провести вычисление функции f(x) в точках, расположенных через равные интервалы по оси x, до тех пор, пока не будут найдены два последовательных значения функции f(xn) и f(xn+1), имеющие противоположные знаки.

На втором этапе производится уточнение корня. Найденный интервал [xn, xn+1], содержащий корень, делится пополам

Затем по разности знаков функции на концах интервала определяем, на каком из полученных двух интервалов находится корень уравнения. Найденный интервал снова делится пополам и т.д.. В результате интервал, на котором находится корень сужается. Процесс повторяется до тех пор, пока f(xср) не станет достаточно близким к нулю. Блок-схема алгоритма метода показана на рис.2.

Рис.2 Блок-схема алгоритма метода половинного деления

Численное интегрирование

К численному интегрированию обращаются, когда нельзя через элементарные функции аналитически записать первообразную интеграла

или когда подобная запись имеет сложный вид.

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

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

1. Методы Ньютона - Котеса. Эти методы требуют, чтобы значения x были заданы с постоянным шагом. Они основаны на полиномиальной аппроксимации подынтегральной функции. Алгоритмы методов просты и легко поддаются программной реализации.

2. Методы Гаусса - Кристоффеля - методы наивысшей алгебраической точности. Эти методы используют неравно отстоящие узлы, расположенные по алгоритму, обеспечивающему минимальную погрешность интегрирования. Требуют большего объема памяти, чем методы первой группы.

3. Методы Монте-Карло. В этих методах узлы выбираются с помощью датчика случайных чисел. Ответ носит вероятностный характер. Методы эффективны при вычислении большой кратности.

Рассмотрим методы первой группы. По этому методу интервал разбивается на n равных отрезков, длина каждого из которых . Для вычисления каждой элементарной площади подынтегральную функцию на данном отрезке заменяют с некоторой степенью точности более простой функцией, интеграл от которой можно вычислить, используя только ординаты на концах отрезка. Метод Ньютона -Котеса дает три формулы для приближенного вычисления определенного интеграла

1. При замене f(x) на каждом отрезке прямой, параллельной оси x (рис.3,а), получим формулу прямоугольников

2. При замене f(x) на каждом отрезке прямой, соединяющей ординаты концов отрезка (рис.3,б), получим формулу трапеций

3. При замене f(x) дугой параболы, проведенной через концы трех ординат (рис.3,в), получим формулу Симпсона. При этом число отрезков должно быть четным

Точность вычисления интеграла тем выше, чем больше n и меньше h, и в пределе при n h 0 указанные формулы дадут точную величину определенного интеграла.

Рис.3 Численное определение интеграла

Блок схема алгоритма численного определения интеграла с заданной степенью точности представлена на рис. 4. На блок-схеме a,b,e,n - исходные данные, е - заданная точность, n - число разбиений, i1 - начальное значение интеграла, его можно задать равным нулю, I - вычисленное значение интеграла.

Для достижения заданной точности число разбиений удваивают до тех пор, пока будет удовлетворено условие |i - i1|<e.

Рис.4 Блок схема алгоритма численного определения интеграла с заданной степенью точности по формуле прямоугольников

Обработка экспериментальных данных средствами MathCAD.

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

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

Интерполяция.

Интерполяцией называют заполнение отрезками кривых промежутков между заданными точками по тому или иному закону. Для проведения интерполяции в первую очередь должна быть задана экспериментальная зависимость в виде набора точек на плоскости. Для этого должны быть заданы два одномерных массива (вектора) - vx и vy, содержащие соответственно значения координат x и y каждой точки. При этом важно, чтобы значения в векторе vx были заданы в порядке возрастания.

Система MathCad позволяет проводить линейную интерполяцию и сплайн-интерполяцию наборов экспериментальных точек.

Простейшим вариантом интерполяции является линейная интерполяция. Она заключается в простом соединении точек между собой отрезками прямых. Для реализации такой интерполяции в MathCad существует встроенная функция linterp(vx,vy,x) , где vx vy - уже известные векторы, содержащие координаты последовательности точек, x - координата точки, в которой нужно вычислить значение интерполирующей функции. Пример построения линейной интерполяции приведен на рис.5

Рис.5 Линейная интерполяция

На практике линейная интерполяция применяется редко.

Из всех видов интерполяции наиболее часто используется интерполяция, где экспериментальные точки попарно соединяются отрезками полиномов. Чаще всего для этого выбирают полиномы третьей степени(поэтому такая кривая и называется кубическим сплайном). Для того чтобы найти коэффициенты этих полиномов, очевидно, недостаточно того условия, что кривая должна проходить через экспериментальные точки. Поэтому на сплайн накладываются дополнительные условия сшивки - первая и вторая производные слева и справа от каждой экспериментальной точки должны быть равны между собой. Но и после этого количество условий остается на два меньше, чем количество неизвестных коэффициентов. Дополнительные два условия должны быть наложены в начальной и конечной экспериментальных точках, поскольку в них нет условий сшивки. Эти условия можно выбрать по-разному. В MaqthCad существуют три различных функции для построения кубических сплайнов с различными дополнительными условиями.

· lspline(vx,vy) - в начальной и конечной точках накладывается условие линейности, т.е. вторая производная от функции равна нулю. Первая буква в названии функции - l , означает linear (линейный).

· pspline(vx,vy) - на первом и последнем интервале кривая является параболой, т.е. полиномиальный коэффициент при x3 равен нулю. Буква p означает parabolic (параболический).

· cspline(vx,vy) - полиномиальные коэффициенты при x3 на первых двух интервалах равны между собой точно так же, как на последних двух интервалах. Буква c означает cubic (кубический).

Результатом каждой из перечисленных функций является вектор, содержащий значения вторых производных от интерполяционной кривой во всех точках, заданных в массиве vx. Для того чтобы исходя из этого вектора построить кривую, нужно воспользоваться встроенной функцией interp(v,vx,vy,x), где vx и vy - массивы экспериментальных точек, v - массив, полученный как результат одной из трех функций, перечисленных выше, x- координата, в которой нужно вычислить значение интерполяционной кривой. Пример интерполяции кубическим сплайном приведен на рис.6.

Рис.6 Интерполяция кубическим сплайном

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

Кубический сплайн является эффективным средством построения интерполяционной кривой в подавляющем большинстве случаев. Но иногда использование кубического сплайна может привести к нежелательным результатам. Чаще всего это происходит в тех случаях, когда данные очень неравномерно распределены вдоль оси x. В таких случаях на кривой кубического сплайна могут появляться острые экстремумы в промежутках между экспериментальными точками. В некоторых подобных случаях получить лучщую интерполяционную кривую помогает использование другого вида интерполяции - В-сплайна. Основное отличие В-сплайна от всех описанных выше методов - сшивка отрезков кривых происходит не в экспериментальных точках, а между ними, в специально заданных точках. В MathCad для реализации интерполяции В-сплайном служит функция bspline(vx.vy,u,n), где vx,vy - векторы, содержащие координаты экспериментальных точек, u - вектор, содержащий координаты точек сшивки, n - порядок полинома. Результатом функции bspline является вектор, который далее следует использовать как аргумент функции interp. В-сплайн в MathCad можно построить из отрезков прямых, парабол или кубических парабол, т.е. допустимые значения параметра n - 1,2 или 3. Количество точек сшивки не является произвольной величиной и должно быть всегда на n-1 меньше, чем количество экспериментальных точек. Также на координаты точек сшивки накладываются следующее условие: первая точка сшивки должна быть не правее первой экспериментальной точки, а последняя не левее последней экспериментальной точки. Остальные точки сшивки могут произвольным образом располагаться внутри отрезка. Пример использования В-сплайна приведен на рис.7.

Экстраполяция.

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

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

Рис.7 Сравнение эффективности кубического сплайна и В-сплайна

Также в MathCad есть одна встроенная функция для проведения экстраполяции - predict(v,m,n), где v - вектор значений функции на том отрезке, где она известна (значения аргумента в данной функции не задаются, и считается, что точки распределены равномерно), m - количество элементов вектора v, на основании которых проводится экстраполяция ( естественно, выбираются точки, ближайшие к правой границе), n - количество точек в рассчитываемом векторе. Пример экстраполяции от осциллирующей функции приведен на рис.8

k,101+i

Рис.8 Экстраполяция с помощью функции predict

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

Регрессия

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

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

Линейная регрессия

Линейная регрессия является наиболее простой, но, тем не менее, используется чаще любого другого вида регрессии. Она заключается в нахождении таких значений параметров a и b, чтобы прямая y = a+bx наилучшим образом аппроксимировала заданной набор точек. Для проведения линейной регрессии по методу наименьших квадратов в MathCad существует функция line(vx,vy). Результатом функции line будет вектор, содержащий значения параметров a и b для построения регрессионной прямой. Пример линейной регрессии представлен на рис.9

Рис.9 Линейная регрессия с помощью функции line

На рис.9 для анализа эффективности линейной регрессии построена псевдоэкспериментальная последовательность точек. Для этого к точным значениям линейной функции прибавлены случайные числа, сгенерированные с помощью функции rnorm.Затем проведена линейная регрессия полученного набора точек. Ее результат можно сравнить с исходной прямой. Как видно из рис., параметры регрессионной прямой заметно отличаются от исходных, но на графике обе прямые проходят достаточно близко.

Для линейной регрессии в MathCad реализован также метод медиан с помощью функции medfit(vx,vy). Результатом этой функции является вектор, аналогичный результату line. Нельзя утверждать, что один из двух методов регрессии более точен. Метод наименьших квадратов является наиболее универсальным, поэтому функция line считается в MathCad основной функцией для проведения линейной регрессии.

Полиномиальная регрессия.

Кроме аппроксимации неизвестной функции с помощью прямой, широкое применение находит и аппроксимация с помощью полиномов различной степени. Для этой цели в MathCad существует функция regress(vx,vy,n). Последний аргумент данной функции задает степень полинома. Можно использовать полином любой степени, но не большей, чем число точек в выборке минус один. При n = 1 получится линейная регрессия. На практике наибольшее применение находит полиномиальная регрессия от второй до пятой степени.

Результатом функции regress является вектор, содержащий коэффициенты аппроксимирующего полинома. Эти коэффициенты располагаются в векторе, начиная с четвертого элемента в порядке возрастания степеней. Первые три элемента данного вектора являются служебными и используются для того, чтобы результат функции regress можно было использовать как первый аргумент функции interp по аналогии со сплайн-интерполяцией. Пример полиномиальной регрессии представлен на рис.10.

Рис.10 Полиномиальная регрессия

На рис.10 для иллюстрации полиномиальной регрессии построена псевдоэкспериментальная последовательность точек. В качестве теоретической функции был использован полином третьей степени с коэффициентами 0,1,-2,1. Как видно из примера, коэффициенты, рассчитанные функцией regress, значительно отличаются от коэффициентов исходного полинома. Тем не менее в области экспериментальных точек обе кривые достаточно близки, но за пределами этой области резко расходятся.

Аппроксимация набора точек различными элементарными функциями

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

· expfit(vx,vy,vg) - регрессия экспоненциальной функцией y = a*eb*x+c.

· sinfit(vx,vy,vg) - регрессия синусоидальной функцией y = a*sin(x+b)+c.

· pwrfit(vx,vy,vg) - регрессия степенной функцией e = a*xb +c.

Перечисленные функции используют трехпараметрическую аппроксимирующую функцию, нелинейную по параметрам. При вычислении оптимальных значений трех параметров регрессионной функции по методу наименьших квадратов возникает необходимость в решении сложной системы из трех нелинейных уравнений. Такая система часто может иметь несколько решений. Поэтому в функциях MathCad, которые проводят регрессию трехпараметрическими зависимостями, введен дополнительный аргумент vg. Данный аргумент - это трехкомпонентный вектор, содержащий приблизительные значения параметров a,b и c, входящих в аппроксимирующую функцию. Неправильный выбор элементов вектора vg может привести к неудовлетворительному результату регрессии. На рис.11 приведен пример проведения экспоненциальной регрессии с помощью функции expfit, регрессия проведена для двух различных значений вектора vg.

Рис.11 Экспоненциальная регрессия

Аппроксимация набора точек функциями, заданными пользователем.

В MathCad существуют средства для проведения регрессии самого общего вида. Это означает, что можно использовать любые функции в качестве аппроксимирующих и находить оптимальные значения любых их параметров, как линейных, так и нелинейных. В том случае, если регрессионная функция является линейной по всем параметрам, т.е. представляет линейную комбинацию жестко заданных функций, провести регрессию можно с помощью встроенной функции linfit(vx,vy,F). Аргумегт F - это векторная функция, из элементов которой должна быть построена линейная комбинация, наилучшим образом аппроксимирующая заданную последовательность точек. Результатом работы функции linfit является вектор линейных коэффициентов. Каждый элемент этого вектора - коэффициент при функции, стоящей на соответствующем месте в векторе F. Таким образом, для того чтобы получить регрессионную функцию, достаточно скалярно перемножить эти два вектора. Пример такой аппроксимации представлен на рис.12.

Рис.12 Использование функции linfit

Сглаживание

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

· medsmooth(vy,n) - реализует алгоритм «бегущих» медиан. Параметр n определяет ширину «окна» сглаживания, т.е. количество точек, которые используются при вычислении сглаженного значения в каждой точке. Параметр n должен быть целым нечетным числом. Данный алгоритм лучше всего подходит для сглаживания наборов точек, в которых лишь некоторые точки резко выбиваются из общей гладкой последовательности (см.рис.13).

· ksmooth(vx,vy,b) - алгоритм Гауссового ядра. В данном алгоритме сглаженное значение в каждой точке вычисляется как весовое среднее от всего набора данных с ядром в виде функции Гаусса. Параметр b - это параметр ширины функции ядра. Данный метод наилучшим образом подходит для фильтрации зашумленного сигнала (см.рис.14).

· supsmooth(vx,vy) - в данном алгоритме значение в каждой точке заменятся на значение регрессионной прямой, построенной с использованием некоторого количества близлежащих точек, причем данное количество выбирается по-разному для каждой точки с помощью адаптивного алгоритма.

Рис.13 Сглаживание сигнала с узкими нерегулярностями

Рис.14 Сглаживание зашумленного сигнала

Дискретное преобразование Фурье

Еще одна операция, широко применяемая при обработке сигналов , - вычисление Фурье-спектра. Для проведения данной операции с сигналами, заданными в виде массивов данных, система MathCad содержит численный алгоритм, называемый быстрым преобразованием Фурье (Fast Fourier Transform FFT). Для реализации данного алгоритма в MathCad существует несколько различных функций.

· В том случае, если набор данных v состоит из 2m элементов, а также все числа в наборе данных действительны, то для проведения Фурье-преобразования следует пользоваться функциями fft(v) или FFT(v). Обе эти функции выполняют Фурье-преобразование, а различие между ними заключается лишь в нормировке результата. Результатом функций fft и FFT будет массив из 2m -1+1 комплексных чисел.

· Если какое-либо из условий, перечисленных выше, не выполняется, т.е. массив содержит комплексные числа или не может быть расширен до размерности 2m, то следует пользоваться функциями cfft(v) и CFFT(v). Результатом функций cfft и CFFT будет массив комплексных чисел той же размерности, что и исходный.

Результатом каждой из перечисленных функций будет массив комплексных чисел. Исходя из этого массива, можно построить амплитудно-частотную (АЧХ) или фазово-частотную (ФЧХ) характеристику сигнала. Для построения АЧХ следует вычислить абсолютное значение каждого элемента в массиве. Пример такого построения приведен на рис.15.

Рис.15 Использование функции fft для построения АЧХ - сигнала

На рис 16 проведена фильтрация сигнала, построенного на рис.15. Конечно, такой метод фильтрации можно использовать, только если уровень шума ощутимо ниже уровня сигнала.

Рис.16 Амплитудная фильтрация сигнала

Страницы: 1, 2, 3, 4