时间戳

首页 > 解决方案 > 连续过度放松

python - 连续过度放松

问题描述

这里我有一些 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。显然，随着欧米茄值的增加，迭代次数应该会减少。这是维基百科提供的算法：

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 的工作算法吗？如果您可以对代码进行一些评论或帮助我如何更改此代码，那就太好了。谢谢！

标签： pythonpython-3.xnumpyiteration

解决方案

这是基于您提供的 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)

推荐阅读