Anonymous user
Gaussian elimination: Difference between revisions
→{{header|Modula-3}}
m (C++ - avoid division by zero) |
|||
Line 3,283:
</pre>
=={{header|Nim}}==
{{trans|Kotlin}}
<lang Nim>const Eps = 1e-14 # Tolerance required.
type
Vector[N: static Positive] = array[N, float]
Matrix[M, N: static Positive] = array[M, Vector[N]]
SquareMatrix[N: static Positive] = Matrix[N, N]
func gaussPartialScaled(a: SquareMatrix; b: Vector): Vector =
doAssert a.N == b.N, "matrix and vector have incompatible dimensions"
const N = a.N
var m: Matrix[N, N + 1]
for i, row in a:
m[i][0..<N] = row
m[i][N] = b[i]
for k in 0..<N:
var imax = 0
var vmax = -1.0
for i in k..<N:
# Compute scale factor s = max abs in row.
var s = -1.0
for j in k..N:
let e = abs(m[i][j])
if e > s: s = e
# Scale the abs used to pick the pivot.
let val = abs(m[i][k]) / s
if val > vmax:
imax = i
vmax = val
if m[imax][k] == 0:
raise newException(ValueError, "matrix is singular")
swap m[imax], m[k]
for i in (k + 1)..<N:
for j in (k + 1)..N:
m[i][j] -= m[k][j] * m[i][k] / m[k][k]
m[i][k] = 0
for i in countdown(N - 1, 0):
result[i] = m[i][N]
for j in (i + 1)..<N:
result[i] -= m[i][j] * result[j]
result[i] /= m[i][i]
#———————————————————————————————————————————————————————————————————————————————————————————————————
let a: SquareMatrix[6] = [[1.00, 0.00, 0.00, 0.00, 0.00, 0.00],
[1.00, 0.63, 0.39, 0.25, 0.16, 0.10],
[1.00, 1.26, 1.58, 1.98, 2.49, 3.13],
[1.00, 1.88, 3.55, 6.70, 12.62, 23.80],
[1.00, 2.51, 6.32, 15.88, 39.90, 100.28],
[1.00, 3.14, 9.87, 31.01, 97.41, 306.02]]
let b: Vector[6] = [-0.01, 0.61, 0.91, 0.99, 0.60, 0.02]
let refx: Vector[6] = [-0.01, 1.602790394502114, -1.6132030599055613,
1.2454941213714368, -0.4909897195846576, 0.065760696175232]
let x = gaussPartialScaled(a, b)
echo x
for i, xi in x:
if abs(xi - refx[i]) > Eps:
echo "Out of tolerance."
echo "Expected values are ", refx
break</lang>
{{out}}
<pre>[-0.01, 1.602790394502114, -1.613203059905562, 1.245494121371439, -0.4909897195846595, 0.06576069617523238]</pre>
=={{header|OCaml}}==
| |||