Аппроксимации рядом Тейлора

Начнем с аппроксимации функции хорошо известным рядом Тейлора степени 8 относительно середины интервала (точки с х=2):

Аппроксимация Чебышева-Паде

Теперь рассмотрим еще более точную рациональную аппроксимацию Чебышева-Паде. Это такая рациональная функция r[m, n](х) с числителем степени т и знаменателем степени п такой же, как и для разложения в ряд Чебышева. Функция r [m, n](х) согласуется с разложением в ряд Чебышева f(x) членом степени m+n. Мы вычислим аппроксимацию Чебышева-Паде степени (4,4), подобную обычной Паде- аппроксимации, успешно выполненной ранее:

Аппроксимация полиномами Чебышева

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

Разложим функцию f(x) на [0, 4] в ряд Чебышева с точностью 1*10-8. Это означает, что все члены с коэффициентами меньше чем эта величина, будут опущены. Такая точность обеспечивается полиномом 13 степени:

Движение частицы в магнитном поле

От реального мира перейдем к микромиру. Пусть микрочастица массой 9* 10-31 кг и зарядом +1,6*10"19 Кл влетает в магнитное поле с индукцией В = 0,1 Тл под углом а=80°. Рассчитаем траекторию движения частицы при начальной скорости Vo= 1*107м/с:

> restart;

Сила Лоренца, действующая на движущуюся частицу F = q*(E+[v, В]). Проекции векторного произведения [v, В] на оси х, у, z:

[v.B]x=vy*Bz-vz*By [v,B]y=vz*Bx-vx*Bz [v,B]z=vx*By-vy*Bz

В соответствии с этим известные из курса физики дифференциальные уравнения, описывающие траекторию полета частицы по осям х, у, z имеют вид:

Эффективная оценка рациональных функций

Полиномы числителя и знаменателя в минимаксной аппроксимации уже выражены в форме Горнера (то есть в форме вложенного умножения). Оценка полиномом степени п в форме Горнера при n-умножениях и n-суммированиях — это наиболее эффективная схема оценки для полинома в общей форме. Однако для рациональной функции степени (т, п) мы можем делать кое-что даже лучше, чем просто представить выражения числителя и знаменателя в форме Горнера. Мы можем нормализовать рациональную функцию так, что полином знаменателя будет со старшим коэффициентом, равным 1. Мы можем также заметить, что вычисление рациональной функции степени (т, п) в форме Горнера требует выполнения все m+n сложений , m+n-1 умножений и 1 деления. Другими словами, общий индекс действия есть:

m+n операций умножения/деления; m+n операций сложения/вычитания.

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

max(m,n) операций умножения/деления; m+n операций сложения/вычитания.

Например, если m = n, тогда эта новая схема требует выполнения только поло-, вины числа действий умножения/деления по сравнению с предшествующим методом. Для рациональной функции MlnimaxApprox вычисление в форме, выраженной выше, сводится к 9 действиям умножения/деления и 8 действиям сложения/вычитания. Число операций умножения/деления можно сократить до 8, нормализуя знаменатель к форме monic. Мы можем теперь вычислить непрерывную (цепную) дробь для той же самой рациональной функции. Вычисление по этой схеме, как это можно видеть из вывода Maple, сводятся только к 4 действиям деления и 8 действиям сложения/вычитания:

> MinimaxApprox := confracform(MinimaxApprox):

> lprint(MinimaxApprox(x));

-.468857770747е-1+1.07858705749/(х+4.41994843227+16.1901737091/ (х+4.29121842830+70.1948525272/(х-10.2912843004+ 4.77536150167/(х+1.23883665458))))

Малосигнальный анализ усилителя на полевом транзисторе

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

Минимаксная аппроксимация

Классический результат теории аппроксимации заключается в том, что минимакс как наилучшая аппроксимация рациональной функции степени (т, п) достигается, когда кривая ошибки имеет m+n+2 равных по величине колебаний. Кривая ошибки аппроксимации Чебышева-Паде имеет нужное число колебаний, но эта кривая должна быть выровнена (по амплитуде выбросов кривой ошибки) с тем, чтобы обеспечить наилучшее минимаксное приближение. Эта задача решается с помощью функции minimax:

Моделирование цепи на туннельном диоде

А теперь займемся моделированием явно нелинейной цепи. Выполним его для цепи, которая состоит из последовательно включенных источника напряжения Es, резистора Rs, индуктивности L и туннельного диода, имеющего N-образную вольтамперную характеристику (ВАХ). Туннельный диод обладает емкостью С, что имитируется конденсатором С, подключенным параллельно туннельному диоду. Пусть ВАХ реального туннельного диода задана выражением:

> restart:

> A:=.3t: а:=10: В:=1*10^(-8): b:=20:

> Id:=Ud->A*Ud*exp(-a*Ud)+B*(exp(b*Ud-D):

Id:=Ud->AUde(-aUd)+Be(bUd-1)

Построим график ВАХ:

> plot(Id(Ud), Ud=-.02..0.76,color=black):

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

Работа цепи описывается системой из двух дифференциальных уравнений:

di/dt=(Es-i(t)*Rs-u(t))/L

du/dt=(i(t)-Id(u(t))/C

Моделирование физических явлений

Расчет траектории камня с учетом сопротивления воздуха

Вы хотите метнуть камень в огород вашего вредного соседа? Разумеется, во время его отсутствия. Давайте промоделируем эту ситуацию, предположив два актуальных случая: дело происходит на Луне и на Земле. В первом случае сопротивления воздуха (как и его самого) нет, а в другом — сопротивление воздуха есть и его надо учитывать. Иначе камень упадет в ваш огород, а не в огород соседа!

Итак, пусть подвернувшиеся под руку камни с массой 500 и 100 г брошены под углом 45° к горизонту со скоростью Vo = 20 м/с. Найдем их баллистические траектории, если сила сопротивления воздуха Fтр=А*V, где А=0,1 Н*с/м. Сравним их с траекториями, получающимися без учета сопротивления воздуха.

Начнем с подключения пакета plots, нужного для визуализации данной задачи:

> restart;

> with(plots):

Warning, the name changecoords has been redefined

Составим параметрические уравнения для проекций скорости на оси координат:

> Vox:=Vo*cos(a1pha):Voy:=Vo*sin(alpha):

Vox:= Focos(a)

Voy :=Vo sin(a)

Мы рассматриваем два случая: камень массой 500 г и камень массой 100 г. Поскольку для каждого случая мы предусматриваем расчет в двух вариантах (с учетом сопротивления воздуха и без такого учета), то мы должны составить 4 системы дифференциальных уравнений (ДУ). Каждая система состоит из двух ДУ второго порядка и вид этих систем известен из курса физики. Ниже представлено задание этих систем ДУ (для первой системы дан вывод ее вида):

Моделирование и расчет электронных схем

Нужно ли применять Maple для моделирования и расчета электронных схем?

Нужно ли применять системы компьютерной математики для анализа, расчета и моделирования электронных схем? Ответ на этот вопрос не так прост, как кажется с первого взгляда. С одной стороны, к услугам пользователя компьютера сейчас имеется ряд программ схемотехнического моделирования, например Micro-CAP, Electronics Workbench, PSpice, Design Labs и др., автоматически составляющих и решающих большие системы уравнений состояния электронных схем и моделирующих работу бесчисленного множества электронных схем без кропотливого «ручного» составления уравнений.

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

Моделирование рассеивания альфа- частиц

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

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

Итак, пусть в нашем теоретическом опыте альфа- частицы с энергией 4 МэВ рассеиваются тонкой золотой фольгой. Рассчитать траекторию частицы, приближающейся к ядру атома Аи. Прицельное расстояние р равно 2*10-15 м. Приступим к решению задачи и зададим вначале систему дифференциальных уравнений для траектории альфа- частицы:

Небольшое введение

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

Многие читатели полагают, что системы компьютерной математики хорошо работают на таких простых примерах, но от них мало толку при решении реальных задач математики, физики или радиоэлектроники. Это, конечно, заблуждение. Дело просто в том, что при решении таких задач руководящая роль пользователя сильно возрастает. Вы должны понимать, что не Maple 7 решает вашу задачу, а вы! И система Maple 7 лишь помогает в этом трудном деле. Так что при неудачах в решении своих специфических задач следует прежде всего пенять на себя и на свое незнание возможностей системы Maple 7, а вовсе не на свою помощницу.

В том, что Maple можно успешно использовать при решении вполне конкретных научных и практических задач, призваны убедить примеры, приведенные ниже. Разумеется, и их нельзя отнести к таким сложнейшим задачам, как проектирование ядерного реактора или расчет траектории полета космического корабля, — не стоит забывать, что такие расчеты делают на суперкомпьютерах, а не на домашнем компьютере, который стоит перед вами. И объем материалов по сопровождению и результатам таких расчетов многократно превосходит объем всей этой книги. Тем не менее в этом уроке вы встретите решение вполне реальных и полезных задач в области математики, физики и радиоэлектроники. Почему не в механике, гидродинамике или в оптике? Да потому, как верно сказал наш народный пророк Козьма Прутков: «нельзя объять необъятное». Приведенные примеры отчасти обусловлены личными пристрастиями автора, но они полезны каждому пользователю, желающему всерьез оценить возможности Maple 7.

Описанные в этом уроке задачи являются реальными документами, созданными и отлаженными в среде Maple 7 и лишь затем перенесенными в рукопись книги. Так что они заодно служат примерами того, как надо оформлять такие документы. В то же время от некоторых «излишеств» оформления (например, закрывающихся и открывающихся секций) мы отказались, дабы не усложнять описание документов явно второстепенными деталями. Начнем этот урок с решения весьма актуальной для многих областей применения математики задачи — аппроксимации сложной функции.

Паде-аппроксимация

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

Преобразование в код Фортрана или С

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

Применение интеграла Дюамеля для расчета переходных процессов

Вернемся к линейным цепям и рассмотрим еще один полезный метод расчета электрических цепей — с помощью интеграла Дюамеля. При нем можно рассчитать временную зависимость выходного напряжения u2(t) цепи по известному входному сигналу ul(t) и переходной характеристике цепи a(t). Возьмем в качестве первого классического примера дифференцирующую RC-цепь и вычислим ее реакцию на экспоненциально нарастающий перепад напряжения.

Представлены заданные зависимости ul(t) и a(t), аналитическое выражение для интеграла Дюамеля (одна из 4 форм) и аналитическое выражение для искомой зависимости u2(t). Пока последнее выражение довольно простое. В конце этого фрагмента документа построены графики зависимостей ul(t), a(t) и u2(t).

Обратите внимание на то, что выражение для u2(t), получаемое с помощью интеграла Дюамеля, стало намного сложнее. Тем не менее получено как аналитическое выражение для реакции цепи u2(t), так и графики ul(t), a(t) и u2(t). Они показаны внизу графика.

Проектирование цифрового фильтра

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

Спроектируем фильтр N+1-ro порядка класса FIR (Finite Impulse Response или с конечной импульсной характеристикой). Каждая из N ячеек временной задержки фильтра удовлетворяет следующей зависимости выходного сигнала у от входного х вида:

Расчет аналогового фильтра на операционном усилителе

Теперь рассмотрим проектирование аналогового полосового фильтра на операционном усилителе, схема которого приведена на Рисунок 17.16.

Разделение изотопов

Рассмотрим еще одну классическую задачу ядерной физики — разделение изотопов (атомов с одинаковым зарядом ядра, но разной массой). Для этого используют различные способы. В частности, это может быть масс-спектроскопический метод. Из точки А вылетают однозарядные ионы (q = е = 1.6*10-19 Кл) разной массы (от 20 до 23 а.е.м.) и под разными углами в пределах от 80 до 100° к оси х в плоскости ху (Рисунок 17.9). Вдоль оси z приложено магнитное поле В=10-2 Тл. Рассчитаем траектории полета частиц. Будем надеяться, что это подскажет способ разделения изотопов.

Приступим к решению данной задачи. Сила Лоренца, действующая на движущуюся частицу, F = q*(E+[v, В]). Проекции векторного произведения [v, В] на оси х, у, z заданы выражениями:

[v.B]x-vy*Bz-vz*By [v,B]y-vz*Bx-vx*Bz [v,B]z=vx*By-vy*Bz

нами показано, что правильный выбор

Шаг 1

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

Шаг 1

Зададим исходные числовые данные (опустив размерности):
> q:=-1.6e-19: massa:=9.1e-31: V:=le7: alpha:=80*Pi/180:
> Vx:=V*cos(alpha): Vy:=V*sin(alpha): Ex:=0:Ey:=0:Ez:=0: Bx:=0.1:By:=0: Bz:=0:
Построим траекторию движения частиц в пространстве:
> with(DEtools):DEplot3d({sys},{x(t),y(t),z(t)},t=0..2e-9, [[x(0)=O,D(x)(0)=Vx,y(0)==0,D(y)(0)=Vy,z(0)=0,D(z)(0)=0]], stepsize=le-ll,orientation=[24.117]):
Полученная траектория представлена на Рисунок 17.8. Она имеет вид спирали в пространстве. При этом скорость движения частицы вдоль оси х неизменна, а вдоль осей у и z имеет характерную колебательную компоненту. Случай явно куда менее тривиальный, чем полет камня, описанный выше.

Шаг 1

Зададим исходные числовые безразмерные данные для расчета:

Шаг 1

Рисунок 17.9. Иллюстрация к методу разделения изотопов
В соответствии с этим дифференциальные уравнения, описывающие траекторию полета частицы по осям х, у, z имеют вид:

Шаг 1

Введем исходные числовые данные для вычислений:
> ql:=2*i;6e-19:q2:=79*1.6e-19:massa:=4*1.67e-27:EO:=8.85e-12: a:=4e-13:
p:=5e-15:T:=4e6*1.6e-19:V0x:=sqrt(2*T/massa):
Создадим графическую структуру решения нашей системы дифференциальных уравнений для нескольких расчетных отклонений линии движения альфа- частицы от центра ядра атома, находящегося на ее пути:
> with(DEtools):ss:=DEplot({sys},{y(t),x(t)},t=0..7e-20.
[[x(0)=-a,D(x)(0)=VOx,y(0)=p,D(y)(0)=0].
[x(0)=-a,D(x)(0)=VOx,y(0)=p*4.D(y)(0)=0],
[x(0)=-a,D(x)(0)=VOx,y(0)=p*8,D(y)(0)=0],
[x(0)=-a,D(x)(0)=VOx,y(0)=p*12,D(y)(0)=0].
[x(0)=-a;D(x)(0)=VOx,y(0)=p*16,D(y)(0)=0],
[x(0)-a.D(x)(0)-VOx.y(0)-p*20,D(y)(0)-0].
[x(0)=-a,D(x)(0)=VOx,y(0)=p*24,D(y)(0)=0],
[x(0)=-a,D(x)(0)=VOx,y(0)=p*28,D(y)(0)=0]],
x(t)=-a..a,scene=[x(t),y(t)],stepsize=le-21,1inecolor=bl ack):
> with(plottools):yy:=circle([0.0],2E-14,color=red,thickness=2):
Warning, the name translate has been redefined
Построим центр ядра (кружок со знаком +) и траектории альфа- частиц:
> ss2:=PLOT(TEXT([0.-0.3e-14],'+'), FONT(HELVETICA, OBLIQUE.14)):
Осталось построить график траекторий движения альфа- частиц вблизи центра атома: i
> with(plots):
Warning, the name changecoords has been redefined
> disp1ay([ss,yy,ss2],tit1e='Pacceивание а-частиц',axes=framed);
График траекторий движения альфа- частиц вблизи ядра представлен на Рисунок 17.11. Этот график настолько нагляден, что не требует пояснения.
Моделирование движения альфа- частиц вблизи малого и «массивного» ядра атома дает наглядное представление о математической и физической сути данного опыта. Надо лишь помнить, что нельзя нацеливать альфа- частицы прямо в центр ядра. Более сложные, чем приведенные, расчеты показывают, что при этом альфа-частица настолько близко подходит к ядру, что надо учитывать новые факторы, возникающие при близком взаимодействии. Они могут привести к тому, что частица будет поглощена ядром- Но это уже тема нового разговора,, выходящего за рамки данной книги.

Шаг 1

Рисунок 17.16. Схема полосового фильтра на интегральном операционном усилителе
Подготовимся к расчету фильтра:
> restart:
Зададим основные уравнения, описывающие работу фильтра на малом сигнале:

Мы можем найти аналитическое представление

Шаг 2

Рисунок 17.8. Траектория движения частицы в магнитном поле
Мы можем найти аналитическое представление для траектории частицы в виде параметрически заданной (с параметром времени t) системы из трех уравнений:

Шаг 2

Выполним решение заданных систем ДУ:

Шаг 2

Зададим исходные числовые данные для расчета:
> q:=1.6e-19:V:=le4:
> Vx:=V*cos(a1pha):Vy:=V*sin(a1pha):Ex:=0:Ey:=0:Ez:=0:Bx:=0: By:=0:Bz:=le-2:
Выполним решение составленной выше системы дифференциальных уравнений:

Шаг 2

Рисунок 17.11. Траектории движения альфа- частиц вблизи ядра атома

Шаг 2

Введем круговую частоту:
> omega := 2*Pi*f;
W := 2пf
Найдем коэффициент передачи фильтра и его фазо- частотную характеристику как функции от частоты:
> gain := abs(eva1c(Vo/Vi)):
> phase := evalc(op(2,convert(Vo/Vi.polar))):
Для просмотра громоздких аналитических выражений для этих параметров замените знаки двоеточия у выражений для gain и phase на знак точки с запятой. Далее введем конкретные исходные данные для расчета:
> R3 :=1000:
> R4 := 3000:
> СЗ :=0.08*10^(-6):
> С4 := 0.01*10^(-6):
Построим АЧХ фильтра как зависимость коэффициента передачи в децибелах (dB) от частоты f в Гц:
> plot(DogWf), 20*log10(gain), f=[10..50000], color=black, title='Коэффициент передачи dB как функция от частоты f в Гц'):
Эта характеристика представлена на Рисунок 17.17. Здесь полезно обратить внимание на то, что спад усиления на низких и высоких частотах происходит довольно медленно из-за малого порядка фильтра.

Шаг 2

б
Рисунок 17.12. Принципиальная (о) и эквивалентная(6) схемы усилителя на полевом транзисторе
Наша цель заключается в расчете характеристик усилителя операторным методом. Подключим нужный нам пакет plots:
> restart:with(plots):
Warning, the name changecoords has been redefined
Из законов Киргофа вытекает, что сумма токов, втекающих в каждый узел и вытекающих из него равна 0. Следовательно, для узлов эквивалентной схемы Рисунок 17.12 можно записать следующую систему уравнений в операторной форме:

Шаг 2

Вычислим FIR- коэффициенты для прямоугольного окна фильтра:
> С :-=(n) -> limit(g,t=n):h := aray(0..N): N2:=N/2:
> for n from 0 to N2 do h[N2-n]:= evalf(C(n)); h[N2+n] := h[N2-n]; od:
Определим массивы входного x(n) и выходного у(n) сигналов:
> х := array(-N..T): y := аггау(0..Т):
Установим значение х(n) равным 0 для времени меньше 0 и 1 для времени >=0:
> for n from -N to -1 do x[n] := 0; od:
> for n from 0 to Т do x[n] := Dirac(n); od:
Вычислим временную зависимость для выходного сигнала:
> for n from 0 to Т do y[n] := sum(h[k]*x[n-k],k=0..N): od:
Построим график импульсной характеристики фильтра, отражающей его реакцию на сигнал единичной площади с бесконечно малым временем действия:
> р := [seq([j/fs,y[j]],j=0..T)3:
> plot(p, time=0..3*N/fs, labels=[time,output], axes=boxed, xtickmarks=4, title-'Иипульсная характеристика фильтра',color=black);
Он показан на Рисунок 17.19. Нетрудно заметить, что эта характеристика свидетельствует об узкополосности фильтра, поскольку его частоты fl и fh различаются несильно. В этом случае полосовой фильтр по своим свойствам приближается к резонансному, хотя само по себе явление резонанса не используется.

Шаг 2

Поскольку заведомо известно, что схема имеет малые значения L и С, мы задали с помощью параметров достаточно малый шаг решения для функции dsolve — stepsize=l(T(-11) (с). При больших шагах возможна численная неустойчивость решения, искажающая форму колебаний, получаемую при моделировании. Используя функции odeplot и displ ay пакета plots, построим графики решения в виде временных зависимостей u(t) и 10*i (t) и линии, соответствующей напряжению Es источника питания:
> gu:=odeplot(F,[t,u(t)],0,tm,color=black,
labels=['tVu(t),10*i(tr]):
> gi:=odeplot(F,[t,10*i(t)],0..tm.color-black):
> ge:=odeplot(F,[t,Es].0..tm.color=red): .
> display(gu.gi,ge);
Эти зависимости представлены на Рисунок 17.26. Из них хорошо видно, что цепь создает автоколебания релаксационного типа. Их форма сильно отличается от синусоидальной.

Шаг 2

Рисунок 17.1. График аппроксимируемой функции
Итак, вычисление f(x) по ее интегральному представлению совершенно не эффективно. Наша цель состоит в разработке процедуры вычислений, которая дала бы 6 точных цифр результата в интервале [0..4] и требовала, по возможности, наименьшего числа арифметических операций для каждого вычисления. Втайне не вредно помечтать о том, чтобы после аппроксимации время вычислений уменьшилось бы хотя в несколько раз. Что получится на деле, вы увидите чуть позже. А пока войдем в дебри аппроксимации.

Шаг 2

Кривая ошибок для аппроксимации полиномом Тейлора строится командой:
> plotd(f- TaylorApprox,0..4,.co1or=black);
и имеет вид, представленный на Рисунок 17.2. Эта кривая нас, прямо скажем, не слишком радует, поскольку погрешность в сотни раз превышает заданную.

Шаг 2

и имеет вид, показанный на Рисунок 17.3.

Шаг 2

Схема Горнера минимизирует число арифметических операций, заменяя операции возведения в степень операциями последовательного умножения.

Шаг 2

Рисунок 17.4. Кривая ошибки при Паде-Чебышева рациональной аппроксимации
> maxChebPadeError :=abs( F(0) - ChebPadeApprox(O) );
maxChebPadeError= .1236746 10-5
Мы достигли впечатляющего успеха и остается сделать еще один шаг в направлении повышения точности аппроксимации.

Шаг 2

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

с магнитным полем показывает, что

Шаг 3

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

Шаг 3

Рисунок 17.17. АЧХ фильтра на операционной усилителе
Далее построим фазо-частотную характеристику фильтра как зависимость фазы в радианах от частоты f в Гц:
> plot ([log10(f),phase, f=10..50000], color=black, title=*Фазо-частотная характеристика фильтра*);
Фазо-частотная характеристика (ФЧХ) фильтра показана на Рисунок 17.18
На ФЧХ фильтра можно заметить характерный разрыв, связанный с превышением фазовым углом граничного значения я. Такой способ представления фазового сдвига общепринят, поскольку его изменения стремятся вписать в диапазон от -я до п.

Шаг 3

Переменные напряжения на узлах схемы находятся из аналитического решения данной системы. При этом заблокируем вывод их аналитических значений, поскольку он очень громоздок. Тем не менее вы можете посмотреть на полученные формулы, поставив знак точки с запятой вместо знака двоеточия в приведенных ниже выражениях:
> solve({eql,eq2,eq3.eq4}б{Vl,V2.V3,Vo}):
Обеспечим присвоение переменным Vo, VI, V2 и V3 найденных из решения системы уравнений значений:
> assign(%):
Теперь найдем операторную передаточную функцию в аналитическом виде:

Шаг 3

Рисунок 17.19. Импульсная характеристика цифрового фильтра
Вычислим АЧХ фильтра, используя прямое преобразование Фурье. Оно после подготовки обрабатываемых массивов реализуется функцией FFT:
> rо := array (1..T+1): io := arrayd. .T+l):
> for n from 0 to Т do ro[n+l] :- y[n]; io[n+l] := 0; od:
> FFT(m,ro,io):
Построим график АЧХ фильтра:
> р :=[seq([j*fs/(T+l),abs(ro[j+l]+io[j+l]*I)3,j=O..T/2)]:
> plot(p, frequency=0..fs/2, tabels=[frequency,gain], tit1e='AЧX фильтра',со1ог=black);
Он представлен на Рисунок 17.20. Нетрудно заметить, что и впрямь АЧХ фильтра напоминает АЧХ резонансной цепи — она имеет вид узкого пика. Вы можете легко проверить, что раздвижением частот fl и fh можно получить АЧХ с довольно плоской вершиной и резкими спадами (говорят, что такая характеристика приближается к прямоугольной).

Шаг 3

Рисунок 17.26. Временные зависимости напряжения на туннельном диоде и тока
Решение можно представить также в виде фазового портрета, построенного на фоне построенных ВАХ и линии нагрузки резистора Rs:
> gv:=plot({Id(Ud),(Es-Ud)/Rs},Ud=-.05..0.75,color=black,
labels=[Ud,Id]):
> gpp:=odeplot(F.[u(t),i(t)],0..tm,color=blue):
> display(gv,gpp);
Фазовый портрет колебаний показан на Рисунок 17.27.

Шаг 3

Рисунок 17.2. Кривая погрешности при аппроксимации рядом Тейлора
Типичное свойство аппроксимации рядом Тейлора состоит в том, что ошибка мала вблизи точки разложения и велика вдали от нее. В данном случае самая большая ошибка имеет место в левой оконечной точке. Чтобы вычислить значение ошибки в точке х =0, что ведет к делению на нуль (см. определение для f(x)), мы должны использовать значение предела:
> maxTaylorError := abs( Limit(f(x), х-0) - ТауlorАрргох(0) );
maxTaylorError := .0015029620
Итак, в самом начале наших попыток мы потерпели полное фиаско. Но отчаиваться не стоит, ибо, как говорят, «даже у хорошей хозяйки первый блин — комом».

Шаг 3

Рисунок 17.3. Кривая погрешности при Паде- аппроксимации степени (4,4)
Как и при аппроксимации рядом Тейлора, ошибка здесь мала вблизи точки разложения и велика вдали от нее. Мы снова видим из графика, что для указанной функции, самая большая ошибка — в левой оконечной точке. Однако максимальная ошибка в Паде- аппроксимации уже на порядок меньше, чем при аппроксимации полиномом Тейлора:

Построим графики траекторий для первого

Шаг 4

Построим графики траекторий для первого случая:

Шаг 4

Эти графики показаны на Рисунок 17.10.

Шаг 4

Рисунок 17.18. ФЧХ фильтра на операционном усилителе

Шаг 4

В соответствии с выбранным операторным методом анализа введем обозначения:

Шаг 4

Рисунок 17.20. АЧХ цифрового полосового фильтра
Теперь приступим к тестированию фильтра. Зададим входной сигнал в виде зашумленного меандра с частотой 500 Гц и размахом напряжения 2 В:
> 1 :=round(fs/2/500):
> for n from 0 by 2*1 to Т do
> for n2 from 0 to 1-1 do
> if n+n2 <= Т then
> x[n+n2] := evalf(-l+rand()/10^12-0.5);
> fi:
> if n+n2+1 <= Т then
> x[n+n2-H] :-=eva1f(l+ranoX)/10^12-0.5);
> fi;
> od:
> od:
Временная зависимость синтезированного входного сигнала представлена на Рисунок 17.21.

Шаг 4

Рисунок 17.27. Фазовый портрет колебаний на фоне ВАХ туннельного диода и линии нагрузки резистора Rs
О том, что колебания релаксационные можно судить по тому, что уже первый цикл колебаний вырождается в замкнутую кривую — предельный цикл, форма которого заметно отличается от эллиптической.
Итак, мы видим, что данная цепь выполняет функцию генератора незатухающих релаксационных колебаний. Хотя поставленная задача моделирования цепи на туннельном диоде успешно решена, в ходе ее решения мы столкнулись с проблемой обеспечения малого шага по времени при решении системы дифференциальных уравнений, описывающих работу цепи. При неудачном выборе шага можно наблюдать явную неустойчивость решения.

Шаг 4

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

Графики траекторий полета камня

Шаг 5

Графики траекторий полета камня с массой 500 г представлены на Рисунок 17.6.

Шаг 5

Рисунок 17.10. Траектории движения частиц
Полученные графики (Рисунок 17.10) наглядно показывают на одну из возможностей разделения изотопов. Как говорится, осталось подставить «стаканчик» в нужное место для ловли нужных изотопов. Разумеется, это только изложение идеи одного из методов разделения изотопов. Увы, на практике приходится использовать сложнейшие и дорогие физические установки для решения этой актуальной задачи.

Шаг 5

Это позволяет найти Н как функцию от частоты f также в аналитическом виде:

Шаг 5

Рисунок 17.21. Синтезированный входной сигнал
Вычислим реакцию фильтра на входной сигнал:
> for n from 0 to T do
> y[n] := sum(h[k]*x[n-k],k=0..N);
> od:
Построим график выходного сигнала:
> р := [seq([j/fs, x[j]], j=0..T)]:q:= [seq([j/fs , y[j]] , j =0..Т)]:
> plot(p,time=0..T/fs/4,1abels=[time,volts],title='Входной сигнал\сolor=black);
> plot(q,tine=0..T/fs/4,1abels=[tirae,volts], titlе='Выходной сигнал",color=black);
Временная зависимость выходного сигнала показана на Рисунок 17.22. Нетрудно заметить, что в конце концов выходной сигнал вырождается в пятую гармонику входного сигнала, но этому предшествует довольно заметный переходной процесс. Он связан с узкополосностью данного фильтра.

Теперь построим графики траекторий для

Шаг 6

Рисунок 17.6. Баллистические траектории камня с массой 500 г
Теперь построим графики траекторий для второго случая:
> display({a3,a4,t1},title='Tpaeкт. полета тела массой 100 г, labels=[x.у], labelfont=[TIMES.ROMAN,14]):
Они представлены на Рисунок 17.7.

Шаг 6

Это тоже довольно громоздкое выражение, и его применение при «ручном» анализе потребовало бы от нас немало изобретательности. Между тем Maple 7 позволяет «в два счета» определить из него амплитудно-частотную (AVM) и фазо- частотную (PhaseAV) характеристики усилителя как функции частоты:
> AVM=-evalc(abs(H)):
> PhaseAV:=evalc(argument(H)):
Преобразуем AVD в логарифмическую характеристику, выражающую усиление в децибелах (dB):
> AVdB:=20*1og10(AVM):
Такая характеристика более привычна для специалистов в радиоэлектронике. Соответственно фазо-частотную характеристику выразим в градусах:
> R2D:=evalf(360/(2*Pi));R2D := 57.29577950
> AVdeg:=R2D*PhaseAV:
Теперь можно перейти к обычным численным расчетам. Зададим конкретные значения компонент эквивалентной схемы усилителя:
> Rl:=100: R2:=100000: R3:=1000: R4:=10000: Cl:=1.*10^(-6): С2:=5*10^(-12): СЗ:=1*10^(-6): mu:=50:
Построим амплитудно-частотную характеристику усилителя:
> gaindata:-NULL:
phasedata:=NULL:
for a from 0 to 8 do:
for i from 2*10^a to l(T(a+l) by 10^a do
gaindata:=gaindata, [1. evalf(subs(f=i,AVdB))];
phasedata:=phasedata, [i, eva1f(subs(f=i,AVdeg))]:
od: od:
> 1oglogp1ot([gaindata]. thickness»2, color=black, style=1ine, axes=boxed,
title=`Коэффициент усиления K(f)`,1abels=['Частота (Hz)VK(d8)']):
Она показана на Рисунок 17.13.

Шаг 6

Рисунок 17.22. Временная зависимость выходного сигнала цифрового фильтра
Вычислим спектры входного и выходного сигналов, подготовив массивы выборок сигналов и применив прямое преобразование Фурье с помощью функции FFT:
> Н := array(l..T+l):1i :=array(1..Т+1):
> for n from 0 to T do ,
> ri[n+l] := x[n]*2/T: ii[rn-l] := 0;
> ro[n+l] := y[n]*2/T; Io[rrfl] := 0;
> od:
> FFT(m.ri,ii):rTT(m,ro,io):
Построим график спектра входного сигнала, ограничив масштаб по амплитуде значением 0,5 В:
> р := [seq([j*fs/(T+l),abs(n[j+l]+ii[j-H]*I)],j=0..T/2)]:
> q := [seq([j*fs/(T-H),abs(ro[j-H]+To[j+l]*I)],j=0..T/2)]:
> plot(p, frequency=0..fs/2,y0..0.5,labe1s=[Частотa.V],title='Частотный спектр входного сигнала',color=black);
Этот график представлен на Рисунок 17.23. Из него хорошо видно, что спектральный состав входного сигнала представлен только нечетными гармониками, амплитуда которых убывает по мере роста номера гармоники. Пятая гармоника на частоте 2500 Гц находится посередине полосы пропускания фильтра, ограниченной граничными частотами фильтра 2300 и 2700 Гц. Заметны также беспорядочные спектральные линии шума сигнала в пределах полосы прозрачности фильтра.
Теперь построим график спектра выходного сигнала:
> p1ot(q, frequency=0..fs/2,y=0..0.5,labe1s=[Частотa,V], title='Частотный спектр выходного сигнала'бcolor=black);
Он представлен на Рисунок 17.24. Хорошо видно эффективное выделение пятой гармоники сигнала и прилегающей к ней узкой полосы шумового спектра.

Баллистические траектории камня при массе

Шаг 7

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

Шаг 7

Рисунок 17.13. Амплитудно-частотная характеристика усилителя
Далее зададим построение фазо-частотной характеристики усилителя:
> 1og1ogplot([phasedata], thickness=2, color=b1ue, style=line, axes=boxed, title='Фаэовый сдвиг (в градусах)`, labels=['Частота (Hz)','Фаза']);
Она представлена на Рисунок 17.14.

Шаг 7

Рисунок 17.23. Спектрограмма входного сигнала

Найдем номинальный коэффициент усиления на

Шаг 8

Рисунок 17.14. Фазо-частотная характеристика усилителя
Найдем номинальный коэффициент усиления на частоте f=1000 (Гц):
> AVmid:=eva1f(subs(f=1000, AVdB)):
AVmid=33.12074854
Имея аналитическое выражение для амплитудно-частотной характеристики, можно составить уравнения для вычисления граничных частот (по спаду усиления на -dAV в dB):
> dAV:=3: #Ослабление (в dB на граничных частотах)
> eq5:=AVmid-dAV=20*log10(AVM):
Теперь можно найти эти частоты — нижнюю и верхнюю:
> flow:=fsolve(eq5,f. f-10..2000):flow:= 23.61659476
> fhigh:=fsolve(eqS,f, f-2000..100*10*6);
fliigh := .5737800225 107
Мы можем построить и более наглядную амплитудно-частотную характеристику с точками, соответствующими граничным частотам:
> with(plottools) :h:=log10(AVnvid-dAV):
aplot:= Loglogplot([gaindata], thickness=2, color=b1ack. style=line, axes=boxed,
title='Частоты flow и fhigh среза', labels=['Частота (Hz)VK(dB)']):
bplot:=line([0.1,h], [7.1,h], color=black, linestyle=3):
cplot:=line([log10(flow),0.58],[logHK flow). 1.6], color=blue, linestyle=3):
dplot:=line([log10(fh1gh).0.58],. [log10(fhigh).1.6],. color=red,. 1inestyle=3):
display([aplot.bplot,cplotJ,dplot]):
Эта характеристика показана на Рисунок 17.15.
На ней проставлены синяя и красная пунктирные вертикали, соответствующие найденным граничным частотам flow и fhigh, а также пунктирная горизонталь, соответствующая коэффициенту усиления на этих частотах. Это позволяет наглядно оценить частотный диапазон работы усилителя.
Таким образом, задача расчета усилителя в малосигнальном режиме полностью решена. Мы получили значение номинального коэффициента усиления, рассчитали нижнюю и верхнюю граничные частоты, получили аналитические выражения для амплитудно-частотной и фазо-частотной характеристик усилителя и построили их наглядные графики.

Шаг 8

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

с выделенными точками граничных частот

Шаг 9

Рисунок 17.15. Амплитудно-частотная характеристика с выделенными точками граничных частот

Сравнение времен вычислений

Теперь определим время, необходимое для вычисления функции/(л:) в 1000 точек, используя первоначальное интегральное определение, и сравним его с временем, требующимся для схемы MinimaxApprox в виде непрерывной дроби. Так как наше приближение будет давать только 6 точных цифр, мы также потребуем 6 точных цифр и от интегрального представления функции:

> Digits :=б: St :=time():

> seq( evalf(f(i/250.0)), i = 1..1000 ):

> oldtime := time() - st;

oldtime-81.805

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

> Digits := 8: st := tirae():

> seq( MinimaxApprox(i/250.0), i = 1..1000 ):

> newtime :» time()- st;

newtime:= .694

Ускорение вычисления при аппроксимации есть:

> SpeedUp := oldtime/newtime;

SpeedUp:=U7.S7464

Мы видим, что процедура вычислений, основанная на MinimaxApprox, выполняется почти в 120 раз быстрее процедуры с использованием исходного интегрального определения. Это просто феноменальный успех, полностью оправдывающий время, потерянное на предварительные эксперименты по аппроксимации и ее оптимизации! Разумеется, при условии, что вы будете применять эту аппроксимацию многократно.

Выбор аппроксимации для сложной функции

Задание исходной функции и построение ее графика

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

Воспользуемся возможностями пакета numapprox, для чего прежде всего подключим его:

> restart:with(numapprox):

[chebdeg,chebmult,chebpade,chebsort,chebyshev, confracform,hermite_pade,hornerform, infnorm,laurent,minimax,pade,remez]

Будем искать приемлемую аппроксимацию для следующей, отнюдь не простой, тестовой функции: