Question

Write a script or function in Python that approximates the solution to the system ??⃗=?⃗⃗ using the Jacobi Method. The inputs should be an nxn matrix A, an n-dimensional vector?⃗⃗, a starting vector ?⃗0,an error tolerance ?, and a maximum number of iterations N. The outputs should be either an approximate solution to the system??⃗=?⃗⃗ or an error message, along with the number of iterations completed. Additionally, it would be wise to build in functionality that allows you to optionally print the current estimated solution value at each iteration.

Answer 1

import numpy as np
from numpy.linalg import *

def jacobi(A, b, x0, tol, maxiter=N):
"""
Performs Jacobi iterations to solve the line system of
equations, Ax=b, starting from an initial guess, ``x0``.

Terminates when the change in x is less than ``tol``, or
if ``maxiter`` [default=N] iterations have been exceeded.

Returns 3 variables:
1. x, the estimated solution
2. rel_diff, the relative difference between last 2
iterations for x
3. k, the number of iterations used. If k=maxiter,
then the required tolerance was not met.
"""
n = A.shape[0]
x = x0.copy()
x_prev = x0.copy()
k = 0
rel_diff = tol * 2

while (rel_diff > tol) and (k < maxiter):

for i in range(0, n):
subs = 0.0
for j in range(0, n):
if i != j:
subs += A[i,j] * x_prev[j]

x[i] = (b[i] - subs ) / A[i,i]
k += 1

rel_diff = norm(x - x_prev) / norm(x)
print(x, rel_diff)
x_prev = x.copy()

return x, rel_diff, k