连续过度放松

Question

这里我有一些 python 脚本，它使用 Gauss-Seidel 方法求解线性方程组：

import numpy as np

ITERATION_LIMIT = 1000

#system

A = np.array([[15., -4., -3., 8.],

          [-4., 10., -4., 2.],

          [-3., -4., 10., 2.],

          [8., 2., 2., 12.]

          ])

# vector b

b = np.array([2., -12., -4., 6.])

print("System of equations:")

for i in range(A.shape[0]):

    row = ["{0:3g}*x{1}".format(A[i, j], j + 1) for j in range(A.shape[1])]

    print("[{0}] = [{1:3g}]".format(" + ".join(row), b[i]))

x = np.zeros_like(b)

for it_count in range(1, ITERATION_LIMIT):

    x_new = np.zeros_like(x)

    print("Iteration {0}: {1}".format(it_count, x))

    for i in range(A.shape[0]):

        s1 = np.dot(A[i, :i], x_new[:i])

        s2 = np.dot(A[i, i + 1:], x[i + 1:])

        x_new[i] = (b[i] - s1 - s2) / A[i, i]

    if np.allclose(x, x_new, rtol=1e-8):

        break

    x = x_new

它输出的是：

Iteration 379: [-21.36409652 -22.09743    -19.9999946   21.75896845]

Iteration 380: [-21.36409676 -22.09743023 -19.99999481  21.75896868]

Iteration 381: [-21.36409698 -22.09743045 -19.99999501  21.7589689 ]

我的任务是利用此方法制作连续过度松弛 (SOR) 方法，该方法使用 omega 值来减少迭代次数。如果omega = 1，则变为 Gauss-Seidel 方法、if < 1- 简单迭代法> 1和< 2- SOR。显然，随着 omega 值的增加，迭代次数应该减少。这是维基百科提供的算法：

Inputs: A, b, omega

Output: phi (roots for linear equations)

Choose an initial guess phi to the solution

repeat until convergence

  for i from 1 until n do

    sigma <- 0

    for j from 1 until n do

      if j ≠ i then

        sigma <- sigma + A[i][j]*phi[j]

      end if

    end (j-loop)

    phi[i] = phi[i] + omega*((b[i] - sigma)/A[i][i]) - phi[i]

  end (i-loop)

  check if convergence is reached

end (repeat)

有人在python上有工作算法吗？如果您可以对代码进行一些评论或帮助我如何更改此代码，那就太好了。谢谢！

Answer 1

这是基于您提供的 Wiki 参考的实现。import numpy as npdef sor_solver(A, b, omega, initial_guess, convergence_criteria):&nbsp; """&nbsp; This is an implementation of the pseduo-code provided in the Wikipedia article.&nbsp; Inputs:&nbsp; &nbsp; A: nxn numpy matrix&nbsp; &nbsp; b: n dimensional numpy vector&nbsp; &nbsp; omega: relaxation factor&nbsp; &nbsp; initial_guess: An initial solution guess for the solver to start with&nbsp; Returns:&nbsp; &nbsp; phi: solution vector of dimension n&nbsp; """&nbsp; phi = initial_guess[:]&nbsp; residual = np.linalg.norm(np.matmul(A, phi) - b) #Initial residual&nbsp; while residual > convergence_criteria:&nbsp; &nbsp; for i in range(A.shape[0]):&nbsp; &nbsp; &nbsp; sigma = 0&nbsp; &nbsp; &nbsp; for j in range(A.shape[1]):&nbsp; &nbsp; &nbsp; &nbsp; if j != i:&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; sigma += A[i][j] * phi[j]&nbsp; &nbsp; &nbsp; phi[i] = (1 - omega) * phi[i] + (omega / A[i][i]) * (b[i] - sigma)&nbsp; &nbsp; residual = np.linalg.norm(np.matmul(A, phi) - b)&nbsp; &nbsp; print('Residual: {0:10.6g}'.format(residual))&nbsp; return phi#An example case that mirrors the one in the Wikipedia articleresidual_convergence = 1e-8omega = 0.5 #Relaxation factorA = np.ones((4, 4))A[0][0] = 4A[0][1] = -1A[0][2] = -6A[0][3] = 0A[1][0] = -5A[1][1] = -4A[1][2] = 10A[1][3] = 8A[2][0] = 0A[2][1] = 9A[2][2] = 4A[2][3] = -2A[3][0] = 1A[3][1] = 0A[3][2] = -7A[3][3] = 5b = np.ones(4)b[0] = 2b[1] = 21b[2] = -12b[3] = -6initial_guess = np.zeros(4)phi = sor_solver(A, b, omega, initial_guess, residual_convergence)print(phi)