第7章土壤水分運(yùn)動(dòng)的數(shù)值解

1、7.1 7.1 概述概述7.2 7.2 差分的概念及分類(lèi)差分的概念及分類(lèi)7.3 7.3 顯式差分、隱式差分和中心差分顯式差分、隱式差分和中心差分7.4 7.4 土壤水分運(yùn)移方程的線(xiàn)性化方法土壤水分運(yùn)移方程的線(xiàn)性化方法7.5 7.5 線(xiàn)性化方法與土壤水分參數(shù)的取值線(xiàn)性化方法與土壤水分參數(shù)的取值7.6 7.6 邊界條件的處理及追趕法求解三對(duì)角方程邊界條件的處理及追趕法求解三對(duì)角方程7.7 7.7 垂直一維非飽和土壤水流計(jì)算流程圖垂直一維非飽和土壤水流計(jì)算流程圖Huyakorn, P. S. and G. F. Pinder. Computational methods in subsurface

2、flow. Academic Press, New York, 1983.薛禹群，謝春紅，水文地質(zhì)學(xué)的數(shù)值法，煤炭工業(yè)出版社，薛禹群，謝春紅，水文地質(zhì)學(xué)的數(shù)值法，煤炭工業(yè)出版社，1980.土壤水分運(yùn)動(dòng)方程的求解方法：土壤水分運(yùn)動(dòng)方程的求解方法： 解析解（解析解（Analytical solution） 數(shù)值解（數(shù)值解（Numerical solution）(1) 有限差分法（有限差分法（Finite Difference method）(2) 有限單元法（有限單元法（Finite Element method）7.1 7.1 概述概述 差分的由來(lái)（差分的由來(lái)（Taylor展開(kāi)）展開(kāi)） if =

3、 (x) exists, thenkkkxkxxxxxx!)()(kkkkxkxxxxxx!) 1()()(1)(2)(0)()(xxxxxx)(0)()(xxxxxxwhen x 0, thenxxxxxxxxxxx)()(lim)()(lim00 向前差分向前差分xxxxx)()( 向后差分向后差分xxxxx)()( 中心差分（中心差分（Eqn.(1) - Eqn. (2)）xxxxxxxx)(02)()(2xxxxxx2)()(i.e.截?cái)嗾`差截?cái)嗾`差0 ( x)截?cái)嗾`差截?cái)嗾`差0 ( x)截?cái)嗾`差截?cái)嗾`差0 ( x2)0 x ABCEFxi-1xixi+1如圖所示：如圖所示： BC表示

4、向前差分；表示向前差分； AB表示向后差分；表示向后差分； AC表示中心差分；表示中心差分；In addition, Eqn.(1) + Eqn.(2) )(0)(2)()(22222xxxxxxxxx)(0)()(2)(2222xxxxxxxxi.e.:其截?cái)嗾`差亦為與其截?cái)嗾`差亦為與 x2同階的無(wú)窮小量。同階的無(wú)窮小量。i=1, z(1)=0i=2i=3.i=n, z(n)=Lz考虛垂直一維問(wèn)題（如右圖所示），考虛垂直一維問(wèn)題（如右圖所示），泛定方程可寫(xiě)為（以泛定方程可寫(xiě)為（以 方程為例）：方程為例）：剖面深度為剖面深度為L(zhǎng)，共剖分出，共剖分出n個(gè)節(jié)點(diǎn)（個(gè)節(jié)點(diǎn)（i = 1, 2, 3, ,

5、 n），），z (1) = 0, z (n) = L，已，已知知t = tk時(shí)刻各節(jié)點(diǎn)含水率時(shí)刻各節(jié)點(diǎn)含水率 的分布，求的分布，求t = t k+1時(shí)刻各節(jié)點(diǎn)含水率時(shí)刻各節(jié)點(diǎn)含水率 的分布，的分布， kkttt1zKzDzt 7.3.17.3.1顯式差分格式顯式差分格式kiizKzDztttkikii121212121iikikikizzzDzDzDz iikikikikizzDzD112121 112121iikikikikizzDzD7.3.17.3.1顯式差分格式顯式差分格式 2112111211121 zzzzDzzDzDziiikikikiiikikikiki21212121iiki

6、kikizzKKzK21212121112111211iiiiiikikiiiikikiikikizzKKzzDzzDt)(21112121iiiizzzz) 1, 3, 2( 2212111211121111niKKzzDzzDzztiiiikikiiiikikiiiikiki7.3.17.3.1顯式差分格式顯式差分格式 若取等步長(zhǎng)若取等步長(zhǎng) z，即，即 z =zi+1 - zi = zi - zi-1，則有：，則有：) 1, 3, 2(212112112121niKKzDDztiikikiikikiikiki7.3.17.3.1顯式差分格式顯式差分格式 顯式差分格式是有條件收斂的顯式差分格

7、式是有條件收斂的，一般應(yīng)滿(mǎn)足，一般應(yīng)滿(mǎn)足max2Dzrt其中：其中：r多取多取1/2。由于土壤接近飽和時(shí)，。由于土壤接近飽和時(shí)，Dmax很很大，故大，故 t一般要求取值很小，耗費(fèi)機(jī)時(shí)。一般要求取值很小，耗費(fèi)機(jī)時(shí)。7.3.17.3.1顯式差分格式顯式差分格式 1kiiKzDzt) 1, 3, 2( 22/12/1111121111121111niKKzzDzzDzztiiiikikiiiikikiiiikiki7.3.27.3.2隱式差分格式隱式差分格式 )(1111iiiizzzztr)(1112iiiizzzztrLetikiikiikiiHGFE11111then ( i = 2, 3,

8、, n - 1)三對(duì)角方程三對(duì)角方程7.3.27.3.2隱式差分格式隱式差分格式 where;211212iiiiDrGDrE11211212iiiiiDrDrGEF2121112iiiikiiKKzztH隱式差分格式是無(wú)條件收斂的。隱式差分格式是無(wú)條件收斂的。7.3.27.3.2隱式差分格式隱式差分格式 kikikiiKzDzKzDzKzDzt12/1217.3.37.3.3中心差分格式中心差分格式 ) 1, 3, 2()(2)(22212111111211111121111niKKzzDzzDzztiiiikikikikiiiikikikikiiiikiki7.3.37.3.3中心差分格式

9、中心差分格式 )(1111iiiizzzztr)(1112iiiizzzztrLetikiikiikiiHGFE11111then ( i = 2, 3, , n - 1)where;211212iiiiDrGDrE7.3.37.3.3中心差分格式中心差分格式 11211212iiiiiDrDrGEF中心差分格式也是無(wú)條件收斂的。中心差分格式也是無(wú)條件收斂的。212111112) 1(iiiikiikiiikiiiKKzztGGEEH7.3.37.3.3中心差分格式中心差分格式 7.4 7.4 土壤水分運(yùn)移方程的線(xiàn)性化方法土壤水分運(yùn)移方程的線(xiàn)性化方法 以基質(zhì)勢(shì)水頭以基質(zhì)勢(shì)水頭h為因變量的一維垂

10、直向土壤為因變量的一維垂直向土壤水分運(yùn)動(dòng)基本方程為：水分運(yùn)動(dòng)基本方程為： K hhhC hK htzzz7.4 7.4 土壤水分運(yùn)移方程的線(xiàn)性化方法土壤水分運(yùn)移方程的線(xiàn)性化方法7.4.17.4.1顯式差分格式顯式差分格式相應(yīng)h方程的差分方程為：11/211/21211kkkkkkkkkiiiiiiiiiChhKhhKhhtz1/21/2kkiiKKz7.4 7.4 土壤水分運(yùn)移方程的線(xiàn)性化方法土壤水分運(yùn)移方程的線(xiàn)性化方法7.4.27.4.2隱式差分格式隱式差分格式相應(yīng)h方程的差分方程為1111111111/211/2111/21/22kkkkkkkkkkiiiiiiiikiiihhKhhKhh

11、KKCtzz7.4 7.4 土壤水分運(yùn)移方程的線(xiàn)性化方法土壤水分運(yùn)移方程的線(xiàn)性化方法7.4.37.4.3中心差分格式中心差分格式相應(yīng)h方程的差分方程為1111111221221122111111122212222kkkkkkkkiiiiiiiikkkkkkiikiiiiihhhhhhCKtzKKhhhhKzz 土壤水分運(yùn)移方程中的各參數(shù)均依賴(lài)于變量土壤水分運(yùn)移方程中的各參數(shù)均依賴(lài)于變量 或或 m，從而使得方程呈現(xiàn)出較強(qiáng)的非線(xiàn)性性，求解前必須將其線(xiàn)從而使得方程呈現(xiàn)出較強(qiáng)的非線(xiàn)性性，求解前必須將其線(xiàn)性化，得到性化，得到n元一次線(xiàn)性代數(shù)方程組，以便于求解。元一次線(xiàn)性代數(shù)方程組，以便于求解。（1）顯式

12、線(xiàn)性化）顯式線(xiàn)性化 計(jì)算過(guò)程中，方程中的各參數(shù)如：計(jì)算過(guò)程中，方程中的各參數(shù)如：C ( m)，K( m) or D( )，K( )等均以時(shí)段初的值（即前一時(shí)段的值）代入。等均以時(shí)段初的值（即前一時(shí)段的值）代入。 適用條件：此法只能用于適用條件：此法只能用于 m或或 變化緩慢的情況；當(dāng)變化緩慢的情況；當(dāng)變化劇烈時(shí)，此法會(huì)導(dǎo)致較大的偏差。變化劇烈時(shí)，此法會(huì)導(dǎo)致較大的偏差。7.5.17.5.1線(xiàn)性化方法線(xiàn)性化方法例如：已知例如：已知k - 1與與k時(shí)刻各節(jié)點(diǎn)的時(shí)刻各節(jié)點(diǎn)的 值，需求值，需求k + 1時(shí)刻各時(shí)刻各節(jié)點(diǎn)的節(jié)點(diǎn)的 值，應(yīng)用顯式外推線(xiàn)性化參數(shù)的方法如下（以值，應(yīng)用顯式外推線(xiàn)性化參數(shù)的方法如下

13、（以導(dǎo)水率導(dǎo)水率K的求法為例，其他各參數(shù)的求法與此類(lèi)似）：的求法為例，其他各參數(shù)的求法與此類(lèi)似）：1111)(kkkkkikikikitttt(*) 求求 ： 方法（方法（a）：）：21kiK)(21121kikiki2121kikiKK方法方法(b)：算術(shù)平均算術(shù)平均:幾何平均幾何平均:調(diào)和平均調(diào)和平均:121)(21kikikiKKK)()(121kikikiKKK)()()()(21121kikikikikiKKKKK 求求 ： 方法（方法（a）：）：kiK21)(21121kikikikikiKK2121方法方法(b)：算術(shù)平均算術(shù)平均幾何平均幾何平均調(diào)和平均調(diào)和平均kikikiKKK

14、121)(21)()(121kikikiKKK)()()()(21121kikikikikiKKKKK 計(jì)算時(shí)，先用顯式線(xiàn)性化的方法，采用隱式差分格式，計(jì)算時(shí)，先用顯式線(xiàn)性化的方法，采用隱式差分格式，求出求出k + 1/2時(shí)刻各節(jié)點(diǎn)的時(shí)刻各節(jié)點(diǎn)的 m或或 （即（即 or ）；根據(jù)）；根據(jù)這些值可求得相應(yīng)的這些值可求得相應(yīng)的 or 等參數(shù)，以此作為計(jì)算時(shí)始末（等參數(shù)，以此作為計(jì)算時(shí)始末（tk to tk+1）的平均值，代）的平均值，代入方程求解入方程求解 ，從而完成方程的線(xiàn)性化工作。，從而完成方程的線(xiàn)性化工作。21kmi21ki2121kmikiKK21kiK1ki7.5.17.5.1線(xiàn)性化方法

15、線(xiàn)性化方法7.5.1.27.5.1.2預(yù)報(bào)校正法預(yù)報(bào)校正法（1）顯式外推法（先求顯式差分方程，得）顯式外推法（先求顯式差分方程，得 作為預(yù)報(bào)作為預(yù)報(bào)值）值） 計(jì)算參數(shù)時(shí)，據(jù)前一時(shí)段始末的函數(shù)值計(jì)算參數(shù)時(shí)，據(jù)前一時(shí)段始末的函數(shù)值 m或或 ，用，用線(xiàn)性外推近似求出本時(shí)段末的函數(shù)值，然后線(xiàn)性外推近似求出本時(shí)段末的函數(shù)值，然后 方法方法(a)：求出計(jì)算時(shí)段始末的平均值，再通過(guò)各參數(shù)與：求出計(jì)算時(shí)段始末的平均值，再通過(guò)各參數(shù)與 m或或 的關(guān)系，求出各參數(shù)；或的關(guān)系，求出各參數(shù)；或 方法方法(b)：先求出計(jì)算時(shí)段始末的參數(shù)值，然后對(duì)參數(shù)?。合惹蟪鲇?jì)算時(shí)段始末的參數(shù)值，然后對(duì)參數(shù)取平均。平均。 平均值的計(jì)算

16、可采用算術(shù)平均、幾何平均、調(diào)和平平均值的計(jì)算可采用算術(shù)平均、幾何平均、調(diào)和平均等方法。均等方法。7.5.17.5.1線(xiàn)性化方法線(xiàn)性化方法7.5.1.27.5.1.2預(yù)報(bào)校正法預(yù)報(bào)校正法1ki（2）線(xiàn)性外推法（已知線(xiàn)性外推法（已知 時(shí)刻，外推時(shí)刻，外推 時(shí)刻）時(shí)刻） 假定含水率在相鄰時(shí)段內(nèi)是線(xiàn)性變化的，根據(jù)前假定含水率在相鄰時(shí)段內(nèi)是線(xiàn)性變化的，根據(jù)前一個(gè)時(shí)段末的含水率，用直線(xiàn)外推近似求出計(jì)算時(shí)段一個(gè)時(shí)段末的含水率，用直線(xiàn)外推近似求出計(jì)算時(shí)段末的含水率：末的含水率：7.5.17.5.1線(xiàn)性化方法線(xiàn)性化方法7.5.1.27.5.1.2預(yù)報(bào)校正法預(yù)報(bào)校正法,1k k 1k 112kkkiii11231

17、22kkkiii 對(duì)于第一時(shí)段，無(wú)法應(yīng)用線(xiàn)性外推法，可用其他對(duì)于第一時(shí)段，無(wú)法應(yīng)用線(xiàn)性外推法，可用其他線(xiàn)性化方法。線(xiàn)性化方法。 7.5.17.5.1線(xiàn)性化方法線(xiàn)性化方法7.5.1.27.5.1.2預(yù)報(bào)校正法預(yù)報(bào)校正法（3 3）顯式預(yù)測(cè)法（兩次求解隱式：第一次預(yù)報(bào)值；）顯式預(yù)測(cè)法（兩次求解隱式：第一次預(yù)報(bào)值；第二次得校正值）第二次得校正值） 首先利用顯式線(xiàn)性化方法，由時(shí)段初的含水率首先利用顯式線(xiàn)性化方法，由時(shí)段初的含水率及相應(yīng)的參數(shù)值，按隱式差分格式求解方程，求及相應(yīng)的參數(shù)值，按隱式差分格式求解方程，求得本時(shí)段末（或時(shí)段中間）的含水率及相應(yīng)的參得本時(shí)段末（或時(shí)段中間）的含水率及相應(yīng)的參數(shù)值（預(yù)報(bào)

18、值），再求解隱式差分方程，得到時(shí)數(shù)值（預(yù)報(bào)值），再求解隱式差分方程，得到時(shí)段末的含水率的值（校正值）。段末的含水率的值（校正值）。計(jì)算時(shí)，先假定本時(shí)段末（計(jì)算時(shí)，先假定本時(shí)段末（tk+1）的）的 m或或 （一般做（一般做法是先取時(shí)段初的函數(shù)值法是先取時(shí)段初的函數(shù)值 or ） ， or 等參數(shù)的迭代初值等參數(shù)的迭代初值 求解線(xiàn)性代數(shù)方程組得求解線(xiàn)性代數(shù)方程組得 or 直至前后兩次迭代值的誤差滿(mǎn)足精度要求為止。直至前后兩次迭代值的誤差滿(mǎn)足精度要求為止。kmiki)1(1kiK)1(1kiD)1(1kiC)1(1kmi)1 (1ki)2(1)2(1)2(1or ,kikikiCDK)2(1)2(1o

19、r kikmi)(1)(1or pkipkmi7.5.1.3迭代法迭代法（a）若取絕對(duì)誤差，判斷條件為：）若取絕對(duì)誤差，判斷條件為：), 2, 1(max1) 1( 1)( 1nipkipkii), 2, 1(max1)1(1)(1nipkmipkmiior7.5.17.5.1線(xiàn)性化方法線(xiàn)性化方法7.5.1.3迭代法迭代法（b）若取相對(duì)誤差，判斷條件為：）若取相對(duì)誤差，判斷條件為：), 2, 1(max2)(1)1(1)(1nipkipkipkii), 2, 1(max2)(1)1(1)(1nipkmipkmipkmiior7.5.17.5.1線(xiàn)性化方法線(xiàn)性化方法7.5.1.3迭代法 1與與

20、2 的取值可根據(jù)實(shí)際問(wèn)題的精度要求而定。的取值可根據(jù)實(shí)際問(wèn)題的精度要求而定。 此法計(jì)算稍復(fù)雜，但精度較高，數(shù)值計(jì)算時(shí)此法計(jì)算稍復(fù)雜，但精度較高，數(shù)值計(jì)算時(shí)經(jīng)常采用。為避免迭代次數(shù)太多，一般可設(shè)置一經(jīng)常采用。為避免迭代次數(shù)太多，一般可設(shè)置一個(gè)判斷語(yǔ)句：當(dāng)?shù)螖?shù)個(gè)判斷語(yǔ)句：當(dāng)?shù)螖?shù)m5時(shí)，將時(shí)間步長(zhǎng)時(shí)，將時(shí)間步長(zhǎng) 減半，再重新計(jì)算。減半，再重新計(jì)算。kkkttt117.5.17.5.1線(xiàn)性化方法線(xiàn)性化方法7.5.1.3迭代法7.5 線(xiàn)性化方法與土壤水分參數(shù)的取值線(xiàn)性化方法與土壤水分參數(shù)的取值7.5.2.1直接法直接法先求出兩結(jié)點(diǎn)的參數(shù)，然后先求出兩結(jié)點(diǎn)的參數(shù)，然后(1)取兩結(jié)點(diǎn)參數(shù)的算術(shù)平均值

21、：取兩結(jié)點(diǎn)參數(shù)的算術(shù)平均值：7.5.2參數(shù)取值參數(shù)取值11212iiiDDD(2)取兩結(jié)點(diǎn)參數(shù)的調(diào)和平均值：取兩結(jié)點(diǎn)參數(shù)的調(diào)和平均值：11122/iiiiiDDDDD(3) 取兩結(jié)點(diǎn)參數(shù)的幾何平均值：取兩結(jié)點(diǎn)參數(shù)的幾何平均值：112iiiDDD7.5 線(xiàn)性化方法與土壤水分參數(shù)的取值線(xiàn)性化方法與土壤水分參數(shù)的取值7.5.2 2間接法間接法先求出兩結(jié)點(diǎn)處含水率的算術(shù)、調(diào)和、幾何平均值，先求出兩結(jié)點(diǎn)處含水率的算術(shù)、調(diào)和、幾何平均值，然后再求出相應(yīng)的參數(shù)值。然后再求出相應(yīng)的參數(shù)值。(1)算術(shù)平均值：算術(shù)平均值：11212iii1122iiDD(2)調(diào)和平均值：調(diào)和平均值：11122/i iiii112

22、2iiDD(3)幾何平均值：幾何平均值：112i ii1122iiDD7.5 線(xiàn)性化方法與土壤水分參數(shù)的取值線(xiàn)性化方法與土壤水分參數(shù)的取值7.5.2 3其他方法其他方法（1）求出的參數(shù)和兩結(jié)點(diǎn)處的參數(shù)和，三點(diǎn)參）求出的參數(shù)和兩結(jié)點(diǎn)處的參數(shù)和，三點(diǎn)參數(shù)取平均，故可稱(chēng)為數(shù)取平均，故可稱(chēng)為“三點(diǎn)式三點(diǎn)式”：11122124iiiiDDDD（2） 11121iiiiiDDd （3）由鄰近結(jié)點(diǎn)的參數(shù)代替，如）由鄰近結(jié)點(diǎn)的參數(shù)代替，如11122,iiiiDD DD前已敘及，采用顯式差分格式可直接求解，采前已敘及，采用顯式差分格式可直接求解，采用隱式或中心差分格式時(shí)總可以得到如下形式的三對(duì)用隱式或中心差分格

23、式時(shí)總可以得到如下形式的三對(duì)角方程組（以角方程組（以 方程為例）：方程為例）：共共n - 2個(gè)方程，再加上上、下邊界條件，即可組成個(gè)方程，再加上上、下邊界條件，即可組成n n 階的階的n元一次方程組，線(xiàn)性化后求解該方程組可獲得元一次方程組，線(xiàn)性化后求解該方程組可獲得原定解問(wèn)題的解。原定解問(wèn)題的解。) 1, 3 , 2( 11111niHGFEikiikiikii7.6.17.6.1邊界條件的處理邊界條件的處理 上、下邊界不論為第一、第二或第三類(lèi)邊界，上、下邊界不論為第一、第二或第三類(lèi)邊界，總可以化為：總可以化為：1121111HGFkknknnknnHFE111(2)(1)7.6.17.6.1

24、邊界條件的處理邊界條件的處理 驗(yàn)驗(yàn) 證：證： 當(dāng)上、下邊界為第一類(lèi)邊界時(shí)當(dāng)上、下邊界為第一類(lèi)邊界時(shí))()()(),()(), 0(121111121kknkkztftftftLtft)(, 0, 1)2.(Eqn)(, 0, 1) 1.(Eqn12n11111knnktfHEFtfHGF7.6.17.6.1邊界條件的處理邊界條件的處理1121111HGFkknknnknnHFE111 當(dāng)上、下邊界為第二、三類(lèi)邊界時(shí)當(dāng)上、下邊界為第二、三類(lèi)邊界時(shí))4()()3(), 0(210tqKzDttqKzDzLzz)(),(122/111112111112312111223knnnknknnkkkktq

25、KzzDtqKzzD1121111HGFkknknnknnHFE111211212123111111112231)(,),(,nknnnnnnnkkKtqHEFzzDEKtqHFGzzDFnknnknnnknnknnknnkkkkkkkkHFEHGFEHGFEHGFEHGF111111111121314313312321321221121121111經(jīng)線(xiàn)性化后，原定解問(wèn)題可化為經(jīng)線(xiàn)性化后，原定解問(wèn)題可化為n n階線(xiàn)性代數(shù)方程組，即：階線(xiàn)性代數(shù)方程組，即：以矩陣的形式表示即為：以矩陣的形式表示即為：nnknknkkknnnnnHHHHHFEGFEGFEGFEGF13211111312111113

26、33222117.6.27.6.2 追趕法求解三對(duì)角方程追趕法求解三對(duì)角方程1原理原理1121111HGFkk(1)(1a) ,( 111111111111HBFAAGABkk7.6.27.6.2 追趕法求解三對(duì)角方程追趕法求解三對(duì)角方程1原理原理2132122112HGFEkkk1122132121212ABEHGAEGFkk13222212kkAGAB(2)(2a)1122212122,ABEHBAEGFA.1111111iiiikiikiiiiiABEHGAEGF1111,iiiiiiiiiiABEHBAEGFA.111kiiiiikiAGAB(ia)(i)ikiikiikiiHGFE1

27、11117.6.27.6.2 追趕法求解三對(duì)角方程追趕法求解三對(duì)角方程111111121nknnknnknnHGFE221111112121nnnnknnknnnnnABEHGAEGF1111111knnnnnknAGAB2211121211,nnnnnnnnnnABEHBAEGFA(n-1)a)(n-1)7.6.27.6.2 追趕法求解三對(duì)角方程追趕法求解三對(duì)角方程nknnknnHFE11111111nnnnknnnnnABEHAEGF11111nnnnnnnnknAEGFABEH(na)(n)7.6.27.6.2 追趕法求解三對(duì)角方程追趕法求解三對(duì)角方程 解出解出 后，利用式（后，利用式（

28、(n-1)a）,，(ia),，(2a), (1a)往前遞推，依次即可解出往前遞推，依次即可解出 ，于，于是得到是得到k+1時(shí)刻各節(jié)點(diǎn)的時(shí)刻各節(jié)點(diǎn)的 值（值（i = 1, 2, ，n）。）。 若若 均為已知（第一類(lèi)邊界條件），則可直均為已知（第一類(lèi)邊界條件），則可直接從式（接從式（2）開(kāi)始求解，依次遞推至式（）開(kāi)始求解，依次遞推至式（n - 1）即可。）即可。1kn11121211 , , , , kkknkn111 , knk7.6.27.6.2 追趕法求解三對(duì)角方程追趕法求解三對(duì)角方程(1) 輸入輸入D、K含水量含水量SITA關(guān)系（建議采用關(guān)系（建議采用函數(shù)函數(shù)子程序調(diào)用）；子程序調(diào)用）；節(jié)

29、點(diǎn)剖分總數(shù)節(jié)點(diǎn)剖分總數(shù)n及對(duì)應(yīng)的空間坐標(biāo)及對(duì)應(yīng)的空間坐標(biāo)z(1:n)（若為等步長(zhǎng)，則（若為等步長(zhǎng)，則輸入空間步長(zhǎng)輸入空間步長(zhǎng)DZ）；）；初始含水量剖面分布初始含水量剖面分布SITA0（1:n）；）；上、下邊界條件上、下邊界條件UBC、LBC（=1，2，or3）及相應(yīng)的邊界）及相應(yīng)的邊界值值QSUBC、QSLBCT（建議采用子程序調(diào)用）；（建議采用子程序調(diào)用）；初始時(shí)間步長(zhǎng)初始時(shí)間步長(zhǎng)DT0、最大時(shí)間步長(zhǎng)、最大時(shí)間步長(zhǎng)DTMAX，計(jì)算結(jié)束時(shí)，計(jì)算結(jié)束時(shí)間間TEND，迭代控制標(biāo)準(zhǔn)，迭代控制標(biāo)準(zhǔn)DET。(2) 賦初值賦初值T(0)=0 初始時(shí)間初始時(shí)間T(1)=DT0 第一次計(jì)算時(shí)間第一次計(jì)算時(shí)間*（以下賦含水量計(jì)算中間值）（以下賦含水量計(jì)算中間值） do i=1, nSITAM(i)=SITA0(i) enddo時(shí)段計(jì)數(shù)時(shí)段計(jì)數(shù)KT=1(3) 判斷計(jì)算時(shí)間是否已超過(guò)需模擬時(shí)間判斷計(jì)算時(shí)間是否已超過(guò)需模擬時(shí)間if (T(KT)=DTMAX) then DT=DTMAX T(KT)=T(KT-1)+DT endif(5) 開(kāi)始計(jì)算系數(shù)矩陣，求解線(xiàn)性代數(shù)方程組開(kāi)始計(jì)算系數(shù)矩陣，求解線(xiàn)性代數(shù)方程組 賦迭代計(jì)次數(shù)賦迭代計(jì)次數(shù) MP1=0(6) if (MP1=5)then (else goto (12)(7) 計(jì)算本時(shí)