Интегрирование дифференциального уравнения.

Сравнение с расчетами, выполненными в 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с.