なぜそれが必要ですか
高調波線形化法は、非線形システムの解析に広く使用されています[1]。 この方法は、2次以上のシステムで自励発振が発生する条件を決定するために使用されます。 高調波線形化では、次の2つの条件を満たす必要があります。 閉じた線形システムは、線形と非線形の2つの部分で構成する必要があります。 線形部品は、高調波に対して良好なフィルタリング特性を備えている必要があります[2]。 自動制御および調整システムには、非線形要素を含むアクチュエータが含まれているため、その分析は非常に緊急の問題です。
基本的なアルゴリズム
高調波信号を非線形要素の入力に到達させます
どこで 
出力信号を第1高調波のみのフーリエ級数に変換します。 取得するもの： 
（1）、（2）を考慮すると、非線形要素の一般方程式は[3]の形式で表すことができます。
ここで、実際の「a」は、高調波の後続のフィルタリングの条件から取得した高調波線形化係数の虚数成分です。 ─循環周波数; -微分演算子。
問題の声明
図1.2に示す2つの非線形要素を比較して、振幅の安定性を調べます。
図1「飽和」タイプの非線形性
（1）を使用して、これまでのところ実部についてのみ、以下を取得します。
図2「オープン」タイプの非線形性
同様の方法で（1）を使用すると、これまでのところ実部についてのみ、
Python実装（4）、（5）
#!/usr/bin/python # -*- coding: utf-8 -*- import matplotlib.pyplot as plt # import matplotlib as mpl # mpl.rcParams['font.family'] = 'fantasy'# mpl.rcParams['font.fantasy'] = 'Comic Sans MS, Aria'# import scipy.integrate as spint # import math # import numpy as np# q=[]# - p=[]# - c=[]# k=1# for m in np.arange(1.0,4,0.1): # def f(x): # return k*(math.sin(x))**2 a = 0# b =math.asin(1/m)# z=spint.quad(f,a,b)[0]# a =math.pi-math.asin(1/m)# b =math.pi# z1=spint.quad(f,a,b)[0]# q.append(round(z+z1,3))# - def f1(x):# return (k/m)*(math.sin(x)) a =math.asin(1/m)# b =math.pi-math.asin(1/m)# z2=spint.quad(f1,a,b)[0] # p.append(round(z+z1+z2,3))# - x=np.arange(1.0,4,0.1)# x2=np.arange(0.0,1.1,0.1)# for i in np.arange(0.0,1.1,0.1):# c.append(p[0]) plt.title(' ', size=14)# plt.xlabel(' ', size=14)# plt.ylabel(' ', size=14)# plt.plot(x, q, color='r', linewidth=2, label='-- ')# plt.plot(x, p, color='b',linewidth=2, label=' ')# plt.plot(x2, c, color='g', linewidth=2, label=' ')# plt.legend(loc='best')# plt.grid(True)# plt.show()#
結果（jpegにはありますが、グラフです！）
ところで、インポートscipy.integrateとの統合はspintとして安定しており、エラーがありません。
グラフから明らかなように、図2に示す非線形要素は、高次の高調波がフィルタリングされると、図1に示すものよりも大きな振幅安定性を提供します。 より詳細な分析では、高調波の割合が小さくなりますが、これは別の記事のトピックです。
もちろん、「他の非線形性についてはどうでしょうか、それはもっと良いのでしょうか？」 図1.2に示す典型的な特性を持つ低周波（最大5 kHz）の専用アンプに対してこれを行いました。 ちなみに、図2の特性はどの参考文献にも見当たりません。
プログラムコードを削減するオプション
たとえば、次のような記号計算を使用できます。
from sympy import * import matplotlib.pyplot as plt import numpy as np x=Symbol('x') m=Symbol('m') z=integrate(sin(x)**2, (x,0,asinh(m)))+integrate(sin(x)**2, (x,pi-asinh(m),pi)) # - z1=integrate(sin(x)**2, (x,0,asinh(m)))+integrate(m*sin(x), (x,asinh(m),pi-asinh(m)) )+integrate(sin(x)**2, (x,pi-asinh(m),pi))# print(z) print(z1)
-sin（asinh（m））* cos（asinh（m））+ asinh（m）
2 * m * cos（asinh（m））-sin（asinh（m））* cos（asinh（m））+ asinh（m）
ソフトウェアの単純化（simplify（z））は何も生成しないため、次の関係を使用します。
私達は得る： 
振幅と制限しきい値の比率を導入し、最終的に次のようになります。
関係（6）は方程式（5）を解き、システム（7）は方程式（4）を解きます。
結論
この記事では、高調波線形化の方法による非線形要素の比較分析の方法について説明します。
参照資料
- ボロノフA.A. 自動制御の理論の基礎。 特別な
線形および非線形システム。 -M .: Energoizdat、1981年、300ページ。 - Goldfarb L.S.、Baltrushevich A.V.、Netushil A.V. 理論
自動制御。 -M .: 1976年、高等学校、430秒 - クリパロE.N. 非線形自動制御システム/ E.N. 噴出。 -L。：エネルギー、1967 .-- 343 p。