Главная страница Случайная страница
Разделы сайта
АвтомобилиАстрономияБиологияГеографияДом и садДругие языкиДругоеИнформатикаИсторияКультураЛитератураЛогикаМатематикаМедицинаМеталлургияМеханикаОбразованиеОхрана трудаПедагогикаПолитикаПравоПсихологияРелигияРиторикаСоциологияСпортСтроительствоТехнологияТуризмФизикаФилософияФинансыХимияЧерчениеЭкологияЭкономикаЭлектроника
⚡️ Для новых пользователей первый месяц бесплатно. А далее 290 руб/мес, это в 3 раза дешевле аналогов. За эту цену доступен весь функционал: напоминание о визитах, чаевые, предоплаты, общение с клиентами, переносы записей и так далее.
✅ Уйма гибких настроек, которые помогут вам зарабатывать больше и забыть про чувство «что-то мне нужно было сделать».
Сомневаетесь? нажмите на текст, запустите чат-бота и убедитесь во всем сами!
Интегрирование дифференциального уравнения.
|
|
Сравнение с расчетами, выполненными в Turbo Delpi
> restart;
> with(linalg): with(plots):
pp: =(x, y)-> [x, y];
Задаем модельный сигнал.
Это финитная функция 
единичной амплитуды, заданная на единичном отрезке-носителе.
> fun: = proc(t) local z; z: =piecewise(t< 0, 0, t< 1/2, 2*t, t< =1, 2-2*t, 0); evalf(z); end;
> plot(fun(t), t=-1..2, thickness=2, color=brown);
Полином
Первый положительный корень полинома определяет период T управляющего сигнала и коэффициент его усиления.
Величина 
определяет длительность сигнала на участке периода.
> p(x): =x^5-7*x-6;
> Koeff: =fsolve(p(x), x, 0..2);
> T: =Koeff;
> tau: =1;
Процедура Period(t, t0, tau, T, f) строит периодическое продолжение сигнала на всю действительную ось, где T - период, 
- длительность сигнала f(t), t0 -- начальный сдвиг (который можно интерпретировать как запаздывание).
> Period: =proc(t, t0, tau, T, f) local x, z;
x: =evalf(t-t0-floor((t-t0)/T)*T);
z: =fun(x/tau); evalf(z);
end;
> plot(Period(x, 0, tau, T, fun), x=-1..3, thickness=2, color=brown);
> #==================================================================
Решаем задачу интегрирования входного сигнала
Экспериментально так, чтобы выходной сигнал был близким к стационарному.
> Koc: =1; Nzac: =7;
> ur: =diff(U(t), t);
> F: =Nzac*(cos((4+Nzac/10)*t+U(t))+Koeff*Period(t, 0, tau, T, f)-Koc*U(t));
> RK: =dsolve({ur=F, U(0)=0.25}, U(t), type=numeric, output=listprocedure);
> fU: =subs(RK, U(t));
> T0: =5; Nt: =50; h: =T0/Nt;
> Tx: =array(0..Nt): U: =array(0..Nt): U_map: =array(0..Nt);
> for j from 0 to Nt do
x: =j*h; z: =fU(x); Tx[j]: =x; U[j]: =z; U_map[j]: =z;
#print(x, z);
od:
> RisU: =zip(pp, Tx, U):
> RU: =plot(RisU):
> display(RU);
> #====================================
> RisU: =zip(pp, Tx, U):
> whattype([RisU]);
> RU0: =plot(RisU, style=point, symbol=cross):
> display(RU0);
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
> fn: =`D: \\wrem\\сава.txt`;
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
> L: =readdata(fn, 2);
Nstrok: =vectdim(L);
> U_n: =array(1..Nstrok);:
T_n: =array(1..Nstrok);
> for j from 1 to Nstrok do
T_n[j]: =L[j, 1];
U_n[j]: =L[j, 2];
#print(j, T_n[j], U_n[j]);
od:
> u1: =zip(pp, T_n, U_n):
> RU1: =plot(u1, style=point, symbol=cross, color=black):
> display(RU, RU1);
Сравнение результатов на одинаковых сетках
> printf(" %s", ` № t U_map U_pas разн \n`);
for k from 0 to Nt do t: =Tx[k]: del: =U_map[k]-U_n[k+1];
printf(" % 3.0f % 6.2f % 8.4f % 8.4f % 8.4f \n", k, t, U_map[k], U_n[k+1], del):
end:;
№ t U_map U_pas разн
0 0.00 0.2500 0.2500 0.0000
1 0.10 0.5902 0.5429 0.0473
2 0.20 0.6686 0.5954 0.0732
3 0.30 0.6499 0.5914 0.0585
4 0.40 0.6294 0.5946 0.0348
5 0.50 0.6687 0.6538 0.0149
6 0.60 0.6394 0.5972 0.0422
7 0.70 0.5007 0.4635 0.0372
8 0.80 0.3497 0.3170 0.0327
9 0.90 0.2349 0.2033 0.0316
10 1.00 0.1875 0.1557 0.0318
11 1.10 0.3274 0.3055 0.0219
12 1.20 0.6228 0.6312 -0.0084
13 1.30 0.7754 0.7561 0.0193
14 1.40 0.7168 0.6719 0.0449
15 1.50 0.5445 0.5026 0.0419
16 1.60 0.3284 0.2978 0.0306
17 1.70 0.0986 0.0785 0.0201
18 1.80 -0.1312 -0.1435 0.0123
19 1.90 -0.2576 -0.2278 -0.0298
20 2.00 -0.2713 -0.2425 -0.0288
21 2.10 -0.1574 -0.1330 -0.0244
22 2.20 0.2102 0.2274 -0.0172
23 2.30 1.0514 1.0792 -0.0278
24 2.40 1.7867 1.7752 0.0115
25 2.50 1.8327 1.6999 0.1328
26 2.60 1.5513 1.4331 0.1182
27 2.70 1.1729 1.0939 0.0790
28 2.80 0.7690 0.7202 0.0488
29 2.90 0.4490 0.4590 -0.0100
30 3.00 0.1933 0.2083 -0.0150
31 3.10 -0.0450 -0.0339 -0.0111
32 3.20 -0.2728 -0.2650 -0.0078
33 3.30 -0.4857 -0.4793 -0.0064
34 3.40 -0.6743 -0.6681 -0.0062
35 3.50 -0.8233 -0.8163 -0.0070
36 3.60 -0.9029 -0.8972 -0.0057
37 3.70 -0.7202 -0.7068 -0.0134
38 3.80 -0.0655 -0.0543 -0.0112
39 3.90 0.9012 0.9790 -0.0778
40 4.00 1.3895 1.2936 0.0959
41 4.10 1.4820 1.3657 0.1163
42 4.20 1.2741 1.1247 0.1494
43 4.30 0.9205 0.8169 0.1036
44 4.40 0.5293 0.4679 0.0614
45 4.50 0.1288 0.0950 0.0338
46 4.60 -0.2709 -0.2909 0.0200
47 4.70 -0.5622 -0.5425 -0.0197
48 4.80 -0.7563 -0.7355 -0.0208
49 4.90 -0.8764 -0.8591 -0.0173
50 5.00 -0.8943 -0.8795 -0.0148
5. Спектральный анализ программы средствами математического пакета MAPLE
> restart;
> with(linalg): with(plots):
pp: =(x, y)-> [x, y];
> f: = proc(t) local z; z: =piecewise(t< 0, 0, t< =1, 1, 0); evalf(z); end;
> plot(f(t), t=-1..2, thickness=2, color=brown);
> T: =2*Pi; #период сигнала
> tau: =2; #длительность сигнала
Для полученной периодической функции вычисляем коэффициенты тригонометрического ряда Фурье
в соответствии с формулами
Коэффициенты Фурье вычисляются в аналитическом виде процедурой Fourier(T, F), где параметрами процедуры являются величина периода T и имя процедуры, вычисляющей значения сигнала на отрезке периодичности. В нашем случае это процедура Period, для использования, которой следует согласовать параметры обеих процедур. Процедура Fourier вычисляет коэффициенты Фурье и передает их как глобальные для процедуры объекты. Интегрируемая функция имеет один параметр - переменную интегрирования. При этом функция должна быть определена на отрезке периодичности, поэтому следует задать процедуру, в которой не должно быть других параметров:
> signal: =proc(x)global tau, T;
f(x/tau);
end;
> plot(signal(x), x=0..T, thickness=2, color=blue);
Таким образом, фактическим параметром, соответствующим названию функции при вызове процедуры Fourier должна быть процедура F.
> Fourier: =proc(T, f) global a0, A, B;
a0: =2/T*evalf(Int(f(x), x=0..T)):;
A: =value(2/T*int(f(x)*cos(2*Pi*k*x/T), x=0..T));
B: =value(2/T*int(f(x)*sin(2*Pi*k*x/T), x=0..T));
end proc:
Вычислим выражения для коэффициентов и выведем их значения:
> Fourier(T, signal):
> a0; evalf(A, 4); evalf(B, 4);
Значения коэффициентов 
и 
можно вычислить подстановкой (процедура subs)
натуральных чисел в полученные выражения
Зафиксируем число коэффициентов (степень)тригонометрического полинома
Фурье и построим совместные графики исходного сигнала и приближающих его
полиномов для нескольких значений степеней. Глобальный параметр k, как и выражения для коэффициентов
Фурье определяются процедурой Fourier
> Trig: =proc(t, n) local z; global k, T, a0, A, B;
z: =a0/2+sum(A*cos(k*2*Pi*t/T)+B*sin(k*2*Pi*t/T), k=1..n);
end:
Построим график триг. полинома пятой степени
> plot(Trig(x, 5), x=-T..2*T);
> grafik: =plot(signal(x), x=-0.2..T+0.2, thickness=3, color=blue, discont=true):
>
Построим совместные графики сигнала и полиномов нескольких степеней
> CL: =brown, red, gold;
> RisTrig: =seq(plot(Trig(t, 2*j), t=-0.1..2*T, color=CL[j]), j=1..3):
> display(RisTrig, grafik);;
ЧАСТЬ 2
> with(plottools):
Выполним спектральный анализ аналогового сигнала. Затем дискретизируем сигнал и также построим его спектр для сравнения и понимания различий в спектрах дискретного и аналогового сигналов. Присоединим процедуры, решающие эту задачу.
> Fourier_T: =proc(F, T0, TN, k:: evaln) local T;
global A0, Ak, Bk;
T: =TN-T0;
A0: =2/T*Int(F(x), x=T0..TN);
Ak: =2/T*int(F(x)*cos(k*x*2*Pi/T), x=T0..TN):
Bk: =2/T*int(F(x)*sin(k*x*2*Pi/T), x=T0..TN):
end proc:
Процедура Trig_polynom подстановкой номеров k в формулы коэффициентов вычисляет эти коэффициенты и сохраняет во вновь созданных массивах a, b [0..N]. Заодно строится график тригонометрического полинома степени N.
> Trig_polynom: =proc(N, T, a0, ak, bk, k:: evaln) local n, Pol, A0, A, B;
global a, b, RisTrig;
a: =array(0..N); b: =array(0..N);
A0: =evalf(a0); a[0]: =A0; b[0]: =0;
A: =seq(evalf(subs(k=n, ak)), n=1..N);
B: =seq(evalf(subs(k=n, bk)), n=1..N);
for n from 1 to N do
a[n]: =A[n]; b[n]: =B[n];
end do;
Pol: =A0/2+sum(A[k]*cos(2*Pi*k*x/T)+B[k]*sin(2*Pi*k*x/T), k=1..N):
RisTrig: =plot(Pol, x=-T/2..3*T/2, color=blue, thickness=2):
RETURN(Pol);
end proc:
Процедура ARR строит спектр амплитуд линейчатой гистограммой.
> ARR: =proc(n:: integer, c) local L, H, ma, mi, k:: integer,
Sim:: array(0..n);
ma: =c[0]; mi: =c[0];
L: =line([0, c[0]], [n, c[n]], thickness=2, color=red);
for k from 1 to n do
if c[k]> ma then ma: =c[k]; end if;
if c[k]< mi then mi: =c[k]; end if;
end do;
H: =ma-mi;
if H=0 then RETURN(L) end;
for k from 0 to n do
if abs(c[k])< H/1000 then
Sim[k]: =ellipse([k, c[k]], 0.2, 0.01*H, color=blue);
else
Sim[k]: =plottools[arrow]([k, 0], [k, c[k]], 0.2, 0.2, 0, color=black);
end if;
end do;
convert(Sim, list);
end:
Процедура Spectr вычисляет модули и аргументы комплексной (синусной) формы тригонометрического полинома.
График спектра фаз - структура Risphi.
> Spectr: =proc(n, a, b, c, Risphi) local k, R, phi;
for k from 0 to n do
c[k]: =evalf(abs(I*a[k]+b[k])):
phi: =evalf(argument(I*a[k]+b[k]));
R[k]: =[eval(k), eval(phi)];
end:;
Risphi: =plot(convert(R, list)):
end:
Определим прямоугольный импульс f(t) единичной длительности.
Построим сигнал (F_for_all) длительности и для него выполним вычисления спектральных характеристик.
Здесь и дальше команды не будут зависеть от варианта сигнала, поскольку название обрабатываемого сигнала теперь у всех одинаково (F_for_all)
Положим величину периода T=2 
(для всех вариантов). Аналитические выражения коэффициентов вычисляются в процедуре Fourier_T.
> tau: =Pi; # длительность сигнала
> T: =2*Pi; # величина периода
> F_for_all: =proc(t) global tau; f(t/tau); end proc:;
> Ris1: =plot(F_for_all(t), t=0..T, color=brown, thickness=2): display(Ris1);
Получим формулы для коэффициентов Фурье, зависящие от номера коэффициента k.
> Fourier_T(F_for_all, 0, T, k):
> a0: =evalf(A0); a_k=Ak; b_k=Bk;
Значения величин С[k] в аналитическом виде можно вычислить (для любого натурального k) по формуле:
> C[k]: =simplify(sqrt(Ak^2+Bk^2));
Для построения спектра амплитуд зафиксируем число коэффициентов и с помощью процедуры Trig_polynom вычислим коэффициенты a[k], b[k] тригонометрического полинома степени N, а также построим его график на отрезке, большем длины периода. Массивы коэффициентов создаются при работе процедуры и являются глобальными, то есть доступными вне процедуры.
> N: =19;
> Trig_polynom(N, T, A0, Ak, Bk, k):
> display(RisTrig, Ris1);
Амплитудный спектр (коэффициенты С[k]) вычисляется процедурой Spectr, которая как побочный эффект строит графическую структуру Risphi1, являющуюся графиком спектра фаз. График спектра амплитуд строит процедура ARR.
> Spectr(N, a, b, c, 'Risphi1'):;
> display(ARR(N, c));
> RisAMP: =display(ARR(N, c)):
> display(Risphi1);
ЧАСТЬ 3
Дискретизация задачи
Пакет MAPLE выполняет вычисления в аналитической форме, если не указано дополнительно требование: получить результат в числовом виде, например команда fsolve. Для получения результатов, совпадающих с вычислениями, выполненными программным путем (среда Turbo Delphi), перейдем к дискретному варианту решения задачи, то есть к вычислениям на конечном множестве точек отрезка [0, T]. Определим массив Y для сеточных значений функции F_for_all(t), соответственно. Определим также массивы a, b для размещения значений коэффициентов тригонометрического полинома, аппроксимирующего функцию f(t) на отрезке длины периода T > .
Параметр Nt, задающий число дискретных отсчетов (0.. Nt) следует по изменять и выбрать каких-нибудь 2-3 варианта. Этот параметр определяет качество аппроксимации сигнала на периоде.
Число коэффициентов тригонометрического полинома выбирается как максимально возможное и равное n = Nt/2. Массивы коэффициентов Фурье определяются по максимальному значению n, но в силу ортогональности системы тригонометрических функций, тригонометрические полиномы меньшей степени имеют те же самые коэффициенты.
Присоединим модифицированную процедуру DTF из пособия, а также процедуру Setka и Spectr_DTF, последняя из которых вычисляет модули и аргументы комплексных коэффициентов ДПФ.
> DTF: =proc (y, N, Y, n) local k, j, p, h;
h: =2*Pi/N;
2.1: for k from 0 to n do
p: =0;
for j from 0 to N-1 do
p: =p+evalf(y[j]*exp(-I*k*j*h));
end;
Y[k]: =evalf(p/N);
end:
end:;
> Setka: =proc(Nt, T, F, Y:: array) local h, j, x, R;
global GrafF;
h: =T/Nt;
for j from 0 to Nt do
x: = evalf(j*h);
Y[j]: = F(x);
R[j]: =[j, eval(Y[j])];
end:
5.1: R[Nt]: =[j, eval(Y[0])];
GrafF: =plot(convert(R, list), 0..Nt, color=brown,
style=point, symbol=circle):
end:
> 6: Spectr_DTF: =proc(n, C, A, phi) local k, R; global Risphi;
6.1: for k from 0 to n do
A[k]: =evalf(abs(C[k])):
phi[k]: =evalf(argument(C[k]));
R[k]: =[eval(k), eval(phi[k])];
end:;
Risphi: =plot(convert(R, list), thickness=2, color=blue, style=point, symbol=box):
end:
Зададим число точек дискретизации N= Nt+1, число коэффициентов n в массивах ДПФ (N- параметр процедур ДПФ и он равен числу дискретных отсчетов сигнала, которые нумеруются индексами 0..N-1).
Для демонстрации свойства симметрии коэффициентов n будем задавать большим Nt. В случае вещественного сигнала (не имеющего мнимой составляющей) в силу свойства симметрии спектральных характеристик (четная - для АЧХ, нечетная - для ФЧХ) полагают n=Nt/2, если требуется только построение их графиков, а не восстановление сигнала обратным дискретным преобразованием Фурье.
Параметры задачи
> Nt: =35; `число дискретных отсчетов -1`:
> n: =Nt; # в обычном случае
n: =57; # для демонстрации свойства симметрии
> C: =array(0..n): phi: =array(0..n): A: =array(0..n):;
Y: =array(0..Nt):
> Setka(Nt, T, F_for_all, Y);
> DTF(Y, Nt+1, C, n):
> Spectr_DTF(n, C, A, phi):;
> display(ARR(n, A));
Таким образом, АЧХ и ФЧХ являются периодическими дискретными последовательностями с периодом N=32 = Nt+1.
АЧХ обладает четной симметрией относительно середины N/2=16, а ФЧХ - нечетной симметрией.
Сравним графики спектров амплитуд аналогового и дискретного сигналов.
> display(RisAMP);
Заключение
В процессе выполнения данной курсовой работы я ознакомилась с языками программирования, как Turbo Delphi, MAPLE.
Написала программы на языках программирования для решения заданного полинома. Сравнила и оценила погрешность выполненных расчетов.
Библиография
1.Кондратьев В.П. Языки программирования. Система Maple.ч. I. Основы работы в системе. Ч. III/. Язык программирования системы. Учебное пособие. Екатеринбург: УрТИСИ ГОУ ВПО «СибГУТИ», 2006.
2.Мудров А.Е. Численные методы для ПЭВМ на языках БЕЙСИК, ФОРТРАН, ПАСКАЛЬ. – Томск. 1991.
3.Немнюгин С. Turbo Pascal. Программирование на языке высокого уровня. Учебник для вузов. 2-е изд.—СПб.: Питер, 2005, 544с.
4.Фаронов В.В. DELPHI. Программирование на языке высокого уровня. Учебник для вузов.—СПб.: Питер, 2005, 640с.
|
|