diff --git a/README.md b/README.md index ca25ff4..ca698c9 100644 --- a/README.md +++ b/README.md @@ -11,7 +11,7 @@ This is a project for the "Laboratorio Computazionale" exam at the University of ## TODO - Parallelization -- Projective coordinates (maybe) +- Homogenization ## Example systems @@ -19,8 +19,8 @@ Here's some tests on 2x2 systems, with the plotted real approximate solutions $$ \begin{align*} -x^2 + y^2 - 4 &= 0 \\ -xy - 1 &= 0 \\ +x^3 + 5x^2 - y - 10 &= 0 \\ +2x^2 - y - 10 &= 0 \\ \end{align*} $$ @@ -32,8 +32,8 @@ $$ $$ \begin{align*} -x^2 + y^2 - 2 &= 0 \\ -xy - 1 &= 0 \\ +x^2 + 2y &= 0 \\ +y - 3x^3 &= 0 \\ \end{align*} $$ @@ -41,15 +41,26 @@ $$ |-------------------|---------------------------------| | ![Solution 2](plots/solutions2.png) | ![Multi-threaded Solution 2](plots/solutions2_8threads.png) | +$$ +\begin{align*} +x^2 + y^2 - 4 &= 0 \\ +xy - 1 &= 0 \\ +\end{align*} +$$ + +| Single-threaded | Multi-threaded (nproc=8) | +|-------------------|---------------------------------| +| ![Solution 3](plots/solutions3.png) | ![Multi-threaded Solution 3](plots/solutions3_8threads.png) | + --- $$ \begin{align*} -x^3 + 5x^2 - y - 10 &= 0 \\ -2x^2 - y - 10 &= 0 \\ +x^2 + y^2 - 2 &= 0 \\ +xy - 1 &= 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) | diff --git a/adapt-step.jl b/adapt-step.jl index 9e5b955..b3a971d 100644 --- a/adapt-step.jl +++ b/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 diff --git a/euler-newton.jl b/euler-newton.jl index 7ec2c04..79dc6c4 100644 --- a/euler-newton.jl +++ b/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 diff --git a/plots/solutions1.png b/plots/solutions1.png index 9837d2c..1bb9e39 100644 Binary files a/plots/solutions1.png and b/plots/solutions1.png differ diff --git a/plots/solutions1_8threads.png b/plots/solutions1_8threads.png deleted file mode 100644 index a70c023..0000000 Binary files a/plots/solutions1_8threads.png and /dev/null differ diff --git a/plots/solutions2.png b/plots/solutions2.png index 639035a..5c83d74 100644 Binary files a/plots/solutions2.png and b/plots/solutions2.png differ diff --git a/plots/solutions2_8threads.png b/plots/solutions2_8threads.png deleted file mode 100644 index e0c8681..0000000 Binary files a/plots/solutions2_8threads.png and /dev/null differ diff --git a/plots/solutions3.png b/plots/solutions3.png index c72d1d4..83e1507 100644 Binary files a/plots/solutions3.png and b/plots/solutions3.png differ diff --git a/plots/solutions3_8threads.png b/plots/solutions3_8threads.png deleted file mode 100644 index a3818ac..0000000 Binary files a/plots/solutions3_8threads.png and /dev/null differ diff --git a/plots/solutions4.png b/plots/solutions4.png new file mode 100644 index 0000000..27d9fad Binary files /dev/null and b/plots/solutions4.png differ diff --git a/report/report.pdf b/report/report.pdf index 3395852..be0ba3c 100644 Binary files a/report/report.pdf and b/report/report.pdf differ diff --git a/report/report.tex b/report/report.tex index 0b405c2..0d9e5a0 100644 --- a/report/report.tex +++ b/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} diff --git a/solve.jl b/solve.jl index 094c71f..f11c4be 100644 --- a/solve.jl +++ b/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")