I want to find the solution to a nonlinear equation of the form u_{xx} + f(u,\lambda) + p= 0 where the unknowns are u(x) and the constant p. There is an addition volume constraint \int u(x)dx=1. For the leading order problem, f(u) is an explicit function of u and some parameters \lambda. After discretization using finite difference, I manage to solve it using packages like NLsolve.jl or even the full family of solutions (with different \lambda) using continuation in BifurcationKit.jl.
Suppose I want to tackle the higher order problem, now f_i=f(u_i) (the parameter \lambda is omitted for a sec just for simplicity) is no longer an explicit function of u (f_i, u_i are the discretized f(u) and u). To obtain f_i, I need to solve a matrix equation Af=b where the entries of the (sparse) matrix A and vector b depend on u. So the whole function I’m trying to solve looks something like this
function f(x, par)
u = @view u[1:end-1]
p = @view u[end]
uxx = finiteDiff(u, 2, par.order, par.dx) # my custom function for F.D.
A, b = assemble_matrices(u, par) # construct A and b using u(x) and given parameters
f = A \ b
residual = similar(x)
residual[1:end-1] = uxx .+ f .+ p
residual[end] = sum(u)/length(u) - 1
return residual
end
I’m aware of things like using non-allocating arrays. The above is just a high-level overview of the function f. Now the problem is that when I use the packages I mentioned above, they usually rely on automatic differentiation, and then very often I get StackOverFlow error inside ForwardDiff when it gets to the f=A\b line. Using finite difference for the Jacobian seems to work, but due to the matrix equation solve it is much slower my leading order problem. I just wonder if there are some general techniques for solving nonlinear equations like this? Do I just need to try my best to simplify the construction of matrix A, use as many non-allocating arrays as possible, or I just have to reformulate the problem to get rid of f=A\b?