77 lines
2.1 KiB
Julia
77 lines
2.1 KiB
Julia
include("Engine.jl")
|
|
|
|
using SparseArrays
|
|
|
|
# this problem is from a sangaku by Irisawa Shintarō Hiroatsu. the article below
|
|
# includes a nice translation of the problem statement, which was recorded in
|
|
# Uchida Itsumi's book _Kokon sankan_ (_Mathematics, Past and Present_)
|
|
#
|
|
# "Japan's 'Wasan' Mathematical Tradition", by Abe Haruki
|
|
# https://www.nippon.com/en/japan-topics/c12801/
|
|
#
|
|
|
|
# initialize the partial gram matrix
|
|
J = Int64[]
|
|
K = Int64[]
|
|
values = BigFloat[]
|
|
for s in 1:9
|
|
# each sphere is represented by a spacelike vector
|
|
push!(J, s)
|
|
push!(K, s)
|
|
push!(values, 1)
|
|
|
|
# the circumscribing sphere is internally tangent to all of the other spheres
|
|
if s > 1
|
|
append!(J, [1, s])
|
|
append!(K, [s, 1])
|
|
append!(values, [1, 1])
|
|
end
|
|
|
|
if s > 3
|
|
# each chain sphere is externally tangent to the two nucleus spheres
|
|
for n in 2:3
|
|
append!(J, [s, n])
|
|
append!(K, [n, s])
|
|
append!(values, [-1, -1])
|
|
end
|
|
|
|
# each chain sphere is externally tangent to the next sphere in the chain
|
|
s_next = 4 + mod(s-3, 6)
|
|
append!(J, [s, s_next])
|
|
append!(K, [s_next, s])
|
|
append!(values, [-1, -1])
|
|
end
|
|
end
|
|
gram = sparse(J, K, values)
|
|
|
|
# make an initial guess
|
|
guess = hcat(
|
|
Engine.sphere(BigFloat[0, 0, 0], BigFloat(15)),
|
|
Engine.sphere(BigFloat[0, 0, -9], BigFloat(5)),
|
|
Engine.sphere(BigFloat[0, 0, 11], BigFloat(3)),
|
|
(
|
|
Engine.sphere(9*BigFloat[cos(k*π/3), sin(k*π/3), 0], BigFloat(2.5))
|
|
for k in 1:6
|
|
)...
|
|
)
|
|
frozen = [CartesianIndex(4, k) for k in 1:4]
|
|
|
|
# complete the gram matrix using Newton's method with backtracking
|
|
L, success, history = Engine.realize_gram(gram, guess, frozen)
|
|
completed_gram = L'*Engine.Q*L
|
|
println("Completed Gram matrix:\n")
|
|
display(completed_gram)
|
|
if success
|
|
println("\nTarget accuracy achieved!")
|
|
else
|
|
println("\nFailed to reach target accuracy")
|
|
end
|
|
println("Steps: ", size(history.scaled_loss, 1))
|
|
println("Loss: ", history.scaled_loss[end], "\n")
|
|
if success
|
|
println("Chain diameters:")
|
|
println(" ", 1 / L[4,4], " sun (given)")
|
|
for k in 5:9
|
|
println(" ", 1 / L[4,k], " sun")
|
|
end
|
|
end |