Commit 8dd4ef8e authored by Lambert Theisen's avatar Lambert Theisen 🔥

Add some files to two region exmaple

parent 9860c0a2
import numpy
dofs = 1425
kn_range = [(i + 1) / 10 for i in range(10)]
numknudsen = len(kn_range)
kn_range1 = [1.0]
numknudsen1 = len(kn_range1)
kn_range2 = [(i + 1) / 10 for i in range(10)]
numknudsen2 = len(kn_range2)
fields = [
"theta", "s_x", "s_y", "p", "u_x", "u_x", "sigma_xx", "sigma_xy", "sigma_yx"
]
numfields = len(fields)
V = numpy.zeros(shape=(numknudsen**2, numfields * dofs))
V = numpy.zeros(shape=(numknudsen1 * numknudsen2, numfields * dofs))
for iKn, kn1 in enumerate(kn_range):
for jKn, kn2 in enumerate(kn_range):
for iKn, kn1 in enumerate(kn_range1):
for jKn, kn2 in enumerate(kn_range2):
# row = numpy.zeros(shape=(numfields * dofs))
field_data = numpy.zeros(shape=(numfields, dofs))
......@@ -25,8 +27,8 @@ for iKn, kn1 in enumerate(kn_range):
# print("next")
row = field_data.reshape(numfields * dofs)
# print(row)
print(iKn * numknudsen + jKn)
V[iKn * numknudsen + jKn] = row
print(iKn * numknudsen1 + jKn)
V[iKn * numknudsen1 + jKn] = row
numpy.savetxt("V_{}.mat".format(numknudsen**2), V)
numpy.savetxt("V_{}.mat".format(numknudsen1 * numknudsen2), V)
print(V)
using DelimitedFiles
using LinearAlgebra
using Plots
using SymPy
kn_range1 = 0.1
kn_range2 = 0.1:0.1:1.0
V10 = readdlm("V_10.mat", ' ', Float64, '\n');
F10 = svd(V10);
F10SN = normalize(F10.S, Inf)
# low rank approximation with rank R
R = 5
V10LR = F10.U[:,1:R] * diagm(F10.S[1:R]) * F10.V[:,1:R]'
# sanity check: Frobenius norm should be S[R+1] is not the case but still low
# ! FIXME
println(norm(V10 - V10LR ,2))
L = F10.U[:,1:R]
OPCompl = I(10) - L * L' # projector to complement of LL'
# ansatz
p(t1, t2, deg) = t2^deg # monomials
# p(t1, t2, deg) = (1/t2)^deg # 1/x
# p(t1, t2, deg) = exp(t2*deg) # exponentials
# p(t1, t2, deg) = log(t1)^deg # logrithms
# p(t1, t2, deg) = (t1*t2)^deg # logrithms
# p(t1, t2, deg) = sin(deg * pi * ((t2-0.1)*100+t1*10) / 100) # logrithms
# p(t1, t2, deg) = (t2-0.1)*100+t1*10 # logrithms
# p(t1, t2, deg) = 1/(t1+t2)
# p(t1, t2, deg) = t2 + 1/(t1)^deg
# p(t1, t2, deg) = 1/(t2)^deg + 1/(t1)^deg
# p(t1, t2, deg) = 1 + t1 + t2
# function p(t1, t2, deg)
# x = Sym("x")
# f(x) = x^float(deg-1)
# G(y) = subs(integrate(f(x), x), x, y)
# return N(G(float(t1)))
# end
degrange = -1:1
P = vcat(
[[p(kn1, kn2, d) for d in degrange]' for kn1 in kn_range1, kn2 in kn_range2]...,
)
# P = V10LR[:,1:30] # optimal
# P = F10.U[:,1:1] # optimal
# P = hcat(F10.U[:,1:1],F10.U[:,1:1]) # linearly dependent
# P = F10.V[1:20,:]' # fails
# P = (L * L')[:,1:20] # fails
Q = Matrix(qr(P).Q)
objective = OPCompl * Q # --> MIN!!!
display(plot([svd(OPCompl).S, svd(objective).S], yaxis=:log, title="Singular Values, max ansatz degree, R=$R", xlabel="index", ylabel="value", label=["(I - LL')" "(I - LL') * Q"], marker = (:dot)))
if false # comparison
V361 = readdlm("V_361.mat", ' ', Float64, '\n');
F361 = svd(V361);
F361SN = normalize(F361.S, Inf)
V400 = readdlm("V_400.mat", ' ', Float64, '\n');
F400 = svd(V400);
F400SN = normalize(F400.S, Inf)
max = 100
plot([F10SN, F400SN[1:max]], yaxis=:log, title="Normlized Singular Values", xlabel="index", ylabel="value", label=["h=0.1" "h=0.05"], marker = (:dot))
end
......@@ -30,7 +30,8 @@ OPCompl = I(100) - L * L' # projector to complement of LL'
# p(t1, t2, deg) = (t2-0.1)*100+t1*10 # logrithms
# p(t1, t2, deg) = 1/(t1+t2)
# p(t1, t2, deg) = t2 + 1/(t1)^deg
p(t1, t2, deg) = 1/(t2)^deg + 1/(t1)^deg
# p(t1, t2, deg) = 1/(t2)^deg + 1/(t1)^deg
p(t1, t2, deg) = 1 + t1 + t2
# function p(t1, t2, deg)
# x = Sym("x")
# f(x) = x^float(deg-1)
......@@ -42,6 +43,11 @@ degrange = -3:3
P = vcat(
[[p(kn1, kn2, d) for d in degrange]' for kn1 in kn_range, kn2 in kn_range]...,
)
# P = V100LR[:,1:30] # optimal
P = F100.U[:,1:5] # optimal
# P = hcat(F100.U[:,1:1],F100.U[:,1:1]) # linearly dependent
# P = F100.V[1:20,:]' # fails
# P = (L * L')[:,1:20] # fails
Q = Matrix(qr(P).Q)
objective = OPCompl * Q # --> MIN!!!
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment