Исследование свойств устройств с переменными магнитодвижущими силами, страница 5

% Формирование анализируемой характеристики:

x=lfile0(nstart:nfinish);                           %  Значения по абсциссе

y=lfile0(nfinish+1:nfinish+ni);               %  Значения по ординате

t=i0;                                                         %  Переприсвоение значений абсцисс

s=y;                                                         %   Переприсвоение значений ординат

s(1)=0;

N=length(t);     k=N;                                %  Длина вектора t

% Организация интерполяции сплайном

sh=1;                                                       % Счетчик

        % Новые значения абсцисс с учетом известного k

t1(1)=t(1);                                              % Переприсвоение 1 точки

for j=2:N

   sti0=(i0(j)-i0(j-1))/k;                           % Шаг

   for jv=1:k

      sh=sh+1;

      t1(sh)=t1(sh-1)+sti0;                        % Точки интерполянта по оси абсцисс

   end

end

% Сплайн-интерполяция характеристики без использования возможности

% сглаживания, но представления результатов в полиномиальной форме

spl=csaps(t,s,1);                                     

% 1-я производная в узлах интерполяции характеристики

dsdt=fnder(spl);

dsdt=fnval(dsdt,t1);

% 2-я производная в узлах интерполяции характеристики

dsdt2=fnder(spl,2);

dsdt2=fnval(dsdt2,t1);

figure(2);

hold on;

plot(t1,dsdt,'b');

hold off;

grid on;

XLABEL('X')

YLABEL('1-я производная от Y')

figure(3);

hold on;

yy2=spline(t1,dsdt2);

%plot(t,yy2,'k'); Хочешь, проверь что получится

fnplt(yy2,'k')

hold off;

grid on;

XLABEL('X')

YLABEL('2-я производная от Y')

figure(1);

hold on;

plot(t,s,'o',t1,fnval(spl,t1),'k');

hold off;

grid on;

XLABEL('Значение "Х"')

YLABEL('Значение "Y"')

LEGEND('- базовые узлы','- исходная характеристика.',...

   'РАССТАНОВКУ НОВЫХ УЗЛОВ ', ...

   'ВЫПОЛНЯТЬ НАЧИНАЯ С НУЛЯ.', ...

   'ЧИСЛО УЗЛОВ НЕ МЕНЯТЬ!',4)

[t,y,BUTTON]=GINPUT(N);                 %Графический ввод предполагаемых узлов

% Обработка новых узлов

y(1)=0;                                                   % Для исключения неточности в начале

s=y;

sh=1;

%Новые значения абсцисс и ординат с учетом k

t1(1)=t(1);

for j=2:N

   sti0=(i0(j)-i0(j-1))/k;                          % Шаг по абсциссе характеристики

   for jv=1:k

      sh=sh+1;

      t1(sh)=t1(sh-1)+sti0;

   end

end

yn=spline(t,s,t1);

spl=csaps(t,s,1);

% 1-я производная в узлах

dsdt=fnder(spl);

dsdt=fnval(dsdt,t1);

% 2-я производная в узлах

ds=fnder(spl,2);

ds2=fnval(ds,t1);

dsdt2=csaps(t1,ds2,1);

figure(1);

hold on;

plot(t,s,'+',t1,fnval(spl,t1),'-.');

hold off;

grid on;

%TITLE ('Волнистость-свидетельство плохого подбора. Продолжить -"Enter"')

XLABEL('Значение "Х"')

YLABEL('Значение "У"')

LEGEND('- базовые узлы','- исходная характеристика.', ...

   '- скорректированные узлы и','- характеристика.', ...

   'Для продолжения нажать "Enter".',4)

pause

figure(2);

hold on;

plot(t1,dsdt,'-.');

hold off;

grid on;

%TITLE('В пределе кривая должна быть монотонной !')

XLABEL('Значение "Х"')

YLABEL('1-я производная от "У"')

LEGEND('- исходная и','- cкоректированная','  характеристики.', ...

'КРИВАЯ ДОЛЖНА ', 'БЫТЬ УНИМОДАЛЬНОЙ!',1)

figure(3);

hold on;

plot(t1,fnval(dsdt2,t1),'-.');

hold off;

grid on;

%TITLE('Смена знака кривой-свидетельство плохого выбор узлов')

XLABEL('Значение "Х"')

YLABEL('2-я производная от "У"')

LEGEND('- исходная и','- cкоректированная','  характеристики.',...

 'КРИВАЯ НЕ ДОЛЖНА ', ' БЫТЬ КОЛЕБАТЕЛЬНОЙ!',4)

% Сохранение на случай, если обработка откладывается «на потом»

save YY.DAT yn /ascii

save YY1.DAT dsdt /ascii

save Y.DAT y /ascii

save X.DAT t /ascii