diff --git a/l2/3.jl b/l2/3.jl new file mode 100755 index 0000000..7f8e185 --- /dev/null +++ b/l2/3.jl @@ -0,0 +1,43 @@ +#!/usr/bin/env julia + +# Jacek Poziemski 272389 + +using LinearAlgebra + +include("hilb.jl") +include("matcond.jl") + +struct Errors + gauss::Float64 + inverse::Float64 +end + +function solve(A::Matrix{Float64}, n::Int)::Errors + x::Vector{Float64} = ones(Float64, n) + b::Vector{Float64} = A * x + gauss::Float64 = norm((A \ b) - x) / norm(x) + inverse::Float64 = norm((inv(A) * b) - x) / norm(x) + + return Errors(gauss, inverse) +end + +println("hilbert matrices") +println() + +for n in 1:30 + A::Matrix{Float64} = hilb(n) + errors = solve(A, n) + println("n: $n; cond: $(cond(A)); rank: $(rank(A)); gauss error: $(errors.gauss); inverse error: $(errors.inverse)") +end + +println() +println("random matrices") +println() + +for n in [5, 10, 20] + for c in [0, 1, 3, 7, 12, 16] + A::Matrix{Float64} = matcond(n, 10.0^c) + errors = solve(A, n) + println("n: $n; c: $c; cond: $(cond(A)); rank: $(rank(A)); gauss error: $(errors.gauss); inverse error: $(errors.inverse)") + end +end diff --git a/l2/hilb.jl b/l2/hilb.jl new file mode 100644 index 0000000..f3e6f29 --- /dev/null +++ b/l2/hilb.jl @@ -0,0 +1,16 @@ +# https://cs.pwr.edu.pl/zielinski/lectures/scna/hilb.jl +function hilb(n::Int) +# Function generates the Hilbert matrix A of size n, +# A (i, j) = 1 / (i + j - 1) +# Inputs: +# n: size of matrix A, n>=1 +# +# +# Usage: hilb(10) +# +# Pawel Zielinski + if n < 1 + error("size n should be >= 1") + end + return [1 / (i + j - 1) for i in 1:n, j in 1:n] +end diff --git a/l2/matcond.jl b/l2/matcond.jl new file mode 100644 index 0000000..5ae4218 --- /dev/null +++ b/l2/matcond.jl @@ -0,0 +1,22 @@ +# https://cs.pwr.edu.pl/zielinski/lectures/scna/matcond.jl +using LinearAlgebra + +function matcond(n::Int, c::Float64) +# Function generates a random square matrix A of size n with +# a given condition number c. +# Inputs: +# n: size of matrix A, n>1 +# c: condition of matrix A, c>= 1.0 +# +# Usage: matcond(10, 100.0) +# +# Pawel Zielinski + if n < 2 + error("size n should be > 1") + end + if c< 1.0 + error("condition number c of a matrix should be >= 1.0") + end + (U,S,V)=svd(rand(n,n)) + return U*diagm(0 =>[LinRange(1.0,c,n);])*V' +end