第三章 似曾相識
曲線有問題……這是一種近乎本能的判斷，來自吳黔過去十幾年處理理論模型和數值結果時形成的經驗。
如果問題出在初始參數，能量曲線通常不會以如此穩定的趨勢增長。眼前的結果更像是某種微小誤差，在一次次計算中不斷積累。
這就好像一道數學題，如果答案是「I am Li Hua.」那你一眼就能看出這是建模出了問題，但如果答案的每個公式都正確，但是結果卻錯了，就要考慮是不是計算出了問題。
顯然，後者解決起來更加困難。
現在就是這樣。吳黔能肯定趙啟明他們的物理建模沒有問題，但是為什么正確的建模會引入這樣大的誤差，他還不是很確定。
其實按照他多年浸淫物理的經驗，應該至少能一眼判斷出錯誤大致出現在哪裡。可現在他的大腦就像被一層迷霧所籠罩，根本沒辦法提取出有用的東西。
這種感覺讓吳黔很難受，就像一個很渴的人突然得到了一瓶冰鎮汽水，卻怎麼都擰不開瓶蓋。
放棄了無意義的空想，吳黔決定等結束之後再去找趙啟明問一問。
交流會結束後，大部分學生陸續離開。吳黔順著人流，走到正在收拾東西的趙啟明身旁。
「我對你剛才介紹的東西很感興趣，能讓我看一下你們的程序嗎？」
趙啟明抬起頭：「你懂編程？大幾的？」
「物理學院，大一八班，吳黔。編程的話，我之前學過一點。」
趙啟明點點頭，打開還沒來得及關上的電腦：「喏，程序就在這裡了，你自己看看吧。不懂的可以問我。」說罷，他繼續收拾起桌上的其他東西，顯然沒指望眼前這個大一新生能看出什麼東西來。
程序是用MATLAB寫的，雖然吳黔感覺自己已經對其細節淡忘了許多，但是靠著代碼的注釋，他還是勉強看懂了趙啟明他們使用的方法。
他們先把雙擺的二階微分方程轉化成了一階方程組，質量、擺長、重力加速度和初始條件也沒有明顯問題。
吳黔繼續向下看。
每經過一個時間步長，程序都會根據當前狀態計算加速度，再更新角速度和角度。
「這裡用的是什麼方法？」
「歐拉法。」趙啟明說道，「最開始圖方便，時間步長設置的比較大。發現問題後我們把步長從零點零一改成了零點零零一，但問題還是存在，只是能量增長得慢了一些。」
歐拉法。
聽到這個名字時，吳黔腦海里的迷霧仿佛晃動了一下。問題很可能就在這裡。
但為什麼？應該用什麼替代？他的大腦依舊拒絕工作。
「我懷疑和數值方法有關。」沉思了一會兒，吳黔還是決定把自己的猜想說出來。畢竟，這個項目看上去是解決任務的最佳選擇。
趙啟明立刻問道：「具體是什麼問題？」
吳黔沉默片刻，沒有逞強。
「現在還不能確定。能把程序發給我一份嗎？我回去檢查一下。」
「可以。不過我們明天下午要向指導老師匯報，如果你能找到原因，最好在那之前聯繫我。」
離開實驗樓後，吳黔沒有回宿舍，而是徑直走向圖書館。
「數值計算……常微分方程……」
【記住本站域名 TW 看書網解悶好，🆂🅷🆄.🆃🆆超流暢 】
他沿著書架尋找，很快抽出一本《數值分析》和一本《常微分方程數值解法》。
對於從未接觸過數值計算的學生而言，僅僅理解這些書里的基本概念，可能就要花費數周。
吳黔在閱覽室角落坐下，翻開了關於歐拉法的章節。一開始，他看得並不快，離散化、局部截斷誤差、整體誤差……這些概念他都好像在哪裡見過，可是大腦就是拒絕認出他們。
但十幾分鐘後，他感覺腦袋中的什麼限制好像被打開了，那些他覺得似曾相識的概念變得越來越熟悉。當第一個概念被重新理解，更多被迷霧籠罩的知識也隨之浮現。
歐拉法利用當前點的斜率預測下一時刻的狀態，每一步都會產生截斷誤差。對於長時間演化的問題，這些誤差可能不斷積累。雙擺對微小擾動又極其敏感，數值誤差會迅速放大，最終表現成虛假的能量增長。
還是用之前冰面小球的例子，小球在模擬中之所以會莫名其妙飛到天上，原因就是在每一刻都有一個隱形的推力作用在了它上面。一開始，這個力造成的影響並不大，但是時間拉長以後，這些隱形推力的合力就把小球扔了出去。
「這就是問題出現的關鍵！」
吳黔翻頁的速度越來越快。最開始，他幾分鐘才能讀完一頁。後來只需要幾十秒。再往後，他的目光只要在公式上略作停留，便能理解推導過程，然後迅速推導出下一步。
兩個小時後，吳黔合上第一本書。問題已經弄清楚了，但指出原因還不算解決。
他翻開第二本書，這本書涉及的是具體的常微分方程的數值解法，相比上一本來說更加符合他的需要。
改進歐拉法、龍格—庫塔方法、辛算法……
一種種常微分方程的數值解法重新回到他的腦海中。
辛算法更適合長時間保持哈密頓系統的結構，但對一個本科生項目而言，直接換成四階龍格—庫塔法更容易實現，也足夠把能量漂移壓到可接受的範圍。
如果說歐拉法只是對方程潦草的猜測一下，龍格—庫塔法就像是一個偏執狂。這個方法會按照階數對方程做多次近似，最後得到一個更加靠譜的結果。
也就是說，採用這個方法之後，誤差會大大減少。雖然不至於消失，但是相對於模擬需要的系統精度來說已經綽綽有餘。
吳黔拿出紙筆，重新推導公式，對照書中的步驟，將解法的關鍵位置寫在草稿紙上。沒過多久，看著寫的滿滿當當的草稿紙，吳黔轉了轉手上的筆……從幾乎記不起歐拉法的性質，到重新掌握常微分方程的數值求解，他只用了三個小時。
「搞定！接下來就是看怎麼把新的數值方法寫進程序了。」
吳黔來到圖書館的電子閱覽區——他的電腦上沒有裝MATLAB。
「回去得把這幾個常用的軟體下回來了。」
MATLAB、Mathematica、LaTeX，再加上一套能用的Linux環境，幾乎都是做物理研究繞不開的工具。想要重新開始科研，首先就得把過去的工作環境一點點搭建回來。
不過2014年的電腦配置遠沒有後來方便。Linux環境要麼使用虛擬機，要麼單獨安裝雙系統，光是配置這些東西，恐怕就要花費不少時間。
總不能每次處理數據，都跑到圖書館來用公共電腦。
吳黔導入下午趙啟明發給他的腳本。MATLAB的具體語法同樣有些陌生，但看過幾段示例後，熟悉感很快回來了。
他將原本的時間推進部分改成四階龍格—庫塔格式，隨後輸入原來的初始條件，按下運行。
幾秒後，兩條能量曲線同時出現在屏幕上。
原程序計算出的總能量依舊不斷上升。修改後的曲線雖然存在細微波動，卻始終穩定在一條極窄的範圍內。
吳黔把模擬時間擴大十倍，曲線依舊沒有出現原本那種明顯的能量漂移。問題找到了，解決方案也完成了。
他把修改後的程序和說明發給趙啟明。
不到一分鐘，對方回復道：「你改好了？」後面跟著一個小人臥槽表情。
「嗯。主要原因是歐拉法產生的誤差不斷積累。縮小步長只能減緩漂移，不能真正解決。」吳黔也回了一張小人跳舞表情。
幾分鐘後，趙啟明發來一張新的截圖。
「真的解決了！」
「我們把模擬時間提高了十倍，能量誤差還是比原來小了兩個數量級。」
「你明天下午有時間嗎？我要去找指導老師，你也來一趟吧。」
就在這條消息彈出的瞬間，淡藍色文字浮現在吳黔眼前。
【檢測到宿主獨立發現並解決現實學術問題。】
【解決方案已獲得相關人員認可。】
【新手任務完成。（新手任務無評級）】
【獲得積分300點。】
【獲得初級抽獎機會一次。】
【學術科技樹解鎖中……】
舒展了一下長久伏案有些酸澀的身體，吳黔長舒一口氣。沒想到這個看上去困難的任務這麼簡單就被解決了。
「寶刀未老啊！」吳黔有些得意的想道。
如果問題出在初始參數，能量曲線通常不會以如此穩定的趨勢增長。眼前的結果更像是某種微小誤差，在一次次計算中不斷積累。
這就好像一道數學題，如果答案是「I am Li Hua.」那你一眼就能看出這是建模出了問題，但如果答案的每個公式都正確，但是結果卻錯了，就要考慮是不是計算出了問題。
顯然，後者解決起來更加困難。
現在就是這樣。吳黔能肯定趙啟明他們的物理建模沒有問題，但是為什么正確的建模會引入這樣大的誤差，他還不是很確定。
其實按照他多年浸淫物理的經驗，應該至少能一眼判斷出錯誤大致出現在哪裡。可現在他的大腦就像被一層迷霧所籠罩，根本沒辦法提取出有用的東西。
這種感覺讓吳黔很難受，就像一個很渴的人突然得到了一瓶冰鎮汽水，卻怎麼都擰不開瓶蓋。
放棄了無意義的空想，吳黔決定等結束之後再去找趙啟明問一問。
交流會結束後，大部分學生陸續離開。吳黔順著人流，走到正在收拾東西的趙啟明身旁。
「我對你剛才介紹的東西很感興趣，能讓我看一下你們的程序嗎？」
趙啟明抬起頭：「你懂編程？大幾的？」
「物理學院，大一八班，吳黔。編程的話，我之前學過一點。」
趙啟明點點頭，打開還沒來得及關上的電腦：「喏，程序就在這裡了，你自己看看吧。不懂的可以問我。」說罷，他繼續收拾起桌上的其他東西，顯然沒指望眼前這個大一新生能看出什麼東西來。
程序是用MATLAB寫的，雖然吳黔感覺自己已經對其細節淡忘了許多，但是靠著代碼的注釋，他還是勉強看懂了趙啟明他們使用的方法。
他們先把雙擺的二階微分方程轉化成了一階方程組，質量、擺長、重力加速度和初始條件也沒有明顯問題。
吳黔繼續向下看。
每經過一個時間步長，程序都會根據當前狀態計算加速度，再更新角速度和角度。
「這裡用的是什麼方法？」
「歐拉法。」趙啟明說道，「最開始圖方便，時間步長設置的比較大。發現問題後我們把步長從零點零一改成了零點零零一，但問題還是存在，只是能量增長得慢了一些。」
歐拉法。
聽到這個名字時，吳黔腦海里的迷霧仿佛晃動了一下。問題很可能就在這裡。
但為什麼？應該用什麼替代？他的大腦依舊拒絕工作。
「我懷疑和數值方法有關。」沉思了一會兒，吳黔還是決定把自己的猜想說出來。畢竟，這個項目看上去是解決任務的最佳選擇。
趙啟明立刻問道：「具體是什麼問題？」
吳黔沉默片刻，沒有逞強。
「現在還不能確定。能把程序發給我一份嗎？我回去檢查一下。」
「可以。不過我們明天下午要向指導老師匯報，如果你能找到原因，最好在那之前聯繫我。」
離開實驗樓後，吳黔沒有回宿舍，而是徑直走向圖書館。
「數值計算……常微分方程……」
【記住本站域名 TW 看書網解悶好，🆂🅷🆄.🆃🆆超流暢 】
他沿著書架尋找，很快抽出一本《數值分析》和一本《常微分方程數值解法》。
對於從未接觸過數值計算的學生而言，僅僅理解這些書里的基本概念，可能就要花費數周。
吳黔在閱覽室角落坐下，翻開了關於歐拉法的章節。一開始，他看得並不快，離散化、局部截斷誤差、整體誤差……這些概念他都好像在哪裡見過，可是大腦就是拒絕認出他們。
但十幾分鐘後，他感覺腦袋中的什麼限制好像被打開了，那些他覺得似曾相識的概念變得越來越熟悉。當第一個概念被重新理解，更多被迷霧籠罩的知識也隨之浮現。
歐拉法利用當前點的斜率預測下一時刻的狀態，每一步都會產生截斷誤差。對於長時間演化的問題，這些誤差可能不斷積累。雙擺對微小擾動又極其敏感，數值誤差會迅速放大，最終表現成虛假的能量增長。
還是用之前冰面小球的例子，小球在模擬中之所以會莫名其妙飛到天上，原因就是在每一刻都有一個隱形的推力作用在了它上面。一開始，這個力造成的影響並不大，但是時間拉長以後，這些隱形推力的合力就把小球扔了出去。
「這就是問題出現的關鍵！」
吳黔翻頁的速度越來越快。最開始，他幾分鐘才能讀完一頁。後來只需要幾十秒。再往後，他的目光只要在公式上略作停留，便能理解推導過程，然後迅速推導出下一步。
兩個小時後，吳黔合上第一本書。問題已經弄清楚了，但指出原因還不算解決。
他翻開第二本書，這本書涉及的是具體的常微分方程的數值解法，相比上一本來說更加符合他的需要。
改進歐拉法、龍格—庫塔方法、辛算法……
一種種常微分方程的數值解法重新回到他的腦海中。
辛算法更適合長時間保持哈密頓系統的結構，但對一個本科生項目而言，直接換成四階龍格—庫塔法更容易實現，也足夠把能量漂移壓到可接受的範圍。
如果說歐拉法只是對方程潦草的猜測一下，龍格—庫塔法就像是一個偏執狂。這個方法會按照階數對方程做多次近似，最後得到一個更加靠譜的結果。
也就是說，採用這個方法之後，誤差會大大減少。雖然不至於消失，但是相對於模擬需要的系統精度來說已經綽綽有餘。
吳黔拿出紙筆，重新推導公式，對照書中的步驟，將解法的關鍵位置寫在草稿紙上。沒過多久，看著寫的滿滿當當的草稿紙，吳黔轉了轉手上的筆……從幾乎記不起歐拉法的性質，到重新掌握常微分方程的數值求解，他只用了三個小時。
「搞定！接下來就是看怎麼把新的數值方法寫進程序了。」
吳黔來到圖書館的電子閱覽區——他的電腦上沒有裝MATLAB。
「回去得把這幾個常用的軟體下回來了。」
MATLAB、Mathematica、LaTeX，再加上一套能用的Linux環境，幾乎都是做物理研究繞不開的工具。想要重新開始科研，首先就得把過去的工作環境一點點搭建回來。
不過2014年的電腦配置遠沒有後來方便。Linux環境要麼使用虛擬機，要麼單獨安裝雙系統，光是配置這些東西，恐怕就要花費不少時間。
總不能每次處理數據，都跑到圖書館來用公共電腦。
吳黔導入下午趙啟明發給他的腳本。MATLAB的具體語法同樣有些陌生，但看過幾段示例後，熟悉感很快回來了。
他將原本的時間推進部分改成四階龍格—庫塔格式，隨後輸入原來的初始條件，按下運行。
幾秒後，兩條能量曲線同時出現在屏幕上。
原程序計算出的總能量依舊不斷上升。修改後的曲線雖然存在細微波動，卻始終穩定在一條極窄的範圍內。
吳黔把模擬時間擴大十倍，曲線依舊沒有出現原本那種明顯的能量漂移。問題找到了，解決方案也完成了。
他把修改後的程序和說明發給趙啟明。
不到一分鐘，對方回復道：「你改好了？」後面跟著一個小人臥槽表情。
「嗯。主要原因是歐拉法產生的誤差不斷積累。縮小步長只能減緩漂移，不能真正解決。」吳黔也回了一張小人跳舞表情。
幾分鐘後，趙啟明發來一張新的截圖。
「真的解決了！」
「我們把模擬時間提高了十倍，能量誤差還是比原來小了兩個數量級。」
「你明天下午有時間嗎？我要去找指導老師，你也來一趟吧。」
就在這條消息彈出的瞬間，淡藍色文字浮現在吳黔眼前。
【檢測到宿主獨立發現並解決現實學術問題。】
【解決方案已獲得相關人員認可。】
【新手任務完成。（新手任務無評級）】
【獲得積分300點。】
【獲得初級抽獎機會一次。】
【學術科技樹解鎖中……】
舒展了一下長久伏案有些酸澀的身體，吳黔長舒一口氣。沒想到這個看上去困難的任務這麼簡單就被解決了。
「寶刀未老啊！」吳黔有些得意的想道。