Anonymous user
Gaussian elimination: Difference between revisions
→version 3: added REXX version 3, changed indentation of the Gauss_elim subroutine.
m (→version 2: added/changed comments and whitespace, increased the number of decimal digits precision, changed highlighting in the output, changed the wording in the REXX sectional header..) |
(→version 3: added REXX version 3, changed indentation of the Gauss_elim subroutine.) |
||
Line 2,278:
exit /*stick a fork in it, we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
Gauss_elim: do j=1 for n; jp=j + 1
do i=jp to n; _=a.j.j / a.i.j
do k=jp to n; a.i.k=a.j.k - _ * a.i.k
end /*k*/
b.i=b.j - _ * b.i
end /*i*/
end /*j*/
x.n=b.n / a.n.n
do j=n-1 to 1 by -1; _=0
do i=j+1 to n; _=_ + a.j.i * x.i
end /*i*/
x.j=(b.j - _) / a.j.j
end /*j*/
say
numeric digits 8
say center('solution', 75, "═");
do o=1 for n; say right('x['o"] = ", 38)
end /*o*/
return</lang>
{{out|output|text= input file: <tt> GAUSS_E.DAT </tt>}}
<pre>
Line 2,349:
x[5] = -0.49098972
x[6] = 0.065760696
</pre>
===version 3===
This is the same as version 2, but in addition, it also shows the residuals.
Code was added to this program version to keep a copy of the original '''A.x.x''' and '''B.n''' arrays (for calculating the residuals).
<br>Also added was rounding the residual numbers to zero if the number of significant decimal digits was
less or equal to 5% of the number of significant fractional decimal digits (in this case, 5% of 1,000 digits for the decimal fraction).
<lang rexx>/*REXX program solves Ax=b with Gaussian elimination and backwards substitution. */
parse arg iFID . /*obtain optional argument from the CL.*/
if iFID=='' | iFID==',' then iFID= 'GAUSS_E.DAT' /*Not specified? Then use the default.*/
pad=left('', 23) /*used for indenting residual numbers. */
do while lines(iFID)\==0 /*read equation sets. */
numeric digits 1000 /*heavy─duty decimal digits precision. */
#=0 /*the number of equations (so far). */
do $=1 while lines(iFID)\==0 /*process the equation*/
z=linein(iFID); if z='' then leave /*Is this a blank line? end─of─data.*/
if $==1 then do; say; say center(' equations ', 75, '▓'); say
end /* [↑] if 1st equation, then show hdr.*/
say z /*display an equation to the terminal. */
if left(space(z), 1)=='*' then iterate /*Is this a comment? Then ignore it.*/
#=# + 1; n=words(z) - 1 /*assign equation #; calculate # items.*/
do e=1 for n; a.#.e= word(z, e); oa.#.e= a.#.e
end /*e*/ /* [↑] process A numbers; save orig.*/
b.#=word(z, n + 1); ob.#=b.# /* ◄── " B " " " */
end /*$*/
if #\==0 then call Gauss_elim /*Not zero? Then display the results. */
say
do i=1 for n; r=0 /*display the residuals to the terminal*/
do j=1 for n; r=r + oa.i.j * x.j /* ┌───◄ don't display a fraction if */
end /*j*/ /* ↓ res ≤ 5% of significant digs.*/
r=format(r - ob.i, , digits() - digits() * 0.05 % 1 , 0) / 1 /*should be tiny*/
say pad 'residual['right(i, length(n) )"] = " left('', r>=0) r /*right justify.*/
end /*i*/
end /*while lines(iFID)¬==0 (1st DO loop)*/
exit /*stick a fork in it, we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
Gauss_elim: do j=1 for n; jp=j + 1
do i=jp to n; _=a.j.j / a.i.j
do k=jp to n; a.i.k=a.j.k - _ * a.i.k
end /*k*/
b.i=b.j - _ * b.i
end /*i*/
end /*j*/
x.n=b.n / a.n.n
do j=n-1 to 1 by -1; _=0
do i=j+1 to n; _=_ + a.j.i * x.i
end /*i*/
x.j=(b.j - _) / a.j.j
end /*j*/ /* [↑] uses backwards substitution. */
say
numeric digits 8 /*for the display, only use 8 digits. */
say center('solution', 75, "═"); say /*a title line for articulated output. */
do o=1 for n; say right('x['o"] = ", 38) left('', x.o>=0) x.o/1
end /*o*/
return</lang>
{{out|output|text= when using the default input file:}}
<pre>
▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓ equations ▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓
* a1 a2 a3 b
* ─── ─── ─── ───
1 2 3 14
2 1 3 13
3 -2 -1 -4
═════════════════════════════════solution══════════════════════════════════
x[1] = 1
x[2] = 2
x[3] = 3
residual[1] = 0
residual[2] = 0
residual[3] = 0
▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓ equations ▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓▓
* a1 a2 a3 a4 a5 a6 b
* ─────── ─────── ─────── ─────── ─────── ─────── ───────
1 0 0 0 0 0 -0.01
1 0.63 0.39 0.25 0.16 0.10 0.61
1 1.26 1.58 1.98 2.49 3.13 0.91
1 1.88 3.55 6.70 12.62 23.80 0.99
1 2.51 6.32 15.88 39.90 100.28 0.60
1 3.14 9.87 31.01 97.41 306.02 0.02
═════════════════════════════════solution══════════════════════════════════
x[1] = -0.01
x[2] = 1.6027904
x[3] = -1.6132031
x[4] = 1.2454941
x[5] = -0.49098972
x[6] = 0.065760696
residual[1] = 0
residual[2] = 0
residual[3] = 0
residual[4] = 0
residual[5] = 0
residual[6] = 0
</pre>
| |||