Главная
страница 1

Системы символьных вычислений.

1. Теоретическое занятие


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

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

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

Напротив, знаменитые великие вычисления ХIХ века содержат большое количество манипуляций с формулами. Наиболее известным является, несомненно, расчет Леверье орбиты Нептуна, который был основан на возмущениях орбиты Урана и привел к открытию Нептуна. Наиболее впечатляющие вычисления с карандашом и бумагой выполнены также в области астрономии: Делоне потребовалось 10 лет для вычисления орбиты Луны и еще 10 лет для ее проверки. Проверка выявила несколько ошибок в знаках и несколько "потерянных" двоек. Результат не является численным, поскольку он состоит в основном из формулы, которая сама по себе занимает 128 страниц 4-й главы его книги. Заметим, что с использованием пакета Mathematica на персональном компьютере типа Pentium 166 проверка данных вычислений, проведенная по архивам Делоне, занимает около 10 минут (не учитывая время на набор формул).

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

Очевидно, численные расчеты не исключают алгебраических вычислений: написание простейших программ требует переписывания формул, на которых основан алгоритм. И мощности компьютеров далеко не решают всех проблем: расчет развития атмосферных процессов с точностью, необходимой метеорологам для 48-часового прогноза, на наиболее мощных компьютерах, доступных в настоящее время (типа CRAY), требует гораздо более 48 часов. Увеличение производительности в 10 раз едва ли улучшит ситуацию, поскольку одновременно с уточнением теории потребуется в 2 раза уменьшить размер сетки (чтобы отличить, например, погоду в Петербурге от погоды в Петергофе), а это потребует 8-кратного увеличения объема вычислений.


Компьютерная алгебра.


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

В настоящее время появились хорошо работающие системы, такие как Maple, Mathematica, Macsyma, Derive, Axiom и некоторые другие. Все упомянутые выше системы, так же как и большинство неупомянутых систем, являются весьма дружественными по отношению к пользователю. Конечно же, и синтаксис языка пользователя у них различный, и библиотеки доступных функций могут меняться от нескольких сотен до тысяч, и внутренние структуры и даже используемые алгоритмы значительно отличаются друг от друга, но все они обладают общими свойствами. Таких принципиальных общих свойств значительно больше, чем различий и, таким образом, после освоения одной из систем компьютерной алгебры переход к другой системе не является сложной проблемой.

Разработка, развитие и даже использование этих систем постепенно выделились в автономную научную дисциплину, относящуюся, очевидно, к информатике. Ее цели лежат в области искусственного интеллекта, несмотря на то, что методы все более и более удаляются от нее. Кроме того, используемые алгоритмы вводят в действие все менее элементарные математические средства. Таким образом, эта дисциплина лежит на стыке нескольких областей, что одновременно обогащает ее и делает более трудной в исследовательском плане. Наименование этой дисциплины в последнее время стабилизировалось как "Calcul formel" во французском языке и "Computer algebra" - в английском. В русском языке используют термины "компьютерная алгебра", "символьные и алгебраические вычисления", " аналитические вычисления" и другие.

Использование систем компьютерной алгебры.


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

Однако, в отличие от языка программирования типа Фортран, в котором синтаксические тонкости требуют тщательного изучения, в то время как принципы работы компилятора можно полностью игнорировать, здесь пользователь должен очень быстро разобраться, "как это работает", в частности, как представляются и обрабатываются данные.

В действительности, хотя обычно трудно предсказать время вычисления и размер результатов, знание принципов работы может дать представление о порядке их величины и при необходимости оптимизировать их. Эти оценки в действительности существенны: для большинства алгебраических вычислений результаты получаются почти моментально, и все идет отлично. Но если это не так, то требуемое время и память возрастают обычно экспоненциально. Таким образом, выполнимость данных вычислений не всегда очевидна, и глупо жертвовать значительными ресурсами, когда неудачу можно предсказать заранее. Например, если требуется найти собственные значения матрицы, то для программы на Фортране нет принципиальной разницы 100x100 или 500x500 эта матрица, так как время выполнения растет практически линейно. В Maple или , Mathematica вычисления с матрицей 5x5 могут занимать 15 секунд, в то время как те же по сути вычисления с матрицей 6x6 займут 15 минут.

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

Обычно, начиная работать с любой из систем символьных вычислений, пользователь достаточно легко решает небольшие и несложные примеры и задачи из учебника. Однако, приступая к решению настоящих проблем, пользователь сталкивается с рядом проблем: то компьютер слишком долго считает, то не хватает памяти, то в ответе получается формула на 5-10 страниц, а то машина выдает и вообще неправильный ответ. После этого встает вопрос - "Стоит ли тратить время на детальное изучение таких "игрушечных" систем и не лучше ли потратить это время на написание самих формул?".

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

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

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

Рассмотрим выражение x^100-1. Вид этого выражения "прост" и для компьютера и для пользователя. Разложим этот полином на множители:
(x-1)*(1+x^4+x^3+x^2+x)*(1+x^20+x^15+x^10+x^5)*(1+x)*(1-x+x^2-x^3+x^4)*
(1-
x^5+x^10-x^15+x^20)*(1+x^2)*(x^8-x^6+x^4-x^2+1)*(x^40-x^30+x^20-x^10+1).
Это выражение так же "просто" для компьютера, как и прежнее. Но если вы получите это выражение в результате своих вычислений, то вряд ли вид этого выражения удовлетворит вас. Данный пример показывает, что, даже получив правильный ответ с помощью компьютера, вы должны еще уметь представить его в удобном для пользователя виде.

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

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

2. Практическое занятие


Программы в Maple состоят из набора простых операторов (предложений), выполняемых последовательно. Предложения могут состоять из чисел, переменных, операторов, строк (ограниченных символами ``) и знаков арифметических и других операций. Целые числа имеют "бесконечную" точность, а числа, не являющиеся целыми, представляются в виде отношения двух целых чисел.Каждая переменная характеризуется типом и именем - набором символов, в которых большие и малые буквы различаются. Как обычно, имя не должно совпадать с существующими уже именами.

Данные синтаксические структуры являются общими для большинства языков высокого уровня.

В Maple выполняемые математические выражения вводятся всегда после символа >, а заканчиваются точкой с запятой или двоеточием, если результат не надо выводить на экран. Чтобы продолжить запись предложения на следующей строке используют комбинацию "Shift+Enter". При нажатии клавиши "Enter" предложение выполняется. Обнаружив ошибку, Maple выводит сообщение о ней в следующей строке.

Кроме обычных знаков математических операций используют :



** или ^ - возведение в степень,

! - факториал,

:= - знак присвоения,

<, >, >= , <= , = - логические операции,

% - результат предыдущей выполненной операции,

@ - сложная функция, т.е. ( f @ g )( x )= f(g(x)),

D - дифференциал и т.д.

Обратные кавычки указывают на строку символов, прямые кавычки (апостроф) используются для того, чтобы освободить переменную, которая уже использовалась. Заметим, что такие "парные" символы как кавычки, возведение в степень и некоторые другие вводятся с задержкой (советуем использовать пробел после введения такого "парного" символа).


Авторы всех систем символьных вычислений следуют принципу Кронекера - "Бог создал целые числа - все остальное дело рук человеческих". Во всех языках представление целых чисел первично, так как в этом случае алгоритмы сложения и умножения становятся, по крайней мере, в принципе, довольно простыми. В силу этого "ленивого" принципа Maple выдает только целые числа или имена без численных ответов до тех пор, пока вы об этом не попросите. Основные команды "принудительного" вычисления

evalb - для булевых выражений,
evalc - для комплексных выражений,
evalf - для операций с плавающей запятой,
evalm - для вычисления матричных выражений.

Лабораторная работа 1.

Разложение в ряд Тейлора


Рассмотрим разложение в ряд Тейлора. Начнем с очистки памяти системы и загрузки библиотеки графических подпрограмм:

> restart: with(plots):

Определим функцию от двух переменных

> g:=x->exp(-0.1*x^2)*sin(5*x);

точку, в окрестности которой мы изучаем разложение,

> x0:=0:

и радиус окрестности в которой мы будем изучать разложение функции в ряд

> r:=5:

Видно, что в заданной окрестности у функции нет особых точек

> plot(g(x),x=x0..x0+r,y=-1..1, color=red);

Видно, что построение графиков в Maple достаточно простое.
Теперь зададим начальную степень разложения в ряд Тейлора

> inDeg :=6:

и посмотрим на выражение для полинома Тейлора, используя встроенную функцию

> f:=taylor(g(x),x=x0,inDeg);

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

> N:= 10:

и на каждом шаге будем увеличивать степень полинома Тейлора на

> step:=6:

Построим теперь последовательность графиков. Для этого напишем цикл

> for i from 0 to N do


deg := inDeg + step*i:


f:=convert(taylor(g(x),x=x0,deg),polynom);

A[i]:= plot(f,x=x0..x0+r, y=-1..1, color=blue, title=cat("N=",deg)):

end do:

Написание цикла в Maple так же вполне очевидно.


Выведем теперь эти графики на экран в последовательности

> B:=display([seq(A[i], i=0..N)], insequence = true):

и добавим к ним изображение исходной функции

> G:= plot(g(x),x=x0..x0+r,y=-1..1, color=red);


display([G,B]);

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


Лабораторная работа 2

Сходимость ряда Тейлора


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

> g:=(x,y)->x*exp(-x^2-y^2):

точку, в окрестности которой мы изучаем разложение,

> x0:=2: y0:=1: zm:=0.1:

и окрестность, в которой мы будем изучать разложение в ряд

> r:=4:

Здесь zm – определяет максимальное значение функции g(x,y).
Убедимся, что в заданной окрестности у функции нет особых точек (нажмите мышкой на рисунок и, не отпуская, покрутите изображение)

> plot3d(g(x,y),x=x0-r..x0+r,y=y0-r..y0+r,grid=[25,15],style=patch,axes=boxed);

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

> inDeg := 2:

На каждом шаге будем увеличивать эту степень на

> step:=4:

Зададим общее число шагов

> N:= 8:

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

> for i from 0 to N do



deg := inDeg + step*i:

A[i] := plot3d(abs(mtaylor(g(x,y),[x=x0, y=y0], deg)-g(x,y)), x=x0-r..x0+r, y=y0-r..y0+r,

view= 0..0.2, orientation=[0,0],title=cat("Погрешность при N=",deg)):

end do:

Выведем графики на экран

> B:=display([seq(A[i], i=0..N)], insequence = true,

style=patch, view=0..0.2,axes=boxed):

добавив к ним изображение исходной точки

> C:=pointplot3d([x0,y0,0],color=red,symbol=circle,symbolsize=12, thickness=5):

Теперь можно посмотреть, как увеличивается радиус сходимости в зависимости от степени полинома

> display([B,C]);
Это проекция трехмерного графика, покрутив которую можно увидеть характерный "стаканчик", диаметр которого постепенно увеличивается с ростом степени полинома Тейлора.
Рассмотрим теперь функцию с сингулярностью, которая будет ограничивать увеличение радиуса сходимости.

> g:=(x,y)->1/(x^2+y^2);



G:=plot3d(g(x,y),x=x0-r..x0+r,y=y0-r..y0+r,view=0..60,grid=[25,15],style=patch,axes=boxed):

display([G,C]);

Красным отмечена точка, в окрестности которой рассматривается разложение ряд Тейлора.


Построим графики погрешности аппроксимации

> for i from 0 to N do



deg := inDeg + step*i:

A[i] := plot3d(abs(mtaylor(g(x,y),[x=x0, y=y0], deg)-g(x,y)), x=x0-r..x0+r, y=y0-r..y0+r,

view= 0..0.2, orientation=[0,0],title=cat("Погрешность при N=",deg)):

end do:

Выведем графики на экран

> B:=display([seq(A[i], i=0..N)], insequence = true,

style=patch, view=0..0.2,axes=boxed):

добавив к ним изображение исходной точки и сингулярности, которая ограничивает область сходимости ряда

> F:=pointplot3d([[0,0,0]], color=blue,symbol=circle,symbolsize=12, thickness=5):

display([B,C,F]);

Наглядно видно, что область сходимости не увеличивается, в отличие от предыдущего случая!



Лабораторная работа 3

Вычисление интегралов


Используя цилиндрические координаты найти объем тела, ограниченного поверхностью и плоскостями z=0, z=2.

В цилиндрических координатах поверхность определяется уравнением



. Нарисуем эту поверхность

> restart:



x:=cos(2*theta)^2*cos(theta):

y:=cos(2*theta)^2*sin(theta):

surf:=plot3d([x,y,z],theta=0..2*Pi,z=-1..3,grid=[60,5]):

и две плоскости

> planes:=plot3d({0,2},x=-1.5..1.5,y=-1.5..1.5):

выведем графики на экран

> with(plots): display3d({surf,planes});

Для того, чтобы определить область интегрирования рассмотрим проекцию на плоскость XY

> polarplot((cos(theta)^2-sin(theta)^2)^2,theta=0..2*Pi);
Теперь очевидно, что область интегрирования имеет вид

0 < z < 2 , 0 < θ < 2π , 0 < r < (cos^2 (θ) -- sin^2 (θ) ^ 2

Таким образом,Объем тела равен

> v:=Int(Int(Int(rho,rho=0..(cos(theta)^2-sin(theta)^2)^2),

theta=0..2*Pi),z=0..2);

Вычислим этот интеграл

> value(v);


3. Практическое занятие

Лабораторная работа 4

Метод наименьших квадратов

Напомним, как найти аппроксимирующую функцию в классе полиномов, т.е. когда требуется функцию g, заданную в n точках , приблизить многочленом степени :



Итак, нам необходимо найти полином, который обеспечивает минимум суммы



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



= k=0,..m

Таким образом, для определения коэффициентов необходимо решить систему (m+1) линейных уравнений, которые называют нормальными уравнениями, с (m+1) неизвестными . B матричной форме эту систему уравнений можно записать в виде с=Ba, где B - симметрическая матрица размерности (m+1).


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

Обычно, при фиксированном напряжении U изменяется сопротивление R и измеряется сила тока J.


Загрузим необходимые библиотеки

> restart:



with(plots): with(stats):

и смоделируем экспериментальную установку, используя известный нам по теории закон Ома,

> f:=r->U/r;

Для этого зададим напряжение

> U:=100:

и число измерений

> N:=10:

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

> for i from 1 to N do

R[i]:=i; J[i]:=f(R[i]);

end do:

Таким образом, мы построили график “идеальных” измерений

> Data:=[ [R[k],J[k]] $k=1..N]:

plot([100/r,Data],r=R[1]..R[N],style=[line,point],
symbol =circle,color=[red,maroon]);

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

> NoiseR:=stats[random, normald](N);

> NoiseJ:=stats[random, normald](N):

> for i from 1 to N do



Rn[i]:=R[i]*(1+0.03*[NoiseR][i]/max(NoiseR));

Jn[i]:=J[i]*(1+0.04*[NoiseJ][i]/max(NoiseJ));

end do:

Здесь сопротивление измеряется с погрешностью, не превосходящей 3%, а ток - 4% процентов. Эта погрешность указывается на реальных физических приборах.


Таким образом, после измерения мы получим "экспериментальные" данные

> DataN:=[ [Rn[k],Jn[k]] $k=1..N]:



A:=plot(DataN,style=point,symbol = circle,color=maroon):

%;
По этим “экспериментальным”данным нам необходимо найти напряжение U. Для обработки данных по методу наименьших квадратов загрузим библиотеку

> with(CurveFitting):

и воспользуемся методом наименьших квадратов, используя команду LeastSquares, в которой мы явно указываем предполагаемый нами вид зависимости J=U/R

> g:=LeastSquares(DataN, r, curve=a/r);



Итак, что искомое напряжение примерно равно U=103. Проверьте, что увеличение числа измерений N приводит к улучшению результата.


Сравним график теоретической зависимости и экспериментальных данных

> B:=plot(g(r),r=Rn[1]..Rn[N],color=red):



display([A,B]);

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

> h:=LeastSquares(DataN, r, curve=a*exp(-0.6*r)+c);

> B:=plot(h(r),r=Rn[1]..Rn[N],color=red):



display([A,B]);

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

> describe[linearcorrelation]([Rn[k] $k=1..N],[Jn[j] $j=1..N]): evalf(%);

-0.7943085620

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



Так как погрешность мы вводим генератором случайных чисел, то вы можете получить другой ответ.



Лабораторная работа 5

Быстрое преобразование Фурье


Рассмотрим функцию f(t) от одной переменной или сигнал, который зависит от времени.
Загрузим необходимые нам библиотеки процедур

> restart: with(linalg): with(plots):

и определим интервал [0,T] на котором мы собираемся изучать поведение этой функции

> T:=0.16;

Определим стационарный сигнал - сумму двух периодических функций, частоты которых

> w1:=50: w2:=200:

не меняются со временем

> f:=t-> sin(2*Pi*w1*t)+sin(2*Pi*w2*t);

p1:=plot(f(t), t=0..T, title=`Сигнал 1`):

%;

Используем метод дискретизации, при котором сигнал f(t) заменяется совокупностью его мгновенных значений, называемых выборками, или отсчетами (samples).


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

> dT:=evalf(Pi/max(w1,w2));

Здесь max(w1,w2) - максимальное значение частоты в спектре сигнала.
Далее мы добавим к этоиу сигналу высокочастотный "белый" шум и, поэтому, заранее выберем шаг дискретизации меньше порогового значения

> dT:=0.001;

Зададим также величину выборки - т.е. число отсчетов N

> m:=10: N:=2^m;

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

> fk:= array(1..N): x:= array(1..N):

и заполним их

> for i to N do

fk[i]:=evalf(f(i*dT)); x[i]:=i*dT;

end do:

Добавим теперь к этим значениям произвольный шум

> r:=randvector(N): fk:=evalm(fk+.02*r):

и сравним полученный нами "реальный" сигнал с "идеальным" сигналом

> maxN:=T/dT: Signal:=[ [x[j],fk[j]] $j=1..maxN]:

> pS:=plot(Signal,color=COLOR(RGB, 0.196, 0.6, 0.8),


title="
Сигнал + Шум"):
display([pS,p1]);

Естественно возникают два вопроса - какую информацию об исходной функции можно получить из "реального" сигнала и как очистить этот сигнал от шума?


Для решения этой задачи применяют быстрое преобразование Фурье

> im:=vector(N, 0):

> FFT(m,fk,im);

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


Вычислим максимальную частоту

> dw:=2*Pi/(N*dT);

отвечающую используемой нами выборке и определим ось частот

> fw := x -> (x-1)*evalf(dw/(2*Pi)):



w := vector(N, fw):

для того, чтобы построить спектр сигнала

> Data:=[ [w[j],abs(fk[j])] $j=1..N]:

> p2:=plot(Data,color=coral, title="Спектр"):



%;

Ярко выраженные максимумы при и указывают на наличие в сигнале двух основных частот.


Используя эту информацию построим фильтр - маску, шириной

> dW:=9:

Фильтр - это функция

> filtr:=k->piecewise(k>w1-dW and kN-w1-dW and kw2-dW and kN-w2-dW and k

которая выглядит следующим образом

> plot(filtr,0..N,color=cyan, thickness=2);

Наложим этот фильтр на спектр сигнала

> for i from 1 to N do



fk1[i]:=filtr(i)*fk[i]:

im1[i]:=filtr(i)*im[i]:

end do:

и посмотрим на полученный результат

> plot([[w[j],abs(fk1[j])] $j=1..N],color=coral, title="Спектр+Фильтр");

Теперь восстановим сигнал с помощью обратного преобразования Фурье

> iFFT(m,fk1,im1);

и сравним очищенный от шумов результат с исходной непрерывной функцией

> SF:=[ [x[j],fk1[j]] $j=1..maxN]:

pF:=plot(SF,color=COLOR(RGB, 0.2, 0.5, 0.4),
title="
Очищенный сигнал"):

display([pS,p1]);

display([pF,p1]);

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


4 Практическое занятие

Задачи:


1. Построить разложение в ряд Фурье функции

> g:=x->Heaviside(x+4)+2*Heaviside(x+1)-2*Heaviside(x-1)-Heaviside(x-4):

на интервале [-4.2..4.2]. Коэффициенты в разложении определить по формулам Эйлера-Фурье.

2. Построить разложение в ряд Фурье функции


> f:=x->piecewise(x<1/2,x, x<=1, (x-1),0);

на интервале [0..1]. Коэффициенты в разложении определить по формулам Эйлера-Фурье.

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

4. Найти объем тела ограниченного поверхностями - сверху, и - снизу, над кругом единичного радиуса

5. Найти среднее значение функции над областью, ограниченной кривыми , ,

6. Используя встроенную команду interp построить интерполяционный полином для функции sin(x) на интервале [0,π] и изучить зависимость качества интерполяции в зависимости от степени полинома.

7. Построить интерполяционный полином в форме Ньютона для функции sin(x) на интервале [0,π] и изучить зависимость качества интерполяции в зависимости от степени полинома.

8. Используя быстрое преобразование Фурье изучить возможность очистки от 3% белого шума сигнала
> nu:=280:

F:=x->5*cos(2*nu*Pi*x)*(exp(-640*Pi*(x-1/8)^2)+

exp(-640*Pi*(x-3/8)^2)+

exp(-640*Pi*(x-4/8)^2)+

exp(-640*Pi*(x-6/8)^2)+

exp(-640*Pi*(x-7/8)^2));

Эту функцию можно рассматривать как аналоговый сигнал, отвечающий цифровому сигналу вида [1011011].


Литература


Дэвенпорт Дж, Сирэ И., Турнье Э. Компьютерная алгебра, Мир, 1991.

Цыганов А.В. Символьные вычисления, Методическое пособие, СПбГУ, 1999, www.exponenta.ru



Цыганов А.В. Квантовая механика с Maple, РХД, 1999, www.exponenta.ru





Смотрите также:
Системы символьных вычислений. Теоретическое занятие
260.22kb.
1 стр.
Балтабаева шолпан аипбергеновна
446.46kb.
4 стр.
Лекция История развития вычислений и структур вычислительных систем Сазанов В. М
549.72kb.
4 стр.
Международное гуманитарное право. Вводное занятие
61.25kb.
1 стр.
Современные тенденции в разработке средств визуализации программного обеспечения параллельных вычислений
34.32kb.
1 стр.
Профессор Петров Ю. П. Обеспечение достоверности и надежности компьютерных вычислений
1118.04kb.
5 стр.
Интерактивные системы доказательства
140.92kb.
1 стр.
Конкурс проектов по применению высокопроизводительных вычислений информационные технологии
28.44kb.
1 стр.
Стиральная машина Eurosoba 1000 Содержание. Установка
203.14kb.
1 стр.
Red Hat и Cisco предложили ос enterprise Linux и виртуализацию для новой унифицированной среды вычислений
44.87kb.
1 стр.
Занятие «Вода путешественница»
65.25kb.
1 стр.
Занятие в учреждении дополнительного образования детей
187.67kb.
1 стр.