Урок 17.

Примеры решения научно-технических задач

 

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

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

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

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

 

 

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

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

Трудно представить себе область более широкую и>й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]

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

График этой функции представлен на рис. 17.1. С первого взгляда — это простой график, но тут как раз тот случай, когда простота обманчива. Вы сразу Заметите, что график строится необычно медленно, поскольку в каждой из множества его точек системе Maple 7 приходится вычислять значение интеграла с подынтегральной функцией, содержащей довольно каверзную гамма- функцию. И делает это Maple 7 по сложному и медленному алгоритму адаптивного численного интегрирования.

Рис. 17.1. График аппроксимируемой функции

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

 

 

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

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

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

Кривая ошибок для аппроксимации полиномом Тейлора строится командой: 

> plotd(f- TaylorApprox,0..4,.co1or=black);

и имеет вид, представленный на рис. 17.2. Эта кривая нас, прямо скажем, не слишком радует, поскольку погрешность в сотни раз превышает заданную.

Рис. 17.2. Кривая погрешности при аппроксимации рядом Тейлора

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

> maxTaylorError := abs( Limit(f(x), х-0) - ТауlorАрргох(0) );

 maxTaylorError := .0015029620

Итак, в самом начале наших попыток мы потерпели полное фиаско. Но отчаиваться не стоит, ибо, как говорят, «даже у хорошей хозяйки первый блин — комом».

 

 

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

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

Кривая ошибки для интервала [0, 4] строится командой:

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

Рис. 17.3. Кривая погрешности при Паде- аппроксимации степени (4,4)

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

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

 

 

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

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

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

Можно проверить для этого примера, что кривая ошибки при аппроксимации рядом Чебышева колеблется. Поскольку ряд Чебышева был оборван на члене 8-й степени (как и полином ряда Тейлора), то максимальная ошибка оказалась равной приблизительно 0,6* 10^-5. Эта величина уже на два порядка меньше, чем ошибка при Паде- аппроксимации, вычисленная выше. Но все же немного не дотягивает до наших требований.

Для последующих вычислений полезно заметить, что мы можем использовать процедуру для нахождения численных значений f(x), которая будет намного эффективнее, чем прямое определение, которое требует численного интегрирования для каждого значениях. А именно определим процедуру численной оценки, основанную на разложении в ряд Чебышева степени 13, так как максимальная ошибка  при такой аппроксимации меньше чем 10~⁸, и обеспечивает для нашей цели достаточную точность. Мы определим полином Чебышева Т(х) из пакета orthopoly и затем для эффективной оценки преобразуем его в форму Горнера:

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

 

 

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

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

 Построим кривую ошибок:

> with(orthopoly, Т):

> plot(F = ChebPadeApprox, 0..4,color=black):

Она представлена на рис. 17.4.

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

Рис. 17.4. Кривая ошибки при Паде-Чебышева рациональной аппроксимации

> maxChebPadeError :=abs( F(0) - ChebPadeApprox(O) );

maxChebPadeError= .1236746 10^-5

Мы достигли впечатляющего успеха и остается сделать еще один шаг в направлении повышения точности аппроксимации.

 

 

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

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

Максимальная ошибка в аппроксимации MinimaxApprox дается значением переменной maxerror. Заметим, что мы наконец достигли нашей цели получения аппроксимации с ошибкой меньшей, чем 1*10^-6:

> maxMinimaxError := maxerror;

maxMinimaxError := .585025375366 10^-6

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

> plot(F = MinimaxApprox,0..4,color=black):

График ошибки, представленный на рис. 17.5, показывает равные по амплитуде колебания.

 Рис. 17.5. График ошибки при минимаксной аппроксимации

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

 

 

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

Полиномы числителя и знаменателя в минимаксной аппроксимации уже выражены в форме Горнера (то есть в форме вложенного умножения). Оценка полиномом степени п в форме Горнера при 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))))

 

 

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

Теперь определим время, необходимое для вычисления функции/(л:) в 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 раз быстрее процедуры с использованием исходного интегрального определения. Это просто феноменальный успех, полностью оправдывающий время, потерянное на предварительные эксперименты по аппроксимации и ее оптимизации! Разумеется, при условии, что вы будете применять эту аппроксимацию многократно.

 

 

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

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

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

 

 

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

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

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

Итак, пусть подвернувшиеся под руку камни с массой 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 системы дифференциальных уравнений (ДУ). Каждая система состоит из двух ДУ второго порядка и вид этих систем известен из курса физики. Ниже представлено задание этих систем ДУ (для первой системы дан вывод ее вида):

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

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

Создадим графические объекты — результаты решения систем ДУ:

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

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

Рис. 17.6. Баллистические траектории камня с массой 500 г

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

> display({a3,a4,t1},title='Tpaeкт. полета тела массой 100 г, labels=[x.у], labelfont=[TIMES.ROMAN,14]):

Они представлены на рис. 17.7.

Рис. 17.7. Баллистические траектории камня при массе 100 г

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

 

 

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

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

> 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 имеют вид:

Зададим исходные числовые данные (опустив размерности):

> 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 имеет характерную колебательную компоненту. Случай явно куда менее тривиальный, чем полет камня, описанный выше.

Рис. 17.8. Траектория движения частицы в магнитном поле

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

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

 

 

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

Рассмотрим еще одну классическую задачу ядерной физики — разделение изотопов (атомов с одинаковым зарядом ядра, но разной массой). Для этого используют различные способы. В частности, это может быть масс-спектроскопический метод. Из точки А вылетают однозарядные ионы (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

Рис. 17.9. Иллюстрация к методу разделения изотопов

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

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

> 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:

Выполним решение составленной выше системы дифференциальных уравнений:

Построим графики решения:

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

Рис. 17.10. Траектории движения частиц

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

 

 

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

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

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

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

Введем исходные числовые данные для вычислений:

> 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. Этот график настолько нагляден, что не требует пояснения.

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

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

 

 

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

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

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

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

 

 

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

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

Рис. 17.12. Принципиальная (о) и эквивалентная(6) схемы усилителя на полевом транзисторе

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

> restart:with(plots):

Warning, the name changecoords has been redefined

Из законов Киргофа вытекает, что сумма токов, втекающих в каждый узел и вытекающих из него равна 0. Следовательно, для узлов эквивалентной схемы рис. 17.12 можно записать следующую систему уравнений в операторной форме:

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

> solve({eql,eq2,eq3.eq4}б{Vl,V2.V3,Vo}):

Обеспечим присвоение переменным Vo, VI, V2 и V3 найденных из решения системы уравнений значений:

> assign(%):

Теперь найдем операторную передаточную функцию в аналитическом виде:

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

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

Это тоже довольно громоздкое выражение, и его применение при «ручном» анализе потребовало бы от нас немало изобретательности. Между тем 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.

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

Далее зададим построение фазо-частотной характеристики усилителя:

> 1og1ogplot([phasedata], thickness=2, color=b1ue, style=line, axes=boxed, title='Фаэовый сдвиг (в градусах)`, labels=['Частота (Hz)','Фаза']);

Она представлена на рис. 17.14.

Рис. 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, а также пунктирная горизонталь, соответствующая коэффициенту усиления на этих частотах. Это позволяет наглядно оценить частотный диапазон работы усилителя.

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

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

 

 

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

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

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

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

> restart:

Зададим основные уравнения, описывающие работу фильтра на малом сигнале:

Введем круговую частоту:

> 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. Здесь полезно обратить внимание на то, что спад усиления на низких и высоких частотах происходит довольно медленно из-за малого порядка фильтра.

Рис. 17.17. АЧХ фильтра на операционной усилителе

Далее построим фазо-частотную характеристику фильтра как зависимость фазы в радианах от частоты f в Гц:

> plot ([log10(f),phase, f=10..50000], color=black, title=*Фазо-частотная характеристика фильтра*);

Фазо-частотная характеристика (ФЧХ) фильтра показана на рис. 17.18

На ФЧХ фильтра можно заметить характерный разрыв, связанный с превышением фазовым углом граничного значения я. Такой способ представления фазового сдвига общепринят, поскольку его изменения стремятся вписать в диапазон от -я до п.

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

 

 

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

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

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

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

> restart:with(p1ots):

Warning, the name changecoords has been redefined

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

> N := 64: # Число секций фильтра (на 1 меньше порядка фильтра)

> fs:= 10000: # Частота квантования

> fl = 2300: # Нижняя граничная частота

> fh = 2700: # Верхняя граничная частота

> m := 10: # 2^m > N - число точек для анализа

Вычислим:

Вычислим 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 различаются несильно. В этом случае полосовой фильтр по своим свойствам приближается к резонансному, хотя само по себе явление резонанса не используется.

Рис. 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 можно получить АЧХ с довольно плоской вершиной и резкими спадами (говорят, что такая характеристика приближается к прямоугольной).

Рис. 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.

Рис. 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. Нетрудно заметить, что в конце концов выходной сигнал вырождается в пятую гармонику входного сигнала, но этому предшествует довольно заметный переходной процесс. Он связан с узкополосностью данного фильтра.

Рис. 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. Хорошо видно эффективное выделение пятой гармоники сигнала и прилегающей к ней узкой полосы шумового спектра.

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

Рис. 17.24. Спектрограмма выходного сигнала цифрового фильтра

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

 

 

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

А теперь займемся моделированием явно нелинейной цепи. Выполним его для цепи, которая состоит из последовательно включенных источника напряжения 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 

Рис. 17.25. ВАХ туннельного диода

Пусть задано Es = 0,35 В, Rs= 15 Ом, С = 10*10-12, L = 30*10-9 и максимальное время моделирования tm=10*10-9. Итак, задаем исходные данные:

> Es:=.35:Rs:=15:C:=10*10^(-12):L:=30*10^(-6):tm:=10*10^(-9):

 Составим систему дифференциальных уравнений цепи и выполним ее решение с помощью функции dsolve:

Поскольку заведомо известно, что схема имеет малые значения 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. Из них хорошо видно, что цепь создает автоколебания релаксационного типа. Их форма сильно отличается от синусоидальной.

Рис. 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.

Рис. 17.27. Фазовый портрет колебаний на фоне ВАХ туннельного диода и линии нагрузки резистора Rs

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

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

 

 

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

Вернемся к линейным цепям и рассмотрим еще один полезный метод расчета электрических цепей — с помощью интеграла Дюамеля. При нем можно рассчитать временную зависимость выходного напряжения 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). Они показаны внизу графика.

 

 

Что нового мы узнали?

В этом уроке мы научились:

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