p = tf('p') % декларируем р как оператор Лапласа
Wlc=tf(2/(0.05*p+1)/(0.02*p+1)/p) % вводим передаточную функцию линейной части системы
nyquist(Wlc) % строим АФЧХ линейной части системы
syms A pi % декларируем переменные коэффициентов гармонической линеаризации
Wne=440*sqrt(A^2-0.25^2)/(pi*A^2) + j*0 % вводим частотную гармоническую ПФ нелинейной части
simplify(-1/Wne) % преобразуем Wne в обратную амплитудно-фазовую характеристику нелинейности с обратным знаком –1/WНЭ(jА)
pretty(-1/Wne) % выводим –1/WНЭ(jА) в дробном виде, выделяем вещественную и мнимую части
UWne1=-(pi*A^2)/(440*(A^2 - 1/16)^(1/2)) % вводим вещественую часть –1/WНЭ(jА)
VWne1=0 % вводим мнимую часть –1/WНЭ(jА)
dataA=[0.25:0.025:10]' % расчитываем значения аргумента коэффициентов гармонической линеаризации
dataUWne1=subs(UWne1,A,dataA) % вычисляем значения вещественой части –1/WНЭ(jА)
dataVWne1=dataUWne1*0 % вычисляем значения мнимой части –1/WНЭ(jА)
% далее вручную строим на имеющемся плоте а АФЧХ линейной части системы
% график –1/WНЭ(jА) по значениям вещественой и мнимой частей –1/WНЭ(jА) командой Add Data
% ищем точку пересечения гафиков Wlc и –1/WНЭ(jА) и значение функции –1/WНЭ(jА) (в примере -0,02844)
solve(UWne1+0.02844) % находим знаения аргумента –1/WНЭ(jА) при –1/WНЭ(jА)= -0,02844 или –1/WНЭ(jА)+0,02844 = 0
pi=3.1416 % вводим значение pi
% еще раз считаем уже с pi = 3,1416, только положительные аргументы
A1=(2*(3741529239777681/152587890625 - (61168041*pi^2)/1562500)^(1/2) + 122336082/390625)^(1/2)/(2*pi)
A2=(122336082/390625 - 2*(3741529239777681/152587890625 - (61168041*pi^2)/1562500)^(1/2))^(1/2)/(2*pi)