計算方法——共軛梯度法求解線性方程組

1、計算方法上機報告計算方法上機報告1 共軛梯度法求解線性方程組1.1 算法原理及程序框圖當線性方程組 Ax = b 的系數(shù)矩陣 A 是對稱正定矩陣是，可以采用共軛梯度法對該方程組進行求解，可以證明，式(1)所示的 n 元二次函數(shù)f ( x) = 1 xT Ax - bT x(1)2取得極小值點 x*是方程 Ax = b 的解。共軛梯度法是把求解線性方程組的問題轉(zhuǎn)化為求解一個與之等價的二次函數(shù)極小化的問題。從任意給定的初始點出發(fā)，沿一組關于矩陣 A 的共軛方向進行線性搜索，在無舍入誤差的假定下，最多迭代 n 次（其中 n 為矩陣 A 的階數(shù)），就可求得二次函數(shù)的極小點，也就求得線性方程組 Ax =

2、 b 的解。其迭代格式為公式(2)。5kx(k +1) = x(k ) +a d (k )(2)共軛梯度法中關鍵的兩點是迭代格式(2)中最佳步長ak 和搜索方向 d (k)的確定。其中ak 可以通過一元函數(shù)f(a) = f(x(k)+ad(k)的極小化來求得，其表達式為公式(3)；取 d (0)= r(0) = b-Ax(0)，則 d(k+1) = r(k+1) +bkd(k)，要求 d(k+1)滿足 (d(k+1) , Ad(k) = 0，可得bk 的表達式(4)。r(k )T d (k )ak = d (k )T Ad (k )(3)b =- r(k +1)T Ad (k )kd (k )

3、T Ad (k )(4)經(jīng)過一系列的證明和簡化，最終可得共軛梯度法的計算過程如下，計算程序框圖如圖 1。(1) 給定初始計算向量 x(0)即精度e >0； (2) 計算 r(0) = b -Ax(0)，取 d(0) = r(0)；(3) fork =0ton-1dor(k )T r(k )ak = d (k )T Ad (k ) ；kx(k +1) = x(k ) +a d (k ) ；r(k +1) = b + Ax(k +1) ；若| r(k +1) |£ e 或k +1 = n ，則輸出近似解 x(k+1)，停止；否則，轉(zhuǎn)；| r(k +1) |2| r(k) |2bk

4、= 2 ；2kd (k +1) = r(k +1) + b d (k ) ；end do圖 1 共軛梯度法求解線性方程組程序框圖1.2 程序使用說明共軛梯度法求解線性方程組的 matlab 程序見附錄 1，該程序可以求解系數(shù)矩陣為對稱正定矩陣的線性方程組。在使用該程序時，可將程序復制到 matlab 命令窗口中直接運行或者復制到編輯窗口中保存運行，運行時刻根據(jù)提示輸入，直至得到結(jié)果。開始運行程序時，會出現(xiàn)提示“請選擇系數(shù)矩陣、右端項以及系數(shù)矩陣階數(shù)的輸入方式：從文件中輸入數(shù)據(jù)輸入 1，從命令窗口輸入數(shù)據(jù)請輸入 2。”此時需要選擇數(shù)據(jù)輸入的方式（方式 1 和方式 2），即文件輸入（選擇 1）或者

5、手動輸入（選擇 2）。當線性方程組 Ax = b 的系數(shù)矩陣的階數(shù)較大時，將該系數(shù)矩陣 A、右端項 b 以及系數(shù)矩陣的階數(shù)n 保存為txt 格式放在D 盤的根目錄下并分別命名為data_A、data_b 和data_n， 并輸入 1 選擇文件輸入。若系數(shù)矩陣 A 的階數(shù)較小，使用手動輸入工作量小時，可在命令框中輸入 2 選擇手動輸入。選擇手動輸入時需要輸入系數(shù)矩陣 A、右端項 b 以及系數(shù)矩陣的階數(shù) n 這三個量，其輸入格式與 matlab 中矩陣、列向量和數(shù)的輸入格式要求相同。A=aijn×n 的輸入格式為a11,a12,a1n;a21,a22,a2n;an1,an2,ann，b=

6、(bi)T 的輸入格式為b1;b2;bn，n 直接輸入對應的數(shù)值即可。定義完需要求解的線性方程組之后，接下來會出現(xiàn)提示“輸入計算要求的精度 epsilon：”，按照提示輸入精度值，如要求的精度為 10-6 時輸入 1e-6 即可。以上數(shù)據(jù)的輸入全部完成，每次按提示輸入完成時按 Enter 鍵繼續(xù)下一步。在程序運行過程中，屏幕上會顯示殘差|r|2 隨著迭代步數(shù)的變化散點圖，程序運行完成之后，matlab 命令窗口會顯示線性方程的解 x，同時該解會被保存到 D 盤根目錄下名為 data_x_result 的 Excel 文件中，方便之后的數(shù)據(jù)處理（注意在求解一個新方程組之前查看 D 盤根目錄，將上

7、次得到的文件刪除，以免影響本次計算的結(jié)果）。1.3 算例計算結(jié)果(1) 數(shù)值分析課本第 29 頁的例題 2.2.2，下面采用共軛梯度法來求解線性方程組 Ax = b，其中æ 9189-27 öæ 1 öç 18450-45 ÷ç 2 ÷A = ç÷ ， b = ç÷ç 901269 ÷ç16 ÷ç -27-459135 ÷ç 8 ÷èøèø由于該方程組的系數(shù)

8、矩陣以及右端項都比較簡單，因此采用從 matlab 命令窗口手動輸入的方式來輸入數(shù)據(jù)，取計算精度為 10-6，運行過程及結(jié)果如圖 2 和圖 3（由于迭代的初始值為隨機產(chǎn)生，因此每次得到的殘量圖會有所不同，但最終都趨于 0）：圖 2 命令窗口顯示的運行結(jié)果2520殘量15105011.522.53迭代步數(shù)圖 3 殘量隨迭代步數(shù)的變化(2) 數(shù)值分析課本第 113 頁的計算實習題 3.2，用共軛梯度法求解線性方程組Ax = b，其中æ -21öæ -1öç 1-21÷ç 0 ÷0 ÷÷÷÷ç÷çA = ç÷ ， b = çç1-21 ÷çç÷çç1-2 ÷ç -1÷èøèø矩陣 A 的階數(shù)取 200 進行求解。由于該線性方程組的系數(shù)矩陣階數(shù)比較大，而且具有一定的規(guī)律，因此首先用matlab 編程將系數(shù)矩陣、右端項以及階數(shù)保存在 D 盤根目錄的三個文件中（生成系數(shù)矩陣，右端項以及階數(shù)的程序見附錄 2），然后運行共軛梯度法程序進行方程