從另一個角度看拉普拉斯變換

08-25

從另一個角度看拉普拉斯變換

來自專欄小潘工程師筆記721 人贊了文章

拉普拉斯變換可以說是現代工程學使用最廣泛的數學工具，它通過數學變換將微積分方程轉化成代數方程，為求解連續空間連續時間的方程提供了可能。但是，一般的教材一上來就是拉普拉斯變換的數學定義，對於其歷史和代表的深刻含義沒有任何介紹，導致很多人一直頭疼不已。今天，我們嘗試一下從不同的角度來看看拉普拉斯變換到底是怎麼回事？先從一個人說起：奧列弗.赫維賽德，一位在科學史地位被嚴重低估的人。

一、奧列弗. 赫維賽德是何許人也

奧列弗. 赫維賽德（Oliver Heaviside）是維多利亞時期英國人，出身於極度貧窮的家庭，聽力部分殘疾，還得過猩紅熱，從未上過大學，完全靠自學和興趣掌握了高等科學和數學。

圖片來自：https://en.wikipedia.org/wiki/Oliver_Heaviside

很多人熟悉赫維賽德是因為MATLAB有一個赫維賽德（Heaviside）函數，它大概長這個樣子，可以看成一個階躍函數，這個函數因為和狄拉克（Dirac）函數之間的千絲萬縷的關係而顯得尤為重要。

我們現在說赫維賽德，當然不是因為這個函數，而是因為我們嘆為觀止、驚為天人，怎麼誇都不過分的麥克斯韋方程組，麥克斯韋本人並沒有見過這個方程組，它在一定程度上應該叫「赫維賽德方程組」。

這四個公式簡直太對稱了！而且它們的含義也很清晰：變化的電場產生磁場，變化的磁場產生電場，電磁波也就是電場和磁場此消彼長、相互轉化、向前傳播的形式。多虧了赫維賽德，麥克斯韋的理論才得以在十九世紀結束之前就真正站穩了腳跟並發揚光大。可以毫不謙虛的說：宇宙間任何的電磁現象，皆可由麥克斯韋方程組解釋，包括光。也正是因為這個方程組完美統一了整個電磁場，讓愛因斯坦始終想要以同樣的方式統一引力場，並將宏觀與微觀的兩種力放在同一組式子中：即著名的「大一統理論」，不幸的是麥克斯韋成功了，而愛因斯坦沒有。

麥克斯韋早在1873年便出版了跨時代巨著《電磁通論》，可惜的是，他英年早逝，他的方程組在生前並沒有得到科學界的關注，其中一個很重要的原因是他的理論描述複雜得令人吃驚，他最初提出的電磁理論公式包含了二十個方程，直接導致了他的理論在首次發表後的10多年時間內，幾乎無人問津。

赫維賽德最偉大的貢獻是簡化了麥克斯韋的原始方程組，通過他天才般的洞察力，挖掘出了蘊含在麥克斯韋方程內部的深刻意義，從而使簡化後麥克斯韋方程組呈現出無與倫比的對稱性，成為歷史上是最漂亮的方程式（沒有唯一）。

而我們今天要說的，是赫維賽德的第二個重要貢獻：運算微積分。學過電磁學的人都知道，在歷史上人們發現好多定理公式，都是用微積分的形式表達的。1880年-1887年之間，赫維賽德在從事電磁場研究的同時，為求解微積分方程，在他的分析計算中引入了微分運算元的概念，這個方法牛逼在什麼地方呢？——它可以將常微分方程轉換為普通代數方程。天才與普通人的區別就是人家是靠「直覺」來解決問題的。赫維賽德是怎麼解微分方程的呢？他把微分、積分運算用一個簡單的運算元來代替。

——微分運算元用 $p$ 來表示，

$pleftrightarrowfrac{d()}{dx}$

——高階微分運算元為，

$p^nleftrightarrowfrac{d^n()}{dt^n}$

——積分運算元符號用 $frac{1}{p}$ 表示，

$frac{1}{p}()leftrightarrowint_{-infty}^{t}()d au$

也就是說，在某種運算元下，積分和微分對應的是倒數關係，至於運算元 $p$ 代表什麼，赫維賽德也沒有多解釋，在缺乏嚴密數學基礎的情況下，人家直接放在文章就用了，還發表了。比如常見的一個二階常微分方程，

$mfrac{d^2x}{dt}+cfrac{dx}{dt}+kx=0$

用赫維賽德的微分運算元變換一下，就變成了代數表達式，

$mp^2+cp+k=0$

赫維賽德之所以這麼做，是因為他的「物理直覺」告訴他這麼做，就是這麼硬。這顯然是一種開外掛的行為，因此也受到當時的主流數學家們們的攻訐，他們認為赫維賽德就是十足的「民科」，文章沒什麼理論依據，自己在那空想呢。當然，赫維賽德也不是弱雞，科學家懟起人來，也是毫不含糊：「因為我不能理解消化過程就拒絕晚餐嗎？不，只要我滿意這個結果。」

好了，扯了那麼遠，有童鞋已經不耐心了：這些和拉普拉斯變換有什麼關係？謎底就是：赫維賽德的微積分運算元，就是拉普拉斯變換的前身。

赫維賽德的運算元驗算雖然缺乏嚴密的數學基礎，往往能給出重要且正確的結果，方法確實有效，無法駁倒。於是在世紀之交，數學家們開始嘗試隨運算元理論進行嚴格化。後來，人在在70年前法國數學家拉普拉斯的一本有關概率論的著作上，找到了這種演算法的依據，但是這本書上提出的並不是現在我們看到的拉普拉斯變換，而是著名的Z變換。什麼？拉普拉斯變換不是拉普拉斯提出的？隨著二戰後拉普拉斯變換的廣泛使用，赫維賽德運算元的作用被弱化了，但是不可否認的是，這是這種「不正規」，僅靠「天才的直覺」而發明的方法，促成了現在拉普拉斯分析法。

二、傅里葉變換（輕量版拉普拉斯變換）

在說拉普拉斯變換以前，我們要先提一下傅里葉變換，這可以看成是輕量版的拉普拉斯變換。傅里葉變換說的是什麼事？說的是自然界的很多現象，都可以用三角函數進行分解。

圖片來源：www.bolvan.ph.utexas.edu

古巴比倫科學家在很早就用了三角函數和逼近的方法，對天體運動進行觀測和預報，1748年，大神歐拉用類似的方法分析可弦的振動（大神就是高產，哪哪都有他的身影），1753年伯努利提出任意物理弦的振動都可以可以表達為三角函數的和，但是他沒給出證明（注意，伯努利是一個家族，3代人中產生了8位科學家，後裔有不少於120位被人們系統地追溯過）。1807年，傅里葉於年在法國科學學會上提交了一篇論文，運用正弦曲線來描述溫度分布，論文里有個在當時具有爭議性的觀點：任何連續周期信號可以由一組適當的正弦曲線組合而成。傅里葉沒有做出嚴格的數學論證。這篇論文的審稿人中，有歷史上著名的數學家拉格朗日和拉普拉斯，當拉普拉斯和其它審稿人投票通過並要發表這篇論文時，拉格朗日堅決反對，認為傅里葉的方法無法表示帶有稜角的信號。法國科學學會屈服於拉格朗日的威望，拒絕了傅里葉的工作。1822年，傅里葉變換隨其著作《熱的解析》發表，但已經是15年之後了。1829年，狄利赫里通過推導其適用範圍，完善了傅里葉變換。

我們知道，三角函數可以通過歐拉公式與復指數函數聯繫起來，

$e^{iomega t}=cos(omega t)+isin(omega t)$

也就是說，復指數函數是與正弦函數、餘弦函數緊密相關的，那是不是某些信號也能通過復指數函數進行分解呢？我們不妨先來研究一下 $e^{iomega t}$ 有什麼特殊的性質。在

galieluo：自然常數「e」，工程中的自然數「1」?

zhuanlan.zhihu.comgalieluo：傅里葉變換後面的到底有什麼小秘密?

zhuanlan.zhihu.comgalieluo：被眾人膜拜的歐拉恆等式是個什麼東東？?

zhuanlan.zhihu.com

三篇文章中，我們得出了以下結論：自然常數 $e$ 可以描述事物的連續變化，虛數單位 $i=sqrt{-1}=e^{ifrac{pi}{2} }$ 可以描述旋轉，復指數函數 $e^{iomega t}$ 代表一個連續旋轉的圓。（這一點不太理解的，可閱讀上面三篇文章，裡面有詳細介紹）。

當頻率為 $omega =1rad /s$ 時，

當 $omega=2rad/s$ 時，

clc;clear;h = animatedline;xl=xlabel(cos(omegat));% yl=ylabel(sin(omegat));% grid on;title(omega = 1rad/s Made by J Pan)axis([-1,1,-1,1]);axis square;N = 100;t=linspace(0,2*pi,N);w=1;x=cos(w*t);y=sin(w*t);a = tic; % start timerfor k = 1:N addpoints(h,x(k),y(k)); hold on quiver(0,0,x(k)*1.1,y(k)*1.1) b = toc(a); % check timer if b > (1/90) drawnow % update screen every 1/30 seconds a = tic; % reset timer after updating endend

因此， ${e^{iomega t}}$ 可以看成是一組連續旋轉的曲線，幅值永遠為1，頻率一直在變。這些旋轉的曲線能幹什麼？

圖片來源網路

你能想像到很多曲線，都可以用這些不同頻率，連續旋轉的圓，通過線性疊加得到，而傅里葉定律，就是對這個結論的數學描述，傅里葉定律說：只要一個函數滿足如狄利赫里條件，都能分解為復指數函數之和，哪怕是如拉格朗日提到的帶有稜角的方波函數。狄利赫里條件為：

(1) 函數在任意有限區間內連續，或只有有限個第一類間斷點；

(2)在一個周期內，函數有有限個極大值或極小值；

(3) $x(t)$ 在絕對可積，即 $int_{0}^{infty}left| x(t) ight|precinfty$ 。

當滿足狄利赫里條件時，傅里葉變換及其逆變換為：

$F(omega)=int_{-infty}^{+infty}f(t)e^{-iomega t}dt$

$f(t)=frac{1}{2pi}int_{-infty}^{+infty}F(omega)e^{iomega t}domega$

於是就可以很好的解釋拉格朗日和傅里葉之間的爭論了——拉格朗日是對的：正弦曲線無法組合成一個帶有稜角的信號，稜角處會有很小高頻波動（吉布斯現象）。但是，我們可以用正弦曲線來非常逼近地表示它，逼近到兩種表示方法不存在能量差別，基於此，傅里葉也是對的。一個從數學家的角度，一個從工程師的角度。

三、拉普拉斯變換（原來就是那麼回事）

傅里葉變換能幫我們解決很多問題，一經問世後便受到廣大工程師們的喜愛，因為它給人們提供了一扇不同的窗戶來觀察世界，從這個窗戶來看，很多事情往往變得簡單多了。但是，別忘了，傅里葉變換有一個很大局限性，那就是信號必須滿足狄利赫里條件才行，特別是那個絕對可積的條件，一下子就攔截掉了一大批函數。比如函數 $f(t)=t^2$ 就無法進行傅里葉變換。這點難度當然拿不到聰明的數學家們，他們想到了一個絕佳的主意：把不滿足絕對的可積的函數乘以一個快速衰減的函數，這樣在趨於 $infty$ 時原函數也衰減到零了，從而滿足絕對可積。

圖片來自：https://www.zhihu.com/question/22102732/answer/369089156

數學描述是：

$lim_{x ightarrow +infty}f(x)e^{-sigma x}=0,sigmain R$

為保證 $e^{-sigma x}$ 一直為衰減函數，我們把 $x$ 定義域縮減到正半軸，這樣可以進行傅里葉變換就變成了：

$F(omega)=int_{0}^{+infty}f(t)e^{-sigma t}e^{-iomega t}dt=int_{0}^{+infty}f(t)e^{-(sigma+iomega) t}dt$

如果假設：

$s=sigma+iomega$

那麼就得到：

$F(s)=int_{0}^{+infty}f(t)e^{-s t}dt$

這就是拉普拉斯變換。

我知道大部分人前面的數學推導沒什麼興趣，接下來就是放彩蛋的時刻了，很多童鞋會說不管傅里葉變換或者拉普拉斯變換是什麼細節，你能說點有意思的，讓人能記憶深刻的信息嗎？

我前面的專欄文章說了， $e^{i omega t}$ 代表按不同頻率旋轉的單位圓，那是在複平面來看的，想像力豐富的同學可以腦補一下，如果把時間軸也加上， $e^{i omega t}$ 長什麼樣子呢？那就是螺旋曲線！

clc;clear;h = animatedline;h1=gcf;view(3);xl=xlabel(cos(omegat));% yl=ylabel(sin(omegat));% zl=zlabel(t);% set(xl,Rotation,30);% set(yl,Rotation,-30);%grid on;title(omega = 1rad/s Made by J Pan)axis([-1,1,-1,1,0,4*pi])N = 200;t=linspace(0,4*pi,N);w=1;x=cos(w*t);y=sin(w*t);a = tic; % start timerfor k = 1:N addpoints(h,x(k),y(k),t(k)); hold on line([0 x(k)],[0 y(k)],[t(k) t(k)],Color,red) b = toc(a); % check timer if b > (1/90) drawnow % update screen every 1/30 seconds a = tic; % reset timer after updating endend

從其他不同的平面看，螺旋曲線就是正弦或者餘弦曲線。

下面我們再來看看 $e^{-st},(s=sigma+iomega)$ 長什麼樣子：

clc;clear;h = animatedline;h1=gcf;view(3);xl=xlabel(cos(omegat));% yl=ylabel(sin(omegat));% zl=zlabel(t);% set(xl,Rotation,30);% set(yl,Rotation,-30);%grid on;title(omega = 1rad/s Made by J Pan)axis([-1,1,-1,1,0,4*pi])N = 200;t=linspace(0,4*pi,N);w=1;sig=-0.2;x=exp(sig*t).*cos(w*t);y=exp(sig*t).*sin(w*t);a = tic; % start timerfor k = 1:N addpoints(h,x(k),y(k),t(k)); hold on line([0 x(k)],[0 y(k)],[t(k) t(k)],Color,red) b = toc(a); % check timer if b > (1/90) drawnow % update screen every 1/30 seconds a = tic; % reset timer after updating endend