# Solving nonlinear equation with an embedded matrix equation

**URL:** https://discourse.julialang.org/t/solving-nonlinear-equation-with-an-embedded-matrix-equation/88240
**Category:** Numerics
**Tags:** nonlinear, nonlinearsolve
**Created:** [October 4, 2022, 4:24pm UTC](https://discourse.julialang.org/t/solving-nonlinear-equation-with-an-embedded-matrix-equation/88240 "2022-10-04T16:24:35Z")
**Posts on this page:** 2
**Page:** 1

<div class="post-metadata">

### Author: ![yhchang96](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/yhchang96/32/36499_2.png) [@yhchang96](https://discourse.julialang.org/u/yhchang96)
#### Post date: [October 4, 2022, 4:24pm UTC](https://discourse.julialang.org/t/solving-nonlinear-equation-with-an-embedded-matrix-equation/88240/1 "2022-10-04T16:24:35Z")

</div>

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

```julia
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`?

---

<div class="post-metadata">

### Author: ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [October 4, 2022, 9:54pm UTC](https://discourse.julialang.org/t/solving-nonlinear-equation-with-an-embedded-matrix-equation/88240/2 "2022-10-04T21:54:52Z")

</div>

Do you have an MWE that you can share? Nothing about `A\b` is bad for automatic differentiation, so what you have described is somewhat of a red herring. Whatever is giving the stackoverflow is something not in the description here. Some questions that would be helpful:

- What are the types of `A` and `b`?
- Can you isolate this `A\b` stackoverflow to be outside the nonlinear solver?
