Skip to content

Commit

Permalink
rootfind with p (#125)
Browse files Browse the repository at this point in the history
* rootfind with p

* add  suggestions

* fix keyword error
  • Loading branch information
lindnemi authored Oct 16, 2023
1 parent cda7cc0 commit da56a6d
Showing 1 changed file with 15 additions and 14 deletions.
29 changes: 15 additions & 14 deletions src/Utilities.jl
Original file line number Diff line number Diff line change
Expand Up @@ -169,14 +169,15 @@ function find_fixpoint(nd, p, initial_guess)
return nl_res.zero
else
println("Failed to find fixpoint. Algorithm did not converge.")
println("Try running nlsolve with other options.")
end
end

export RootRhs

"""
struct RootRhs
RootRhs(of)
RootRhs(nd)
A utility function that provides a root function, to find valid initial
conditions. The operator ``M^\\dagger M - 1`` projects onto the kernel of ``M``.
Expand All @@ -187,36 +188,36 @@ by returning the projection of ``f(x)`` onto the kernel as residue.
The use of the pseudoinverse here means this does not play nice with sparse
matrices. Beware when using for large systems.
"""
struct RootRhs
rhs
mpm
struct RootRhs{T,S}
rhs::T
mpm::S
end
function (rr::RootRhs)(x)
function (rr::RootRhs)(x, p=nothing, t=nothing)
f_x = similar(x)
rr.rhs(f_x, x, nothing, 0.0)
rr.mpm * f_x .- f_x
rr.rhs(f_x, x, p, t)
return rr.mpm * f_x .- f_x
end

function RootRhs(of)
mm = of.mass_matrix
function RootRhs(nd)
mm = nd.mass_matrix
@assert mm !== nothing
mpm = pinv(mm) * mm
RootRhs(of.f, mpm)
return RootRhs(nd.f, mpm)
end


export find_valid_ic


"""
find_valid_ic(of, ic_guess)
find_valid_ic(nd, ic_guess)
Try to find valid initial conditions for problems involving mass matrices.
Uses RootRhs as residual function.
"""
function find_valid_ic(of, ic_guess)
rr = RootRhs(of)
nl_res = nlsolve(rr, ic_guess)
function find_valid_ic(nd, initial_guess; p=nothing)
rr = RootRhs(nd)
nl_res = nlsolve(x -> rr(x, p), initial_guess)
if converged(nl_res) == true
return nl_res.zero
else
Expand Down

2 comments on commit da56a6d

@lindnemi
Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@JuliaRegistrator
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Error while trying to register: "Tag with name v0.8.1 already exists and points to a different commit"

Please sign in to comment.