# Solve a matrix equation by Jacobi's method (Revised)

Mathematical Background

Let A be an N by N matrix of real numbers, b a vector of N real numbers and x a vector N unknown real numbers. A matrix equation is Ax = b.

Jacobi's method is as follows: decompose A = D + R, where D is the matrix of diagonals, and R is the remaining entries.

if you make an initial guess solution x0, an improved solution is x1 = inverse(D) * (b - Rx) where all multiplications are matrix-vector multiplication and inverse(D) is the matrix inverse.

Problem Specification

• Input: Your complete program should accept as input the following data: the matrix A, the vector b, an initial guess x0, and an 'error' number e.
• Output: The program must output the minimum number of iterations such that the latest solution differs by the true solution, by at most e. This means each component of the vectors in absolute magnitude differ by at most e. You must use Jacobi's method for iterations.

How the data is inputted is your choice; it could be your own syntax on a command line, you could take input from a file, whatever you choose.

How the data is outputted is your choice; it can be written to a file, displayed in the command line, written as ASCII art, anything, as long as it is readable by a human.

Further Details

You are not given the true solution: how you calculate the true solution is up to you entirely. You may solve it by Cramer's rule for example, or computing an inverse directly. What matters is that you have a true solution to be able to compare to iterations.

Precision is an issue; some people's 'exact solutions' for comparison may differ. For the purposes of this code golf the exact solution must be true to 10 decimal places.

To be absolutely clear, if even one component of your present iteration solution exceeds its corresponding component in the true solution by e, then you need to keep iterating.

The upper limit to N varies based on what hardware you're using and how much time you're willing to spend running the program. For the purposes of this code golf, assume maximum N = 50.

Preconditions

When your program is called, you are free to assume that the following holds at all times:

• N > 1 and N < 51, i.e. you will never be given a scalar equation, always a matrix equation.
• All inputs are over the field of real numbers, and will never be complex.
• The matrix A is always such that the method converges to the true solution, such that you can always find a number of iterations to minimise the error (as defined above) below or equal to e.
• A is never the zero matrix or the identity matrix, i.e. there is one solution.

Test Cases

A = ((9, -2), (1, 3)), b = (3,4), x0 = (1,1), e = 0.04


The true solution is (0.586, 1.138). The first iteration gives x1 = (5/9, 1), differing by more than 0.04 from the true solution, by at least one component. Taking another iteration we find, x2 = (0.555, 1.148) which differs by less than 0.04 from (0.586, 1.138). Thus the output is

2


A = ((2, 3), (1, 4)), b = (2, -1), x0 = (2.7, -0.7), e = 1.0


In this case the true solution is (2.2, -0.8) and the initial guess x0 already has error less than e = 1.0, thus we output 0. That is, whenever you do not need to make an iteration, you simply output

0


Submission Assessment

This is code golf, with all standard loopholes hereby disallowed. The shortest correct complete program (or function), i.e. lowest number of bytes wins. It is discouraged to use things like Mathematica which wrap up a lot of the necessary steps into one function, but use any language you want.

• You really should wait to get more feedback for it, especially given the recent closed post. PPCG challenges usually share common structure in specifications, what usually contribute to them being easy to understand, rather than tiresome and ambiguous. Try to look at some of the reasonably upvoted challenges and mimic the format. Nov 15, 2017 at 15:05
• @Uriel I realise this, but I feel I've been exhaustive in my specification, and the format while not exactly fitting the majority of questions, can be read linearly, and guide the reader accordingly. The format should also keep the content of the problem itself in mind. Nov 15, 2017 at 15:06
• "The shortest correct complete program" sounds like you only allow programs and not functions. I'd add "/function".
Nov 15, 2017 at 15:20
• +1 formatting makes or breaks my brain's ability to concentrate on a question Nov 15, 2017 at 15:22
• @user1997744 Yup, makes sense. I believe the default is that any other code, like other functions or python imports, is allowed, but also included in the bytecount. Nov 15, 2017 at 15:42

# APL (Dyalog), 786865 49 bytes

Exactly the type of problem APL was created for.

-3 thanks to Erik the Outgolfer. -11 thanks to ngn.

Anonymous infix function. Takes A b e as left argument and x as right argument. Prints result to STDOUT as vertical unary using 1 as tally marks, followed by 0 as punctuation. This means that even a 0-result can be seen, being no 1s before the 0.

{⎕←∨/e<|⍵-b⌹⊃A b e←⍺:⍺∇D+.×b-⍵+.×⍨A-⌹D←⌹A×=/¨⍳⍴A}


Try it online!

Notice how the code reads very similarly to the problem specification:

{} on the given A, b, and e, and and the given x,
⎕← print
∨/ whether there is any truth in the statement that
e< e is smaller than
|⍵- the absolute value of x minus
b⌹ b matrix-divided by
⊃A b e the first of A, b, and e (i.e. A)
←⍺ which are the left argument
: and if so,
⍺∇ recurse on
D+.× D matrix-times
b- b minus
⍵+.×⍨ x, matrix multiplied by
A- A minus
⌹D the inverse of D (the remaining entries)
← where D is
A× A where
=/¨ there are equal
⍳ coordinates for
⍴A the shape of A (i.e. the diagonal)

### Step-by-step explanation

The actual order of execution right-to-left:

{} anonymous function where ⍺ is A b e and ⍵ is x:
A b c←⍺ split left argument into A, b, and e
⊃ pick the first (A)
b⌹ matrix division with b (gives true value of x)
⍵- differences between current values of x and those
| absolute values
e< acceptable error less than those?
∨/ true for any? (lit. OR reduction)
⎕← print that Boolean to STDOUT
: and if so:
⍴A shape of A
⍳ matrix of that shape where each cell is its own coordinates
=/¨ for each cell, is the vertical and horizontal coordinates equal? (diagonal)
A× multiply the cells of A with the that (extracts diagonal)
⌹ matrix inverse
D← store in D (for Diagonal)
⌹ inverse (back to normal)
A- subtract from A
⍵+.×⍨ matrix multiply (same thing as dot product, hence the .) that with x
b- subtract that from b
D+.× matrix product of D and that
⍺∇ apply this function with given A b e and that as new value of x

• The output should be the number of iterations required for an accuracy of e. Nov 15, 2017 at 15:21
• -1: This is not a solution. You need x0 since the whole point is to know how many steps it takes to get to a desired accuracy from a particular starting point. Nov 15, 2017 at 15:22
• @user1997744 Ah, I misunderstood the problem. Sorry.
Nov 15, 2017 at 15:23
• @user1997744 Better?
Nov 15, 2017 at 16:03
• @user1997744 Not arithmetic operation, just ability to read unary, where indeed 0 is nothing.
Nov 15, 2017 at 16:43

# Python 3, 132 bytes

f=lambda A,b,x,e:e<l.norm(x-dot(l.inv(A),b))and 1+f(A,b,dot(l.inv(d(d(A))),b-dot(A-d(d(A)),x)),e)
from numpy import*
l=linalg
d=diag


Try it online!

Uses a recursive solution.

• @Adám I'm not sure I quite understand. I interpreted it as f not having the name within the code block, which I have now fixed; however, if it is a different issue entirely, it may still be a problem. Nov 15, 2017 at 22:50
• @Adám That answer seems to corroborate what I currently have; it is a function with helper code that is capable of working as a unit after its definition. Nov 15, 2017 at 23:40
• Ah, ok. Never mind then. I don't know Python, so I was just curious. Good job!
Nov 16, 2017 at 7:18
• Isn't the stopping criterion "This means each component of the vectors in absolute magnitude differ by at most e"? Basically the max-norm, not L2-norm. Nov 16, 2017 at 19:11
• @NikoNyrh Fixed. Nov 16, 2017 at 19:52

# R, 138 bytes

function(A,x,b,e){R=A-(D=diag(diag(A)))
g=solve(A,b)
if(norm(t(x-g),"M")<e)T=0
while(norm((y=solve(D)%*%(b-R%*%x))-g,"M")>e){T=T+1;x=y}
T}


Try it online!

thanks to NikoNyrh for fixing a bug

It's also worth noting that there is an R package, Rlinsolve that contains a lsolve.jacobi function, returning a list with x (the solution) and iter (the number of iterations required), but I'm not sure that it does the correct computations.

• Isn't the stopping criterion "This means each component of the vectors in absolute magnitude differ by at most e"? Basically the max-norm, not L2-norm. Nov 16, 2017 at 19:11
• @NikoNyrh you are correct! fortunately, the norm function provides that for me as well for no additional bytes. Nov 16, 2017 at 19:38

## Clojure, 212198 196 bytes

#(let[E(range)I(iterate(fn[X](map(fn[r b d](/(- b(apply +(map * r X)))d))(map assoc % E(repeat 0))%2(map nth % E)))%3)](count(for[x I :while(not-every?(fn[e](<(- %4)e %4))(map -(nth I 1e9)x))]x)))


Implemented without a matrix library, it iterates the process 1e9 times to get the correct answer. This wouldn't work on too ill-conditioned inputs but should work fine in practice.

Less golfed, I was quite happy with the expressions of R and D :) The first input % (A) has to be a vector, not a list so that assoc can be used.

(def f #(let[R(map assoc %(range)(repeat 0))
D(map nth %(range))
I(iterate(fn[X](map(fn[r x b d](/(- b(apply +(map * r x)))d))R(repeat X)%2 D))%3)]
(->> I
(take-while (fn[x](not-every?(fn[e](<(- %4)e %4))(map -(nth I 1e9)x))))
count)))