為何LC振蕩電路的物理模擬與傳函模擬結果不符？

12-29

在simulink裡面模擬。一組LC濾波電路，電感 l=1.014e-2H,電容c=1e-3F，使得其剛好在工頻電壓處發生諧振。模擬電路圖如下：

電壓源為1v，50hz的工頻電壓。
波形圖如下：

可以看到發生了諧振，電容兩端的電壓被放大到一千多。

如果把它化成控制框圖的形式，並進行模擬，如下所示：

波形如下;

同樣，信號源採用幅值為1，50hz的正玄信號，這個控制框圖應該是與上面的電路圖等效的，同樣也是發生諧振了，但是為何波形是不斷放大的呢？
更為有趣的是，我把模擬時間設置到100s後，發現它還有一個放大縮小的周期。

結論：

因為LC振蕩電路的固有頻率與輸入頻率有微小的偏差，所以系統的一般解會是兩個頻率的正弦波的疊加。
由於頻率的偏差很微小，所以發生了"拍頻"現象。Beat (acoustics)
兩種模擬結果不同可能是由於模擬系統的初值不同所導致的，而由於系統臨界穩定，不同的初值不會收斂到同一個穩態振蕩解。

簡而言之，我猜測電路模擬是從」穩態「所對應的初值開始計算的；而傳函模擬是從零初值開始計算的。由於名義系統是臨界穩定的，所以後者的模擬到達不了「穩態」，這就是為什麼傳函模擬得不到穩定的正弦波。

利用伯德圖的信息或者系統的解析解，其實可以手動調節系統的初值使得系統直接從穩態開始計算，不過就不能用傳函模塊而得轉用狀態空間模塊了。

以下是詳細的系統分析，嫌麻煩的話可以跳過直接看模擬結果。設 $y=v_c$ ，那麼系統方程是

$LCddot{y}+y=sin(100pi t)$

系統的固有頻率是 $omega_n=frac{1}{sqrt{LC}}$ ，輸入正弦波的頻率是 $omega=100pi$ 。倘若 $omega_n=omega$ ，那麼系統的解會趨於無窮大。然而根據題中所給的數據，兩者還是有很明顯的誤差，前者是314.1593，後者是314.0371。這導致了系統存在一個有限的穩態（受迫）響應。此處也可以通過伯德圖來證實：

s=tf("s"); L=1.014e-2 C=1e-3 [A,phi] = bode(1/(L*C*s^2+1),50*2*pi)

即，假設輸入為 $u(t)=sin(omega t)$ ，存在一個穩態輸出 $y(t)=Asin(omega t+phi)$ ，其中 $A=1285.5$ ， $phi=-pi$ 。此處和 @李崇提供的解析解 $y(t)=frac{sin(omega t)}{LC(omega_n^2-omega^2)}$ 是吻合的。畫出此時的解的圖形：

可以看出這個解對應著電路模擬的結果。計算此解對應的系統初值： $y(0)=0,dot{y}(0)=-4.0386e+05$ 。值得注意的一點是，此處所需的系統初值並不是零初值。另一方面，由於系統是臨界穩定的，存在一個穩態解不代表系統會從任意初值趨近於該穩態解。因此，倘若系統初值並非此特殊初值，系統的響應還會疊加上系統自身的模態。比如說，假如系統從零初值開始，那麼系統的解析解可以用MATLAB symbolic toolbox求出（感謝匿名用戶和 @Falccm 指出解析解的數值問題）。

L=1.014e-2 C=1e-3 syms y(t) eqn = diff(y,t,2) == 1/(L*C)*(-y+sin(100*pi*t)) Dy = diff(y,t); cond = [y(0)==0, Dy(0)==0]; ySol(t) = simplify(dsolve(eqn,cond)) vpa(ySol(t),5)

得出 $y(t)$ 的表達式為：

yt=@(t) 1286.0*sin(314.04*t) - 1285.5*sin(314.16*t)

可以看出，此時系統的解裡面，除了輸入頻率 $omega$ 所對應的模態（角頻率314.16），還有系統自身模態 $omega_n$ （314.04）。由於 $omega$ 和 $omega_n$ 很接近所導致的拍頻現象Beat (acoustics)，系統的輸出看起來會有一個 $frac{omega-omega_n}{2}=0.0611$ 的頻率，對應著102.9 s的周期，即問題描述中正弦波「放大縮小」兩次的周期。將此情況下的系統軌跡畫圖可以看出與題中傳函的模擬結果一致：

plot(0:0.0001:100,yt(0:0.0001:100))

可以看出此處的解析解符合題主傳函模擬的結果。

為了進一步驗證，寫了一個小模擬程序。結果是這樣的：

（鏈接: https://pan.baidu.com/s/1pLarN1h 密碼: hdvu）

假如使用 $y(0)=0,dot{y}(0)=0$ 初值進行計算：

假如使用 $y(0)=0,dot{y}(0)=-4.0386e+05$ 計算：

可以看出模擬的結果精確地符合理論結果（前面給出的解析解）。

結論：是兩種模擬的系統初值不同所導致的，傳函的模擬是從零初值開始計算的，而對於電路模擬，我猜測它是直接從穩態開始計算的。

關於臨界穩定和不穩定系統的「穩態」，我之前也寫過一個答案：

不穩定系統的頻率響應（伯德圖）的物理意義是什麼？ - Kaixiang Wang 的回答

另外略微反對一下 @李崇的答案。

從以上的分析中可以看出，模擬結果精確吻合系統解析解，因此並不是數值計算精度的問題。
系統初值為零時的解並不是穩態解；穩態解需求一個非常特殊的系統初值。
我之前也懷疑電路模擬的時候simulink自動加入了電阻，不過這樣必然會改變系統的頻率響應，而當前的頻率響應精確地反應了模擬中輸出的幅值，所以我感覺這一點可能也不對。

大家的討論非常精彩。剛剛看到，補充幾點吧：

對於線性時不變系統，輸出的模態是由輸入的模態和系統的模態決定的。也就是說，輸出中不會無緣無故地產生輸入和系統中不存在的模態。

確切的說，輸出的模態是輸入的模態和系統的模態的線性疊加。

模態疊加的權重是由初值以及系統的參數決定的。某個模態的權重可以為零，也可以為無窮大。

模態其實是極點。輸入的模態可以由系統的零點抵消，同樣的，系統的模態也可以由輸入的零點抵消。

希望有人能做出一些典型的兩個正弦信號相疊加的形狀，乃至三個正弦信號相疊加等等。這樣的話，有助於大家辨認一些信號。比如第二個圖中的好像很奇怪的信號，放大之後就看的更清楚了。其實應該就是兩個不同頻率的正弦信號的線性疊加。至於解析解中出現了的正弦的乘法，可以用積化和差公式繼續化簡，可能就看得更清楚了。

---

對信號求過傅立葉級數，看各個分量就更清楚了。

有似乎不很重要的一點需要指出：如果我沒算錯的話，MATLAB symbolic toolbox恐怕要背個鍋。

不知道 @Kaixiang Wang的這段話有沒有童鞋仔細看：

可以看出，此時系統的解裡面，除了 $omega$ 所對應的模態（角頻率628.2），還有系統自身模態 $omega_n$ （314.0）以及一個奇怪的模態（0.1221）。

那個角頻率628.2的東西究竟是什麼？這可不是100pi，這分明是200pi好伐！？

而且，如果按照 @李崇給出的那個通式，結合 @Kaixiang Wang的初始條件兩個0，會發現最終結果其實只有兩項就夠了，形式應該是K[sin(w*t)-sin(100pi*t)]，算一下會發現K是1285.53，w是314.03715左右，100pi是314.15927左右，後者減去前者就會得到那個傳奇的0.1221。

重寫 @Kaixiang Wang算出來的這個結果：

yt=@(t)1286.0*sin(314.0*t) - 1.0*sin(314.0*t).*(0.25*cos(628.2*t) + 1286.0*cos(0.1221*t)) + 1.0*cos(314.0*t).*(0.25*sin(628.2*t) - 1286.0*sin(0.1221*t))

用前面設定的符號，寫得清楚一點應該是這樣的：

$Ksin(omega t)-sin(omega t)[frac{1}{4}cos(200 pi t)+Kcos(100 pi t -omega t) ]+cos(omega t)[frac{1}{4}sin(200 pi t)-Ksin(100 pi t -omega t) ]$

好吧我們拆解一下，會發現兩組既不是積化和差也不是和差化積，這叫三角函數加法公式。

（好吧我知道「加法公式」這個名字很俗，其實我第一反應也會管這個叫積化和差或者和差化積，然而誰叫我後說話，來打我啊）

這個事是這個樣子滴，我們再挑戰一下知乎的公式編輯器好了：

原式=

$=Ksin(omega t)-K[cos(omega t)sin(100 pi t-omega t)+sin(omega t)cos(100 pi t-omega t)-frac{1}{4}[sin(omega t)cos(200 pi t)-cos(omega t)sin(200 pi t)]$

$=Ksin(omega t)-Ksin(100 pi t - omega t+omega t)-frac{1}{4}sin(omega t -200 pi t)$

到這裡，前兩項我們喜聞樂見，最後一項是個什麼鬼？注意相比較K=1285.53來說，1/4很小很小，而且這玩意頻率雖然比前兩個大但也差不多（大約是314.28138），因此在圖上是看不出來滴。

這裡可能要找 @小姜，請問這個莫名其妙的模態是幾個意思啊……

同時，解釋得很好的0.1221，其實可以由前兩項搞一個和差化積得出來，應該是這個樣子滴：

$=2Kcos(frac{omega t + 100 pi t}{2})sin(frac{omega t - 100 pi t}{2})=2 imes 1285.53 cos(314.0982t)sin(-0.06106t)$

2*1285.53=2571.06才是幅值的本尊，而0.06106對應周期大概102.9s才是周期的本尊。

收工，不熟悉知乎敲公式真浪費時間啊……

本文已加入matlab工具箱內幕+三角函數基本公式名稱複習豪華套餐系列。

確實和系統的初值關聯很大，由此引申出來的控制理論問題也值得仔細討論和深究。多多關注 @李崇@Kaixiang Wang 他們關於初值和是否發散的分析，Kaixiang Wang的分析很有道理，完全印證了結果的由來。

鑒於最近較忙，理論分析就沒時間深究了，如果只是單純的答題主的這個具體問題，貼兩張圖就好了。

記得勾選設置零初值

圖沒毛病呀

這是solver和absolute tolerance的問題。試試ode23，並且把absolute tolerance改小。

這個結果確實和理論不符，主要是數值計算模擬的精度問題。

我前一陣子遇到過這個問題，多謝 @小姜大神給予的幫助

----------------------------------------------------------

經@Kaixiang Wang提醒了做了模擬

系統圖

ode45，variable step, absolute tolerance為auto

ode45，variable step, absolute tolerance為1e-5

ode23，variable step, absolute tolerance為auto

ode23，variable step, absolute tolerance為1e-5

這個基本上就是數值計算精度的鍋。這是臨界穩定系統，共振頻率上放大幅度大，absolute tolerance設為auto的或者小的話，ode solver的精度不夠導致最後結果不能收斂，形不成穩態

------------------------------------------------------------

@Kaixiang Wang

其實這就是個二階無阻尼系統的震蕩問題，可以描述為

$addot{x}+bx=Fsin(omega t)$

其中系統自身的諧振頻率為

$omega_n=sqrt{b/a}$

系統的解為

$x(t)=frac{Fsin(omega t)}{a(omega_n^2-omega^2)}$

由此可見，當輸入頻率等於諧振頻率的時候，理論上來說分母會等於於0，會把輸入頻率給放大無窮倍，也就是 @Kaixiang Wang說的穩定模態不存在的問題。

但是計算機在做浮點數運算的時候，首先會拒絕做除0的問題，會指定兩個浮點數的差小到一定閾值的時候就認為兩個數相等了。

比如說在c語言當中

double a,b;

.............

if(a==b)

{

..............

}

這種寫法是不可以的，因為這個表達式絕難成立。

而是一般會用

#define ERROR 1e-7

double a,b;

.............

if(abs(a-b)&{

..............

}

這種的寫法。而在matlab這種精度高的數學計算軟體裡面，這個閾值還是可能會變的。

回歸到這個

$addot{x}+bx=Fsin(omega t)$

問題裡面，理論上當輸入頻率嚴格吻合了諧振頻率的時候，輸出是要無界的。但是計算機數值精度導致：

（1）避免除0

（2） $(omega_n^2-omega^2)$ 難以等於0.

我懷疑我的結果和別人不符，是因為我輸入信號的頻率是手動輸入了小數點後四位，而別人用的是自動計算生成的具有小數點後幾十位精度的。

所以說來之前的答案應該都是有問題的。

-------------------------------------------------------------------

再想說一下電路模擬收斂到穩態的問題。我沒用過simulink這個模塊，不知道是否電路裡面增加了等效的電阻，如果有電阻，系統就會變成

$ddot{x}+frac{omega_n}{Q} dot{x}+omega_n^2 x=Fsin(omega t)/a$

的形式。Q是電路的品質因數，在振蕩器電路裡面都是希望Q越高越好，Q越高對於輸入放大的倍數越大，能量損耗越低。在現實中理想LC震蕩電路是不存在的，都是會有電阻的

另外關於初值的問題。當一個系統極點都在虛軸上時候，如果系統的初始值不為0，那麼它就會以自身的諧振頻率成正弦形式振蕩起來。電路裡面經常會用這樣的閉環反饋結構來生成正弦波，控制迴路裡面會有一個automatic gain controller來設定其初始值，並在波形衰減的時候再提供stimulus。關於初始值不為0，又有外部正弦輸入的情況下會是怎麼樣，我沒考慮過，不清楚

----------------------------------------------------------------------------

另外再說下simulink計算精度的問題。還是用

$ddot{x}+frac{omega_n}{Q} dot{x}+omega_n^2 x=Fsin(omega t)/a$

這個系統，當Q小的時候，理論上來說諧振頻率上的輸入是一個有界的90度滯後的正弦波，當Q開始增大的時候，數值精度不夠的問題就體現出來了，正弦波的幅值會開始變化，這是因為數值精度不夠陷入limit cycle了。

-------------------------------------------------------------------------------

最新總結

又去看了機械振動的書，發現前面又很多地方說的不對。總結一下：

考慮一個線性系統