J: Gauss-Jordan Elimination
The task of coding a method for solving Gauss-Jordan equations of a linear system of algebraic equations is the exercise I have chosen to advance in the learning process. The system is Ax = b, where A is an n-by-n matrix, b and unknown x are n-vectors. First, I started with the simplest form with control structures:
gj0 =: dyad :0 NB. usage is same to %.
y=.y,.b
for_d.i.#y do.
for_r.i.#y do.
if.r=d do.continue.end. NB. do not eliminate d'th row
t=.%/ (<"1(r,d),:(d,d)) { y
for_c.d}.>:i.#y do.
y=.(((<r,c){y)-(t*(<d,c){y)) (<r,c)} y
end.
y=.0 (<r,d)} y NB. ensure zero
end.
end. NB. now A is diagonal but not identity matrix, so:
x=.{:"1 y NB. x = b
for_r.i.#y do.
x=.((r{x)%(<r,r){y) r} x NB. divide by coefficients on diagonal
end.
)
Ab =: (".;._2) 0 :0
0.25 _0.16 _0.38 0.17
0.19 _0.22 _0.02 0.41
0.13 0.08 _0.08 _0.13
0.13 _0.1 _0.32 0.65
)
b =: 0.37 0.01 0.01 1.51
(,.".&.>)('A';'b';'gj0 A,.b';'b %. A')
ββββββββββ¬βββββββββββββββββββββββ
βA β0.25 _0.16 _0.38 0.17β
β β0.19 _0.22 _0.02 0.41β
β β0.13 0.08 _0.08 _0.13β
β β0.13 _0.1 _0.32 0.65β
ββββββββββΌβββββββββββββββββββββββ€
βb β0.37 0.01 0.01 1.51 β
ββββββββββΌβββββββββββββββββββββββ€
βb gj0 A β_1 3 _2 2 β
ββββββββββΌβββββββββββββββββββββββ€
βb %. A β_1 3 _2 2 β
ββββββββββ΄βββββββββββββββββββββββ
Right! Then I decided to get rid of as many control structures as possible:
gj1 =:dyad :0
y=.y,.b
for_d.i.#y do.
for_r.d ({.,]}.~[:>:[) i.#y do. NB. for indices without d
t=.%/ (<"1(r,d),:(d,d)) { y
y=.((r{y)-(t*d{y)) r}y NB. no need to iterate for each column
y=.0 (<r,d)} y
end.
end.
({:"1 y)%(+/}:"1 y) NB. b divide by sum of each item of A (drop zeroes)
)
b gj1 A
_1 3 _2 2
OK, now I can try to translate for_r.
-loop into an unspoken form ... but it looks like it will look more cumbersome and I think I am wrong, but what research without errors? I really want to code the Gauss-Iord method tacitly:
- J-coding exercise
- see if performance is better
- try to understand the code in a few weeks :)
Please help me write it to the end or point out a better approach.
source to share
Thanks to Eelvex who encouraged me to look into addons/math/misc/linear.ijs
, I completed the task with this nice code:
gj=: monad :0
I=. i.#y
for_i. I do. y=. y - (col - i=I) */ (i{y) % i{col=. i{"1 y end.
)
gj Ab
1 0 0 0 _1
0 1 0 0 3
0 0 1 0 _2
0 0 0 1 2
It took a while to understand the verb pivot
in linear.ijs
, but the pencil method helps.
source to share