feat: improve adapt_step function

README.md

@@ -11,12 +11,36 @@ This is a project for the "Laboratorio Computazionale" exam at the University of
 ## TODO
 
 - Parallelization
-- Projective coordinates (maybe)
+- Homogenization
 
 ## Example systems
 
 Here's some tests on 2x2 systems, with the plotted real approximate solutions
 
+$$
+\begin{align*}
+x^3 + 5x^2 - y - 10 &= 0 \\
+2x^2 - y - 10 &= 0 \\
+\end{align*}
+$$
+
+| Single-threaded   | Multi-threaded (nproc=8)        |
+|-------------------|---------------------------------|
+| ![Solution 1](plots/solutions1.png) | ![Multi-threaded Solution 1](plots/solutions1_8threads.png) |
+
+---
+
+$$
+\begin{align*}
+x^2 + 2y  &= 0 \\
+y - 3x^3 &= 0 \\
+\end{align*}
+$$
+
+| Single-threaded   | Multi-threaded (nproc=8)        |
+|-------------------|---------------------------------|
+| ![Solution 2](plots/solutions2.png) | ![Multi-threaded Solution 2](plots/solutions2_8threads.png) |
+
 $$
 \begin{align*}
 x^2 + y^2 - 4 &= 0 \\
@@ -26,7 +50,7 @@ $$
 
 | Single-threaded   | Multi-threaded (nproc=8)        |
 |-------------------|---------------------------------|
-| ![Solution 1](plots/solutions1.png) | ![Multi-threaded Solution 1](plots/solutions1_8threads.png) |
+| ![Solution 3](plots/solutions3.png) | ![Multi-threaded Solution 3](plots/solutions3_8threads.png) |
 
 ---
 
@@ -39,17 +63,4 @@ $$
 
 | Single-threaded   | Multi-threaded (nproc=8)        |
 |-------------------|---------------------------------|
-| ![Solution 2](plots/solutions2.png) | ![Multi-threaded Solution 2](plots/solutions2_8threads.png) |
-
----
-
-$$
-\begin{align*}
-x^3 + 5x^2 - y - 10 &= 0 \\
-2x^2 - y - 10 &= 0 \\
-\end{align*}
-$$
-
-| Single-threaded   | Multi-threaded (nproc=8)        |
-|-------------------|---------------------------------|
-| ![Solution 3](plots/solutions3.png) | ![Multi-threaded Solution 3](plots/solutions3_8threads.png) |
+| ![Solution 4](plots/solutions4.png) | ![Multi-threaded Solution 4](plots/solutions4_8threads.png) |

adapt-step.jl

@@ -1,17 +1,18 @@
 module AdaptStep
   using LinearAlgebra
+  using TypedPolynomials
 
   export adapt_step
 
   # Adaptive step size
-  function adapt_step(x, x_old, step, m)
-    Δ = LinearAlgebra.norm(x - x_old)
-    if Δ > 0.1
+  function adapt_step(H, x, t, step, m)
+    Δ = LinearAlgebra.norm([h(variables(H(t))=>x) for h in H(t-step)])
+    if Δ > 1e-8
       step = 0.5 * step
       m = 0
     else
       m+=1
-      if (m == 5)
+      if (m == 5) && (step < 0.05)
         step = 2 * step
         m = 0
       end

euler-newton.jl

@@ -16,14 +16,11 @@ module EulerNewton
     xh = x + Δx * step_size
 
     # Corrector step
-    JHh=differentiate(H(t+step_size), vars)
-    for _ in 1:10
+    JHh=differentiate(H(t-step_size), vars)
+    for _ in 1:5
       JH = [jh(vars=>xh) for jh in JHh]
-      Δx = JH \ -[h(vars=>xh) for h in H(t+step_size)]
+      Δx = JH \ -[h(vars=>xh) for h in H(t-step_size)]
       xh = xh + Δx
-      if LinearAlgebra.norm([h(vars=>xh) for h in H(t+step_size)]) < 1e-8
-        break
-      end
     end
 
     return xh

report/report.tex

@@ -179,19 +179,23 @@ by the sequence of points $(z_i)_{i\in\N}$ defined by the recurrence relation
 $$ z_{i+1}=z_i+h\cdot f(z_i,t_i) ,$$
 where $h$ is the step size.
 In our case, we have $$f(z,t)=-\left(\frac{\partial H}{\partial z}(z,t)\right)^{-1}\frac{\partial H}{\partial t}(z,t)$$ and $t_0=1$, since we track from $1$ to $0$. For the same
-reason, we set $$t_i=t_{i-1}-h.$$ We will also use a variable step size, based on the output of each iteration.
+reason, we set $$t_{i+1}=t_i-h.$$
 \subsubsection{Corrector: Newton's method}
 Since we want to solve $$H(z,t)=0,$$ we can use Newton's method to improve the approximation $\widetilde{z_i}$ obtained by Euler's method to a solution of such equation.
 This is done by moving towards the root of the tangent line of $H$ at the current approximation, or in other words through the iteration
-$$ z_{i+1}=z_i-\left(\frac{\partial H}{\partial z}(z_i,t_i)\right)^{-1}H(z_i,t_i) ,$$
-where this time $z_0=\widetilde{z}_i$ and $t_0=t_i$ as obtained in the Euler step.
+$$ z_{i+1}=z_i-\left(\frac{\partial H}{\partial z}(z_i,t_{i+1})\right)^{-1}H(z_i,t_{i+1}) ,$$
+where this time $z_0=\widetilde{z}_i$, with $\widetilde{z}_i$ and $t_{i+1}$ obtained from the $i$-th Euler step.
 
-Usually, only a few steps of Newton's method are needed; we will use a fixed maximum of $10$ steps,
-stopping the iterations when the desired accuracy is reached, for instance when the norm of $H(z_i,t_i)$ is less than $10^{-8}$.
+Usually, only a few steps of Newton's method are needed; we will use a fixed number of 5 iterations.
 At this point, we use the final value of the Newton iteration as the starting value for the next Euler step.
 \subsubsection{Adaptive step size}
+In order to improve the efficiency of the algorithm, we will use an adaptive step size, which will be based on the norm of the residual of the Newton iteration.
+If the desired accuracy is not reached, for instance when the norm of $H(z_i,t_i)$ is bigger than $10^{-8}$,
+then we halve the step size; if instead we have 5 "successful" iterations in a row, we double the step size.
 
 \section{Parallelization}
+\subsection{Multithreading}
+\subsection{MPI}
 
 \section{Implementation}
 \subsection{Julia code}

solve.jl

@@ -4,24 +4,25 @@ using TypedPolynomials
 # Local dependencies
 include("start-system.jl")
 include("homotopy.jl")
-include("homogenize.jl")
+#  include("homogenize.jl")
 include("euler-newton.jl")
 include("adapt-step.jl")
 include("plot.jl")
 using .StartSystem
 using .Homotopy
-using .Homogenize
+#  using .Homogenize
 using .EulerNewton
 using .AdaptStep
 using .Plot
 
 # Main homotopy continuation loop
-function solve(F, (G, roots) = start_system(F), maxsteps=10000)
+function solve(F, (G, roots) = start_system(F), maxsteps = 1000)
   #  F=homogenize(F)
   H=homotopy(F,G)
   solutions = []
+  step_array = []
 
-  @time Threads.@threads for r in roots
+  Threads.@threads for r in roots
     t = 1.0
     step_size = 0.01
     x0 = r
@@ -29,30 +30,43 @@ function solve(F, (G, roots) = start_system(F), maxsteps=10000)
     steps = 0
 
     while t > 0 && steps < maxsteps
-      x = en_step(H, x0, t, step_size)
-      (m, step_size) = adapt_step(x, x0, step_size, m)
-      x0 = x
+      x0 = en_step(H, x0, t, step_size)
+      (m, step_size) = adapt_step(H, x0, t, step_size, m)
       t -= step_size
       steps += 1
     end
     push!(solutions, x0)
+    push!(step_array, steps)
   end
 
-  return solutions
+  return (solutions, step_array)
 end
 
 # Input polynomial system
 @polyvar x y
+C = [x^3 - y + 5x^2 - 10, 2x^2 - y - 10]
+Q = [x^2 + 2y, y - 3x^3]
 F = [x*y - 1, x^2 + y^2 - 4]
 T = [x*y - 1, x^2 + y^2 - 2]
-C = [x^3 - y + 5x^2 - 10, 2x^2 - y - 10]
 
-sF = filter(u -> imag(u[1]) < 0.1 && imag(u[2]) < 0.1, solve(F))
-sT = filter(u -> imag(u[1]) < 0.1 && imag(u[2]) < 0.1, solve(T))
-sC = filter(u -> imag(u[1]) < 0.1 && imag(u[2]) < 0.1, solve(C))
+(sC, stepsC) = solve(C)
+(sQ, stepsQ) = solve(Q)
+(sF, stepsF) = solve(F)
+(sT, stepsT) = solve(T)
+
+println("C: ", stepsC)
+println("Q: ", stepsQ)
+println("F: ", stepsF)
+println("T: ", stepsT)
+
+sC = filter(u -> imag(u[1]) < 0.1 && imag(u[2]) < 0.1, sC)
+sQ = filter(u -> imag(u[1]) < 0.1 && imag(u[2]) < 0.1, sQ)
+sF = filter(u -> imag(u[1]) < 0.1 && imag(u[2]) < 0.1, sF)
+sT = filter(u -> imag(u[1]) < 0.1 && imag(u[2]) < 0.1, sT)
 
 # Plotting the system and the real solutions
 ENV["GKSwstype"]="nul"
-plot_real(sF, F, 4, 4, "1")
-plot_real(sT, T, 4, 4, "2")
-plot_real(sC, C, 6, 12, "3")
+plot_real(sC, C, 6, 12, "1")
+plot_real(sQ, Q, 2, 2, "2")
+plot_real(sF, F, 4, 4, "3")
+plot_real(sT, T, 4, 4, "4")