Параметры системы и интервал исследования

1.  Параметры системы и интервал исследования

 

egsp_data;

fprintf('-------------------------------------------------------------\n');

fprintf('1. Интервал исследования\n');

fprintf('------------------------\n');

fprintf('tmin = %e, c;\n',tmin);

fprintf('tmax = %e, c;\n',tmax);

fprintf('Nt = %i;\n',Nt);

fprintf('\n');

2.  Формирование системы базисных функций

 

settime(T);

setsize(Nt);

Ai = mkint;

Ad = inv(Ai);

Ae = eye(Nt);

fprintf('-------------------------------------------------------------\n');

fprintf('2. Формирование системы функций Уолша\n');

fprintf('-------------------------------------\n');

%pr_matrix(Ai,'Оператор интегрирования Ai');

disp('Оператор интегрирования Ai');

disp(Ai(1:8,1:8));

%pr_matrix(Ad,'Оператор дифференцирования Ad');

disp('Оператор дифференцирования Ad');

disp(Ad(1:8,1:8));

3.  Расчет операторов левой линейной части

 

fprintf('-------------------------------------------------------------\n');

fprintf('3. Операторы левой линейной части\n');

fprintf('---------------------------------\n');

% оператор ПФ W1(s) - электрической части

Aw1 = inv(RS*(Ty*Ae+Ai))*(Ky1*KU*Ai);

%pr_matrix(Aw1,'Оператор Aw1');

disp('Оператор Aw1');

disp(Aw1(1:8,1:8));

% оператор ПФ W2(s) - электромагнитного преобразователя и часть расходов

Aw2 = inv(CS*(Tem^2*Ae+2*Tem*dzem*Ai+Ai^2))*(KFi*Kqh*Ai^2);

%pr_matrix(Aw2,'Оператор Aw2');

disp('Оператор Aw2');

disp(Aw2(1:8,1:8));

% оператор ПФ W3(s) - движения золотника и часть расходов

Aw3 = inv(Kqp1*(Cp+Cg)*(Tz^2*Ae+2*Tz*dzz*Ai+Ai^2))*(Az*Ai^2);

%pr_matrix(Aw3,'Оператор Aw3');

disp('Оператор Aw3');

disp(Aw3(1:8,1:8));

% оператор ПФ W4(s) - местной обратной связи

Aw4 = Az*Ad;

%pr_matrix(Aw4,'Оператор Aw4');

disp('Оператор Aw4');

disp(Aw4(1:8,1:8));

% оператор левой линейной части

Aw34 = inv(Ae+Aw4*Aw3)*Aw3;

Aw_1 = Aw34*Aw2*Aw1;

%pr_matrix(Aw_l,'Оператор левой части Aw_l');

disp('Оператор левой части Aw_l');

disp(Aw_1(1:8,1:8));

4.  Расчет операторов правой линейной части

 

fprintf('-------------------------------------------------------------\n');

fprintf('4. Операторы правой линейной части\n');

fprintf('----------------------------------\n');

% оператор ПФ W5(s) - уравнения расходов

Aw5 = inv(Kqp*(Tg*Ae+Ai))*(Ap*l*Ai);

%pr_matrix(Aw5,'Оператор Aw5');

disp('Оператор Aw5');

disp(Aw5(1:8,1:8));

% оператор ПФ W6(s) - нагрузка

Aw6 = J*Ai^2;

%pr_matrix(Aw6,'Оператор Aw6');

disp('Оператор Aw6');

disp(Aw6(1:8,1:8));

% оператор ПФ W7(s) - трение

Aw7 = Kf*Ad;

%pr_matrix(Aw7,'Оператор Aw7');

disp('Оператор Aw7');

disp(Aw7(1:8,1:8));

% оператор ПФ W8(s) - местная обратная связь

Aw8 = Ap*l*Ad;

%pr_matrix(Aw8,'Оператор Aw8');

disp('Оператор Aw8');

disp(Aw8(1:8,1:8));

% оператор правой линейной части

Aw67 = inv(Ae+Aw7*Aw6)*Aw6;

Aw671 = inv(Ae+Ksh*Aw67)*Aw67;

Aw_r = Kz*inv(Ae+Aw8*Aw671*Aw5)*(Aw671*Aw5);

%pr_matrix(Aw_r,'Оператор правой части Aw_r');

disp('Оператор правой части Aw_r');

disp(Aw_r(1:8,1:8));

5.  Спектральная характеристика входного сигнала

 

fprintf('-------------------------------------------------------------\n');

fprintf('5. Спектральная характеристика входного сигнала\n');

fprintf('-----------------------------------------------\n');

u = zeros(1,Nt)+1;

Cu = fwht(u');

%pr_matrix(Cu,'Cu');

disp('Оператор Cu');

disp(Cu(1:8));

6.  Расчет выходного сигнала методом последовательных приближений

 

fprintf('-------------------------------------------------------------\n');

fprintf('6. Расчет выходного сигнала\n');

fprintf('---------------------------\n');

Cd_old = zeros(Nt,1);

Ce = Cu-Cd_old;

Cx = Aw_1*Ce;

x = iwht(Cx)';

xf = egsp_f(x,xm);

Cxf = fwht(xf');

Cd_new = Aw_r*Cxf;

Ceps = Cd_new-Cd_old;

Eeps = sqrt(Ceps'*Ceps);

fprintf('Нулевое приближение\n');

fprintf('Eeps = %e;\n',Eeps);

d = iwht(Cd_new)';

figure; clf;

plot(t,d);

xlabel('t, c');

ylabel('delta(t)');

Niter = 0;

while Eeps >= 1e-8

Niter = Niter+1;

Cd_old = Cd_new;

Ce = Cu-Cd_old;

Cx = Aw_1*Ce;

x = iwht(Cx)';

xf = egsp_f(x,xm);

Cxf = fwht(xf');

Cd_new = Aw_r*Cxf;

Ceps = Cd_new-Cd_old;

Eeps = sqrt(Ceps'*Ceps);

end

fprintf('%i-ое приближение:\n',Niter);

fprintf('Eeps = %e;\n',Eeps);

d = iwht(Cd_new)';

%my_plot2(t,d,'t, c','delta(t)');

plot(t,d);

xlabel('t, c');

ylabel('delta(t)');

grid on;

toc;