Embedded Files
方程式& code 好讀版：https://www.notion.so/TLM-adjoint-code-9ab5c28237174e9a8ccfb142c7f80c32
(強烈建議 ↑ 點此閱讀)
本篇主要介紹如何根據一個現有的模式的code寫出相應的tangent linear model (TLM，切線性模式) 跟 adjoint model (ADJ，伴隨模式)，關於TLM跟ADJ的簡介請見→這篇。
模式說明
模式說明
這邊以一個非線性的偏微分方程式KdV方程為例
右式等號右邊三項分別為平流項(A)、色散項(D)跟擴散項(S)，γ 和 ν 為色散和擴散係數。把這個方程式寫成在平流項A和色散項D在時間上用leap-frog積分、擴散項S為前方差分、且A、D、S三項在空間上皆為中央差分的數值模式可以寫成 下面的樣子
下標 k和 i 分別表示時間和空間的 index
上面的方程組的matlab code 如下
其中 u0, u1, u2皆為長度=Ndim的向量，Ndim為空間上狀態變數的長度，u0, u1, u2分別為 t=k-2, k-1 ,k 時的狀態變數。
這邊因為使用 leap-frog scheme，所以要得到 t=k 時的狀態變數 u2 時，會用到上一步 (t=k-2) 的狀態變數 u0，也就是迴圈裡的最後一行 u2 = u0 + 2 * dt * (A + D + S)。並且由於是非線性模式，所以需要現在時間的狀態變數來計算 A、D、S；A 和 D 的時間積分用 leap-frog，所以用現在時間 (t=k-1) 的狀態變數 u1 去計算；而 S 則使用前方差分所以直接用u0去計算(註)。空間上 A 跟 D 使用 3 階跟 2 階的中央差分、S 則為 1 階的中央差分。
註：如果是使用簡單的前方差分做數值積分的話A、D、S就會都使用u1，最後一行則寫成 u2 = u1 + dt * (A + D + S)。不過這個 KdV 方程式好像只用前方差分會不穩定，所以得用leap-frog。
TLM
TLM
總之不是很懂這個模式在幹麻也沒關係，只要知道狀態變數是 u0, u1, u2就好。這邊要說明的 TLM 跟 ADJ 的寫法是只要看著code一行一行進行就可以寫出來的(即使不知道每一行的意義也沒關係)。
首先在這個code裡關係到狀態變數向前積分的項為A、D、S跟迴圈裡最後一行的u2=u0+2dt(A+D+S)，主要需要加以改寫的也是這四行。
以A那一行的code為例，把它寫成式子的樣子
上面的式子便是A這一項在TLM的code裡會有的樣子，表示有一已知微小變量du1時，A這一項的相應微小變量。
其他D、S跟最後一行的u2=u0+2*dt*(A+D+S)也都照同樣的方式經過轉換，最後TLM寫出來就是
(對其他行轉換的過程有興趣的話可以參考這個檔案)
ADJ
ADJ
至於adoint code的寫法，同樣以平流項A那一行為例，除了上面導出來的TLM式子以外，需要再另外加上狀態變數du1_{ip}、du1_i、du1_{im}的式子組成一組方程組(理由等一下會講)
上式中變數前面的a表示adjoint變數，算是一種習慣寫法。
上面的矩陣寫成式子就是右邊的樣子
上面左邊的式子們就是原本模式裡A那一行寫成adjoint code時會有的樣子。
觀察得到adjoint code的過程可以發現，一開始在TLM的式子以外加上du1_ip、du1_i、du1_im 的三個式子是為了在得到 adjoint 的式子時得到 au1_{ip}=au1_{ip}+...的形式，這個式子的意思是往回算的的微小差異au1_{ip}會等於原本的微小差異 au1_{ip}加上後面一串東西。
同樣地，把code中的其他行也依照上面的程序完成每一行的adjoint code後就可以寫成整個模式的adjoint model。這邊注意寫adjoint model時順序要反過來，也就是要先寫u2=u0+2*dt*(A+D+S)的adjoint，再依序寫S、D、A的adjoint，且在最一開始要把大家都先歸零，其matlab code如下
在上面的adjoint model 的code裡，input au2為現在時間(t=k)時的一個已知微小變量，仔細觀察可以發現，adjoint code在迴圈裡做的事就是用這個已知的微小變量去計算A、D、S的微小變量 aA、aD、aA，然後再分別用aA、aD、aA去計算出狀態變數 u1、u0 的微小變量au0、au1。au0、au1即為這個adjoint model 的 output，因為在得到u2 (t=k)時同時用了u1 (t=k-1)和 u0 (t=k-2)兩個時間的狀態變數，所以在倒回去時也會同時得到 t=k-1和 t=k-2兩個時間的 adjoint 變數 (au0和au1)。並且由於這個方程式是非線性的，所以在倒回去的同時也需要 t=k-1時的狀態變數 u1這個input。
結語
結語
由上面推導的過程應該可以很清楚理解到就算是同一個方程式，如果寫成數值模式時的寫法不同(例如積分方法不同)，其相應的TLM和 adjoint code也會不同，這也是需要使用 ADJ 的 4DVAR方法其中一個比較為人詬病(?)的地方，也就是只要模式作任何更動，其相應的 ADJ 跟 TLM 也要重寫 (要不是有這點困難，4DAVR 方法可能會更為普及吧)。
在這篇文章中用了一個比較複雜的 KdV方程式的例子 (單純因為筆者最近剛好寫了一個KdV的adjoint code，直接拿來當範例了)，但做法的程序都是相同的，只要掌握狀態變數的微分、寫成矩陣、轉置的過程，應該都可以寫出相應的adjoint model。有興趣的人不妨可以試試看自己寫 Lorenz 63 的model以及相應的 TLM 跟 ADJ (當然也要根據模式積分的方式去寫出自己的TLM 跟 ADJ 喔)。
Page updated
Google Sites
Report abuse