ax
=b經過一定的變換成x=
bx+f
，然後從初始向量出發，計算xk
+1=b
xk+f
，經過一定的次數後得到xk
+1會收斂於真正的值。問題來了？如何得到x=bx+f這種形式？如何證明收斂？接下來的幾個演算法都是圍繞這個問題。
# -*- coding: utf-8 -*-
import numpy as np # a package to calculate
from scipy.sparse import identity # to get a identity matrix
import time # calculate time of diffient method
np.random.seed(2) # set seed to make x0 unchanged
defcreate_tridiagonal_matrices
(a, b, c, n):
# 建立一種特殊的三對角矩陣，主對角線元素b，比主對角線低一行的對角線上a，比主對角線高一行的對角線上c
a = np.zeros((n, n))
if (b <= c or b <= a or b < a + c or n < 2): # 判斷是否符合三對角矩陣的基本定義,以及n的個數，1行的三對角毫無意義
print
"this is not a create_tridiagonal_matrices,please check a,b,c "
a[0][0] = b;
a[0][1] = c
a[n - 1][n - 1] = b;
a[n - 1][n - 2] = a
for j in range(1, n - 1):
a[j][j - 1] = a
a[j][j] = b
a[j][j + 1] = c
return a
defjacobi
(a, b):
n = a.shape[0]#初始化dlu為0
d = np.zeros(a.shape, dtype="float")
l = np.zeros(a.shape, dtype="float")
u = np.zeros(a.shape, dtype="float")
for i in range(n):#根據規則a=d-l-u計算
for j in range(i + 1, n):
u[i][j] = -a[i][j]
l[j][i] = -a[j][i]
for i in range(n):
d[i][i] = a[i][i]
x_k1 = np.zeros((n, 1))#初始化x_k1
#x_k2 = np.random.randn(n, 1) # 初始化x_k2，無所謂
d_inv = np.linalg.inv(d)
b = np.dot(d_inv, l + u)#計算相應的b和f
f = np.dot(d_inv, b)
#x_k2 = np.dot(b, x_k1) + f
while (np.sqrt(np.sum(np.power(np.dot(a,x_k1)-b, 2)))/np.sqrt(np.sum(np.power(b, 2))) > 1e-6): #迭代法不斷的趨近
#x_k1 = x_k2
x_k1 = np.dot(b, x_k1) + f
return x_k1
defgauss_seidel
(a, b):
#思路類似，在迭代的過程中立即用到了上乙個解
# d=np.diag(np.diag(a),0)
n = a.shape[0]
d = np.zeros(a.shape, dtype="float")
l = np.zeros(a.shape, dtype="float")
u = np.zeros(a.shape, dtype="float")
for i in range(n):#根據規則a=d-l-u計算
for j in range(i + 1, n):
u[i][j] = -a[i][j]
l[j][i] = -a[j][i]
for i in range(n):
d[i][i] = a[i][i]
# 初始化x_k1，k2
x_k1 = np.zeros((n, 1))
#x_k2 = np.random.randn(n, 1)
d_minus_l_inv = np.linalg.inv(d - l)
b_g = np.dot(d_minus_l_inv, u) #算出相應b
f_g = np.dot(d_minus_l_inv, b)#算出相應f
#x_k2 = np.dot(b_g, x_k1) + f_g
while (np.sqrt(np.sum(np.power(np.dot(a, x_k1) - b, 2))) / np.sqrt(np.sum(np.power(b, 2))) > 1e-6):
#x_k1 = x_k2
x_k1 = np.dot(b_g, x_k1) + f_g
return x_k1
defs_o_r
(a, b, w):
#思路類似，增加了w來加速迭代過程
# d=np.diag(np.diag(a),0)
n = a.shape[0]
d = np.zeros(a.shape, dtype="float")
l = np.zeros(a.shape, dtype="float")
u = np.zeros(a.shape, dtype="float")
for i in range(n):
for j in range(i + 1, n):
u[i][j] = -a[i][j]
l[j][i] = -a[j][i]
for i in range(n):
d[i][i] = a[i][i]
x_k1 = np.zeros((n, 1))
d_minus_wl_inv = np.linalg.inv(d - w * l)
b_w = np.dot(d_minus_wl_inv, (1 - w) * d + w * u)#算出相應b
f_w = w * np.dot(d_minus_wl_inv, b)#算出相應f
while (np.sqrt(np.sum(np.power(np.dot(a, x_k1) - b, 2))) / np.sqrt(np.sum(np.power(b, 2))) > 1e-6):
x_k1 = np.dot(b_w, x_k1) + f_w
return x_k1
deftimecmp
(a,b,fun1,fun2,fun3):
#通過執行100次的平均時間較各個函式的效率
time1=np.zeros((10,1))
time2 = np.zeros((10, 1))
time3 = np.zeros((10, 1))
for i in range(10):
t1=time.clock()
jacobi(a, b)
time1[i][0]=time.clock()-t1
t2 = time.clock()
gauss_seidel(a, b)
time2[i][0] = time.clock() - t2
t3 = time.clock()
s_o_r(a, b, 1.5)
time3[i][0] = time.clock() - t3
return np.mean(time1),np.mean(time2),np.mean(time3)
n=100
a = create_tridiagonal_matrices(-1, 2, -1, n)
print
"\n 生成的係數行列式是", a
b = np.ones((n, 1))
print
"\n 生成的b是", b
print
"\n 呼叫函式算出來解為", np.linalg.solve(a, b)
t1 = time.clock()
x = jacobi(a, b)
print
"計算花費時間",time.clock()-t1
print
'jacobi計算出來的解:',x
t1 = time.clock()
x = gauss_seidel(a, b)
print
"\n計算花費時間",time.clock()-t1
print
'gauss_seidel計算出來的解:',x
print x
t1 = time.clock()
x = s_o_r(a, b, 0.5)
print
"\n計算花費時間",time.clock()-t1
print
's_o_r計算出來的解:',x
print
"over"
print
"jacobi,gauss_seidel,s_o_r 在10次執行的平均時間",timecmp(a,b,jacobi,gauss_seidel,s_o_r)
C 解方程組之Jacobi迭代法
迭代過程 首先將 方程組中的 係數矩陣 a分解成三部分，即 a l d u，如圖1所示，其中 d為對角陣，l為下三角矩陣，u為上三角矩陣。之後確定迭代格式，x k 1 b x k f 這裡 表示的是上標，括號內數字即迭代次數 如圖2所示，其中 b稱為迭代矩陣，雅克比迭代法中一般記為 j。k 0,1，...
計算方法 迭代法 牛頓法求解方程組
迭代法 給定乙個 方程 f x 0 可以用多種方式來構造它的等價方程 x p x 取定根的乙個近似值 x0,構造序列 xk 1 p xk k 0,1 2 迭代法演算法 1.給定初始值 x0 和精度要求 e 以及最大迴圈次數 k 2.計算 xk 1 p xk 3.xk 1 xk e 如果成立則說明達到...
雅克比迭代法與高斯塞德爾迭代法求解方程組（C語言）
分別用雅可比 迭代法與高斯塞德爾迭代法解下列方程組 雅可比迭代法 include include define eps 1e 6 define max 100 雅可比迭代法 void jacobi float a,int n,float x y i a i n 1 n s a i n 1 i eps...