提示
程序只是对实际数学公式的一个反映,因此在下面我都先给数学公式再给函数,当然,函数是已经封装好的可以直接用的。
Jacobi迭代法 链接到标题
先来说说他的步骤
- 设定初始值x0
- 设定迭代次数k
- 设定精度eps
- 设定迭代次数max_iter
- 设定系数矩阵A和常数项列向量b
其数学迭代公式为 $$x^{(k+1)}i=\frac{\displaystyle b_i-\sum{j=1,j\neq i}^n a_{ij}x^{(k)}j}{a{ii}}$$
def jacobi(A, b, x0, max_iter=10):
n = len(b)
x = x0.copy()
x_new = np.zeros_like(x0)
history = [x0.copy()]
for _ in range(max_iter):
for i in range(n):
s = np.dot(A[i, :], x) - A[i, i] * x[i]
x_new[i] = (b[i] - s) / A[i, i]
x = x_new.copy()
history.append(x.copy())
return x, np.array(history)
function x = jacobi(A, b, k)
% A: 系数矩阵
% b: 常数项列向量
% k: 最大迭代次数
n = length(b);
x = zeros(n, 1);
x_new = zeros(n, 1);
for iter = 1:k
for i = 1:n
sum = b(i);
for j = 1:n
if i ~= j
sum = sum - A(i,j) * x(j);
end
end
x_new(i) = sum / A(i,i);
end
x = x_new;
end
end
Gauss-Seidel迭代法 链接到标题
其数学迭代公式为 $$x^{(k+1)}i=\frac{\displaystyle b_i-\sum{j=1,j\neq i}^n a_{ij}x^{(k+1)}j-\sum{j=1}^n a_{ij}x^{(k)}j}{a{ii}}$$
def gauss_seidel(A,b,x0,max_iter=10):
n=len(b)#获取方程组的阶数
x=x0.copy()#复制初始值
x_new=np.zeros_like(x0)#初始化新的解
history=[x0.copy()]#记录迭代过程
for _ in range(max_iter):#这个带_的循环的意思是循环max_iter次,与不带_的循环不同,不需要用到循环变量
for i in range(n):
s=np.dot(A[i,:],x)-A[i,i]*x[i] #计算除了第i个元素之外的所有元素的乘积。
x[i]=(b[i]-s)/A[i,i]#计算新的解
#上面的两行代码其实是书本上端1公式2.5
history.append(x.copy())#记录迭代过程
return x,np.array(history)#返回最终解和迭代过程
SOR迭代法 链接到标题
SOR迭代法的数学迭代公式是Gauss-Seidel迭代法的一个变种,增加了一个松弛因子$\omega$,其数学迭代公式含$\omega$ 为 $$x^{(k+1)}i=(1-\omega)x^{(k)}i+\frac{\displaystyle \omega}{a{ii}}\left(b_i-\sum{j=1,j\neq i}^n a_{ij}x^{(k+1)}j-\sum{j=1}^n a_{ij}x^{(k)}_j\right)$$
def sor_solver(A, b, x_exact, omega, tol=5e-6, max_iterations=1000):
n = len(b)
x = np.zeros_like(b)
iteration = 0
iteration_process = []
for iteration in range(max_iterations):
x_new = np.copy(x)
for i in range(n):
sigma = 0
for j in range(n):
if j != i:
sigma += A[i, j] * x_new[j]
x_new[i] = (1 - omega) * x[i] + (omega / A[i, i]) * (b[i] - sigma)
iteration_process.append(np.copy(x_new))
if np.linalg.norm(x_new - x_exact, ord=np.inf) < tol:
break
x = x_new
return x, iteration + 1, iteration_process
function x = SOR(A, b, k, omega)
% A: 系数矩阵
% b: 常数项列向量
% k: 最大迭代次数
% omega: 松弛因子
n = length(b);
x = zeros(n, 1);
for iter = 1:k
for i = 1:n
sum1 = b(i);
sum2 = 0;
% 上三角部分
for j = 1:i-1
sum1 = sum1 - A(i,j) * x(j);
end
% 下三角部分
for j = i+1:n
sum2 = sum2 - A(i,j) * x(j);
end
x(i) = (1 - omega) * x(i) + (omega / A(i,i)) * (sum1 + sum2);
end
end
end