连续过度放松

Question

这里我有一些 python 脚本，它使用 Gauss-Seidel 方法求解线性方程组：import numpy as npITERATION_LIMIT = 1000#systemA = np.array([[15., -4., -3., 8.],          [-4., 10., -4., 2.],          [-3., -4., 10., 2.],          [8., 2., 2., 12.]          ])# vector bb = 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, omegaOutput: phi (roots for linear equations)Choose an initial guess phi to the solutionrepeat 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 reachedend (repeat)有人在python上有工作算法吗？如果您可以对代码进行一些评论或帮助我如何更改此代码，那就太好了。谢谢！

Answer 1

这是基于您提供的 Wiki 参考的实现。

import numpy as np

def sor_solver(A, b, omega, initial_guess, convergence_criteria):

  """

  This is an implementation of the pseduo-code provided in the Wikipedia article.

  Inputs:

    A: nxn numpy matrix

    b: n dimensional numpy vector

    omega: relaxation factor

    initial_guess: An initial solution guess for the solver to start with

  Returns:

    phi: solution vector of dimension n

  """

  phi = initial_guess[:]

  residual = np.linalg.norm(np.matmul(A, phi) - b) #Initial residual

  while residual > convergence_criteria:

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

      sigma = 0

      for j in range(A.shape[1]):

        if j != i:

          sigma += A[i][j] * phi[j]

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

    residual = np.linalg.norm(np.matmul(A, phi) - b)

    print('Residual: {0:10.6g}'.format(residual))

  return phi

#An example case that mirrors the one in the Wikipedia article

residual_convergence = 1e-8

omega = 0.5 #Relaxation factor

A = np.ones((4, 4))

A[0][0] = 4

A[0][1] = -1

A[0][2] = -6

A[0][3] = 0

A[1][0] = -5

A[1][1] = -4

A[1][2] = 10

A[1][3] = 8

A[2][0] = 0

A[2][1] = 9

A[2][2] = 4

A[2][3] = -2

A[3][0] = 1

A[3][1] = 0

A[3][2] = -7

A[3][3] = 5

b = np.ones(4)

b[0] = 2

b[1] = 21

b[2] = -12

b[3] = -6

initial_guess = np.zeros(4)

phi = sor_solver(A, b, omega, initial_guess, residual_convergence)

print(phi)