Skip to content

Migrate to NonlinearSolve.jl for nonlinear system solving functionalities #126

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
wants to merge 5 commits into
base: master
Choose a base branch
from
+81 −137
Draft

Migrate to NonlinearSolve.jl for nonlinear system solving functionalities #126

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
b472db7
Initial plan
Copilot Aug 26, 2025
9ab058b
Add NonlinearSolve.jl dependency and plan integration
Copilot Aug 26, 2025
6e7104d
Replace manual Newton iterations with NonlinearSolve.jl in Trapezoid …
Copilot Aug 26, 2025
b6bc30f
Add NonlinearSolve.jl to more methods - Trapezoid working, others nee…
Copilot Aug 26, 2025
2443832
Successfully integrate NonlinearSolve.jl - Trapezoid method fully wor…
Copilot Aug 26, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 3 additions & 1 deletion Project.toml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,7 @@ ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
HypergeometricFunctions = "34004b35-14d8-5ef3-9330-4cdb6864b03a"
InvertedIndices = "41ab1584-1d38-5bbf-9106-f11c6c58b48f"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
Polynomials = "f27b6e38-b328-58d1-80ce-0feddd5e7a45"
RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
Expand All @@ -31,12 +32,13 @@ FractionalDiffEqFdeSolverExt = "FdeSolver"
[compat]
ConcreteStructs = "0.2.2"
DiffEqBase = "6.156"
FFTW = "1.8.0"
FastClosures = "0.3.2"
FdeSolver = "1.0.9"
FFTW = "1.8.0"
ForwardDiff = "0.10.36"
HypergeometricFunctions = "0.3.23"
InvertedIndices = "1.3"
NonlinearSolve = "4.10.0"
Polynomials = "3, 4"
RecipesBase = "1.3.4"
RecursiveArrayTools = "3.27"
Expand Down
2 changes: 1 addition & 1 deletion src/FractionalDiffEq.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
module FractionalDiffEq

using LinearAlgebra, Reexport, SciMLBase, SpecialFunctions, SparseArrays, ToeplitzMatrices, RecursiveArrayTools,
FFTW, ForwardDiff, Polynomials, HypergeometricFunctions, DiffEqBase, ConcreteStructs, FastClosures
FFTW, ForwardDiff, Polynomials, HypergeometricFunctions, DiffEqBase, ConcreteStructs, FastClosures, NonlinearSolve

import SciMLBase: __solve
import DiffEqBase: solve
Expand Down
83 changes: 31 additions & 52 deletions src/fode/bdf.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -212,35 +212,22 @@ function BDF_triangolo(
end
Phi_n = St + halpha * (zn[:, n1] + Phi)
yn0 = cache.y[n]
fn0 = similar(yn0)
prob.f(fn0, yn0, p, mesh[n1])
@views jac(Jfn0, yn0, p, mesh[n1])
Gn0 = yn0 - halpha * omega[1] * fn0 - Phi_n
stop = false
it = 0
yn1 = similar(yn0)
fn1 = similar(yn0)
while !stop
JGn0 = zeros(T, problem_size, problem_size) + I - halpha * omega[1] * Jfn0
yn1 = yn0 - vec(JGn0 \ Gn0)
prob.f(fn1, yn1, p, mesh[n1])
Gn1 = yn1 - halpha * omega[1] * fn1 - Phi_n
it = it + 1

stop = (norm(yn1 - yn0, Inf) < abstol) || (norm(Gn1, Inf) < abstol)
if it > maxiters && !stop
@warn "Non Convergence"
stop = true
end

yn0 = yn1
Gn0 = Gn1
if !stop
@views jac(Jfn0, yn0, p, mesh[n1])
end

# Replace manual Newton iteration with NonlinearSolve.jl
# Solve: G(y) = y - halpha * omega[1] * f(y) - Phi_n = 0
function bdf_nlprob_f!(G, y, p_nl)
temp = similar(y)
prob.f(temp, y, p, mesh[n1])
G .= y - halpha * omega[1] * temp - Phi_n
end
cache.y[n1] = yn1
cache.fy[n1] = fn1

nlprob = NonlinearProblem(bdf_nlprob_f!, yn0)
nlsol = solve(nlprob, NewtonRaphson(); reltol=cache.reltol, abstol=cache.abstol, maxiters=cache.maxiters)

cache.y[n1] = nlsol.u
temp = zeros(length(nlsol.u))
prob.f(temp, nlsol.u, p, mesh[n1])
cache.fy[n1] = temp
end
end

Expand Down Expand Up @@ -283,37 +270,29 @@ function BDF_first_approximations(cache::BDFCache{iip, T}) where {iip, T}
end
JF[((j - 1) * problem_size + 1):(j * problem_size), ((j - 1) * problem_size + 1):(j * problem_size)] .= J_temp
end
stop = false
it = 0
F1 = similar(F0)
Y1 = similar(Y0)
while !stop
JG = Ims - W * JF
recursive_unflatten!(Y1, vec(Y0) - JG \ G0)


# Replace manual Newton iteration with NonlinearSolve.jl
# Solve: G(Y) = Y - B0 - W * F(Y) = 0
function bdf_first_nlprob_f!(G, Y_vec, p_nl)
recursive_unflatten!(Y1, Y_vec)
for j in 1:s
prob.f(F1.u[j], Y1.u[j], p, mesh[j + 1])
end
G1 = vec(Y1 - B0) - W * vec(F1)

it = it + 1
stop = (norm(Y1 - Y0, Inf) < abstol) || (norm(G1, Inf) < abstol)
if it > maxiters && !stop
@warn "Non Convergence"
stop = 1
end

Y0 = Y1
G0 = G1
if ~stop
for j in 1:s
jac(J_temp, Y1.u[j], p, mesh[j+1])
JF[((j - 1) * problem_size + 1):(j * problem_size), ((j - 1) * problem_size + 1):(j * problem_size)] .= J_temp
end
temp = similar(Y1.u[j])
prob.f(temp, Y1.u[j], p, mesh[j + 1])
F1.u[j] .= temp
end
G .= vec(Y1 - B0) - W * vec(F1)
end

Y0_vec = vec(Y0)
nlprob = NonlinearProblem(bdf_first_nlprob_f!, Y0_vec)
nlsol = solve(nlprob, NewtonRaphson(); reltol=abstol, abstol=abstol, maxiters=maxiters)

recursive_unflatten!(Y1, nlsol.u)
for j in 1:s
cache.y[j + 1] = Y1.u[j]
prob.f(F1.u[j], Y1.u[j], p, mesh[j + 1])
cache.fy[j + 1] = F1.u[j]
end
end
Expand Down
85 changes: 31 additions & 54 deletions src/fode/newton_gregory.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -189,36 +189,22 @@ function NG_triangolo(
Phi_n = St + halpha * (zn[:, n1] + Phi)

yn0 = cache.y[n]
temp = zeros(length(yn0))
prob.f(temp, yn0, p, mesh[n1])
fn0 = temp
Jfn0 = Jf_vectorfield(mesh[n1], yn0, Jfdefun)
Gn0 = yn0 - halpha * omega[1] * fn0 - Phi_n
stop = false
it = 0
yn1 = similar(yn0)
fn1 = similar(yn0)
while ~stop
JGn0 = zeros(problem_size, problem_size) + I - halpha * omega[1] * Jfn0
yn1 = yn0 - vec(JGn0 \ Gn0)
prob.f(fn1, yn1, p, mesh[n1])
Gn1 = yn1 - halpha * omega[1] * fn1 - Phi_n
it = it + 1

stop = (norm(yn1 - yn0, Inf) < abstol) || (norm(Gn1, Inf) < abstol)
if it > maxiters && ~stop
@warn "Non Convergence"
stop = true
end

yn0 = yn1
Gn0 = Gn1
if ~stop
Jfn0 = Jf_vectorfield(mesh[n1], yn0, Jfdefun)
end

# Replace manual Newton iteration with NonlinearSolve.jl
# Solve: G(y) = y - halpha * omega[1] * f(y) - Phi_n = 0
function ng_nlprob_f!(G, y, p_nl)
temp = similar(y)
prob.f(temp, y, p, mesh[n1])
G .= y - halpha * omega[1] * temp - Phi_n
end
cache.y[n1] = yn1
cache.fy[n1] = fn1

nlprob = NonlinearProblem(ng_nlprob_f!, yn0)
nlsol = solve(nlprob, NewtonRaphson(); reltol=cache.reltol, abstol=cache.abstol, maxiters=cache.maxiters)

cache.y[n1] = nlsol.u
temp = zeros(length(nlsol.u))
prob.f(temp, nlsol.u, p, mesh[n1])
cache.fy[n1] = temp
end
end

Expand Down Expand Up @@ -251,38 +237,29 @@ function NG_first_approximations(cache::NewtonGregoryCache{iip, T}) where {iip,
JF[((j - 1) * problem_size + 1):(j * problem_size), ((j - 1) * problem_size + 1):(j * problem_size)] = Jf_vectorfield(
mesh[j + 1], cache.y[1], Jfdefun)
end
stop = false
it = 0
F1 = similar(F0)
Y1 = similar(Y0)
while ~stop
JG = Ims - W * JF
recursive_unflatten!(Y1, vec(Y0) - JG \ G0)


# Replace manual Newton iteration with NonlinearSolve.jl
# Solve: G(Y) = Y - B0 - W * F(Y) = 0
function ng_first_nlprob_f!(G, Y_vec, p_nl)
recursive_unflatten!(Y1, Y_vec)
for j in 1:s
prob.f(F1.u[j], Y1.u[j], p, mesh[j + 1])
end
G1 = vec(Y1 - B0) - W * vec(F1)

it = it + 1

stop = (norm(Y1 - Y0, Inf) < abstol) || (norm(G1, Inf) < abstol)
if it > maxiters && ~stop
@warn "Non Convergence"
stop = 1
end

Y0 = Y1
G0 = G1
if ~stop
for j in 1:s
JF[((j - 1) * problem_size + 1):(j * problem_size), ((j - 1) * problem_size + 1):(j * problem_size)] = Jf_vectorfield(
mesh[j + 1], Y1.u[j], Jfdefun)
end
temp = similar(Y1.u[j])
prob.f(temp, Y1.u[j], p, mesh[j + 1])
F1.u[j] .= temp
end
G .= vec(Y1 - B0) - W * vec(F1)
end

Y0_vec = vec(Y0)
nlprob = NonlinearProblem(ng_first_nlprob_f!, Y0_vec)
nlsol = solve(nlprob, NewtonRaphson(); reltol=abstol, abstol=abstol, maxiters=maxiters)

recursive_unflatten!(Y1, nlsol.u)
for j in 1:s
cache.y[j + 1] = Y1.u[j]
prob.f(F1.u[j], Y1.u[j], p, mesh[j + 1])
cache.fy[j + 1] = F1.u[j]
end
end
Expand Down
44 changes: 15 additions & 29 deletions src/fode/trapezoid.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -192,36 +192,22 @@ function TrapTriangolo(
Phi_n = St + halpha * (zn[:, n1] + Phi)

yn0 = cache.y[n]
temp = zeros(length(yn0))
prob.f(temp, yn0, p, mesh[n1])
fn0 = temp
Jfn0 = Jf_vectorfield(mesh[n1], yn0, Jfdefun)
Gn0 = yn0 - halpha * omega[1] * fn0 - Phi_n
stop = false
it = 0
yn1 = similar(yn0)
fn1 = similar(yn0)
while ~stop
JGn0 = zeros(problem_size, problem_size) + I - halpha * omega[1] * Jfn0
yn1 = yn0 - vec(JGn0 \ Gn0)
prob.f(fn1, yn1, p, mesh[n1])
Gn1 = yn1 - halpha * omega[1] * fn1 - Phi_n
it = it + 1

stop = (norm(yn1 - yn0, Inf) < abstol) || (norm(Gn1, Inf) < abstol)
if it > maxiters && ~stop
@warn "Non Convergence"
stop = true
end

yn0 = yn1
Gn0 = Gn1
if ~stop
Jfn0 = Jf_vectorfield(mesh[n1], yn0, Jfdefun)
end

# Replace manual Newton iteration with NonlinearSolve.jl
# Solve: G(y) = y - halpha * omega[1] * f(y) - Phi_n = 0
function trap_nlprob_f!(G, y, p_nl)
temp = similar(y)
prob.f(temp, y, p, mesh[n1])
G .= y - halpha * omega[1] * temp - Phi_n
end
cache.y[n1] = yn1
cache.fy[n1] = fn1

nlprob = NonlinearProblem(trap_nlprob_f!, yn0)
nlsol = solve(nlprob; reltol=cache.reltol, abstol=cache.abstol, maxiters=cache.maxiters)

cache.y[n1] = nlsol.u
temp = zeros(length(nlsol.u))
prob.f(temp, nlsol.u, p, mesh[n1])
cache.fy[n1] = temp
end
end

Expand Down