>> A = [9,1,1;2,10,3;3,4,11] A = 9 1 1 2 10 3 3 4 11 >> b = rand(3,1) b = 0.9501 0.2311 0.6068 Here's matlab's approximation of the true solution of Ax=b: >> inv(A)*b ans = 0.1030 -0.0063 0.0294 >> help jacobi [x,diffs] = jacobi(A,b,x_0,tol) takes a square matrix A and a vector b, a first guess x_0 and a tolerance tol and returns an approximation of the solution of Ax=b using Jacobi iteration. diffs(m) is || x^(m)-x^(m-1)|| make a random first guess for the vector: >> x0 = rand(3,1) x0 = 0.4860 0.8913 0.7621 Jacobi iterate until the tolerance of 10^(-5) is met >> [x,diffs] = jacobi(A,b,x0,10^(-5)); >> format short e >> diffs diffs = Columns 1 through 6 2.0000e-005 1.1940e+000 5.8804e-001 2.3941e-001 1.1502e-001 4.8789e-002 Columns 7 through 12 2.2700e-002 9.8602e-003 4.5064e-003 1.9834e-003 8.9763e-004 3.9793e-004 Columns 13 through 17 1.7913e-004 7.9725e-005 3.5783e-005 1.5961e-005 7.1521e-006 test how close our approximate solution x is to the true solution: >> A*x-b ans = 1.3433e-005 2.8130e-005 3.5134e-005 iterate to smaller stopping tolerance: >> [x,diffs] = jacobi(A,b,x0,10^(-15)); >> diffs diffs = Columns 1 through 6 2.0000e-015 1.1940e+000 5.8804e-001 2.3941e-001 1.1502e-001 4.8789e-002 Columns 7 through 12 2.2700e-002 9.8602e-003 4.5064e-003 1.9834e-003 8.9763e-004 3.9793e-004 Columns 13 through 18 1.7913e-004 7.9725e-005 3.5783e-005 1.5961e-005 7.1521e-006 3.1940e-006 Columns 19 through 24 1.4300e-006 6.3901e-007 2.8595e-007 1.2783e-007 5.7188e-008 2.5570e-008 Columns 25 through 30 1.1438e-008 5.1145e-009 2.2876e-009 1.0230e-009 4.5753e-010 2.0461e-010 Columns 31 through 36 9.1511e-011 4.0925e-011 1.8303e-011 8.1856e-012 3.6608e-012 1.6372e-012 Columns 37 through 42 7.3222e-013 3.2746e-013 1.4645e-013 6.5500e-014 2.9296e-014 1.3101e-014 Columns 43 through 46 5.8634e-015 2.6298e-015 1.1831e-015 5.3083e-016 It took 46 steps but got much closer to the true answer: >> A*x-b ans = -9.9920e-016 -2.0817e-015 -2.5535e-015 Now we'll do it with Gauss Seidel iteration. >> help gauss_seidel [x,diffs] = gauss_seidel(A,b,x_0,tol) takes a square matrix A and a vector b, a first guess x_0 and a tolerance tol and returns an approximation of the solution of Ax=b using Jacobi iteration. diffs(m) is || x^(m)-x^(m-1)|| >> clear diffs >> [x,diffs] = gauss_seidel(A,b,x0,10^(-15)); >> diffs diffs = Columns 1 through 6 2.0000e-015 1.0812e+000 1.8864e-001 3.2436e-002 3.6415e-003 3.3129e-004 Columns 7 through 12 2.8722e-005 1.3981e-006 2.6071e-007 2.8609e-008 2.6241e-009 2.2252e-010 Columns 13 through 18 1.0643e-011 2.0925e-012 2.2452e-013 2.0747e-014 1.7208e-015 8.7604e-017 Gauss Seidel only took 18 steps to reach the same tolerance, where Jacobi took 46 steps.