16

Quasi-Newton methods construct an approximate Hessian for an arbitrary smooth objective function $f(x)$ using values of $\nabla f$ evaluated at the current and previous points. At each iteration of the method the quasi-Newton approximate Hessian is updated using the gradient evaluated at the latest iterate, $x^{(k)}$. These approximate Hessians aren't ...


12

Short answer If you only want second order accuracy and no embedded error estimation, chances are that you'll be happy with Strang splitting: half-step of reaction, full step of diffusion, half step of reaction. Long answer Reaction-diffusion, even with linear reaction, is famous for demonstrating splitting error. Indeed, it can be much worse, including "...


11

Your example is a pretty good indication that the two derivatives (with respect to $x$ and with respect to $u$) do not commute :) (In fact, they're very different beasts -- one is a Fréchet derivative, the other a weak derivative.) Rather, you should consider $a(u) = (\partial_x u)^2$ as the composition $a = f\circ g$, of two functions $$ f: v(x)\mapsto v(...


10

According to Wolfram Alpha, $x^5+3(x-1)=0$ has no closed-form solution, so you can forget about a nice closed-form expression. :) I see nothing wrong with Newton's method; it should be quick and accurate, and with some analysis like the one you are sketching I think you can identify a safe starting point and prove global convergence. It might even be faster ...


9

I assume, that you have conducted a space discretization, so that you are about solving the (vector-valued) ODE $$ \dot u_h(t) = F_h(t,u_h(t)), \text{ on [0,T] }, u_h(0) = \alpha. $$ via a numerical scheme $\Phi$ that advances the approximation $u_h^n$ at the current time instance $t=t^n$ to the next value $u_h^{n+1}$ at $t=t^{n+1}:=t^n+\tau$. Then your ...


8

Yes: See Higham's book "Accuracy And Stability of Numerical Algorithms", second edition, chapter 25: Nonlinear Systems and Newton's Method. In particular, see the section on the "limiting residual" in terms of a condition number for the Jacobian. It may well be that since your system is ill-conditioned, that you quickly hit the limiting residual and your ...


7

This is somehow unexpected, but my recent experience with solving a system of nonlinear equations is that treating them as the right hand side of a system of ordinary equations and then evolve the system with an ODE solver can be considerably faster than with the usual Newton-Raphson iteration. It sounds like you're doing some sort of pseudotransient ...


7

The issues you're running into now are not a failing of Newton-Raphson, but a question of coupling. You're doing iterated sequential coupling -- solving each equation sequentially and then iterating until (hopeful) convergence. No solver choice in place of NR is going to fix this lack of convergence, as long as you are doing iterated sequential coupling. ...


7

It's a bit easier to see if you write your equation in the a semi-discretised system of the form $u^{\prime}(t) = F(u(t))$ and with the application of the $\theta$-method and approximating $u^{\prime}(t) \approx (w^{n+1} - w^{n})/\tau$ this gives, $$w^{n+1} - w^{n} - (1-\theta) \tau F(w^n) - \theta\tau F(w^{n+1}) = 0$$ with unknown vector $w^{n+1}$ and ...


7

This is sometimes known as Buzano's inequality (http://www.jstor.org/stable/2159168). In general, if $\|x\|=1$, $P=xx^{\top}$ is a projection operator, so a simple application of Cauchy-Schwarz leads to $$ 2|b^{\top}xx^{\top}a|-|b^{\top}a| \leq |b^{\top}(2P-I)a| \leq \|a\|\,\|b\|, $$ which gives Buzano's inequality $$ |b^{\top}xx^{\top}a| \leq \frac{\|a\|\,\|...


7

Some of the confusion comes about because mathematicians often consider artificial functions that happen to have size $O(1)$ whereas practitioners use functions and variables that just happen to be whatever they are, often vastly different from one. One of the considerations that is important in practice is to make sure that the choice of physical units does ...


6

There are multiple ways of implementing Dirichlet boundary conditions when using the finite element method. For each unknown $i$ of the system belonging to the Dirichlet boundary, you can zero out row $i$ of the Jacobian and entry $i$ of the right-hand side. At that rate, you can reformulate the problem for just the unknowns in the interior, which reduces ...


6

A few specific examples: Quasi-Newton methods avoid computing second derivatives by using an approximation to the Hessian which is updated at each iteration by a low rank (rank one or rank two) update. This makes factoring the Hessian (or equivalently keeping the inverse in product form) easy to do computationally. There are also limited memory Quasi-...


6

If you need Jacobian matrix information for a numerical method, you should calculate the Jacobian matrix of the discretized form of the equations, since that will be consistent with the discretized equations you are solving.


6

Off the top of my head, there are a number of things that may be going wrong. You might want to verify the following Are the function values or derivatives you are computing numerically stable? Is the rounding error in their evaluation smooth? You can verify this visually by plotting an extremely small interval around a known root such that the function ...


6

The approximation error is $\mathcal{O}(\varepsilon)$, where $\varepsilon$ is the parameter in the finite difference calculation. It will be good to know how this approximation affects the number of CG iterations compared to having the complete Hessian? Are there theoretical/empirical understandings about this? Yes. In short, there are problems where CG ...


6

A low rank update to an $n$ by $n$ matrix $A$ is an update of the form $A=A+UU^{T}$ where $U$ is matrix with $n$ rows but very few columns (typically just one or two.) If the matrix $UU^{T}$ has $k$ columns, then its rank is at most $k$, so this is a low-rank update to the $A$ matrix. Low rank updates are important in computational linear algebra and ...


6

To calculate the reaction forces at a node, Abaqus (or any structural FE code) simply sums the internal forces for all elements attached to that node. The reaction forces are the negative of that sum. For an Abaqus user element, the internal forces for the element are returned from subroutine UEL in the RHS array. The returned stiffness matrix (Jacobian), ...


6

The reason is that GMRES can only be used for solving linear equations, i.e. equations of the form $Ax=b$, where $A$ is some matrix and $x,b$ are vectors. What GMRES does, essentially, is it approximates multiplication by the matrix $A^{-1}$ using a matrix polynomial of $A$. In this case (I assume) $f(y^{n+1},t)$ is not necessarily linear in the vector $y^{...


5

Your observed quadratic convergence indicates that the Jacobian is likely correct. Have you looked at the solutions for your under-resolved configurations? Galerkin optimality uses the operator norm, which contains only the symmetric part, thus the solution of the discrete system could be quite different from the projection of the exact solution. This ...


5

I'm not sure if this is the answer, but consider the case as $\alpha \rightarrow \infty$. In this case you're trying to solve $\iint_{\Omega} \cosh(P) dxdy = 0$ And since cosh looks like this: You can see that integrating it over any domain will never give you zero. If the diffusion of $P$ ($\nabla^2 P$) is not large enough to counteract this limit then ...


5

I suspect the main problem is the magnitude of the values. If you divide through by $\cosh$ to make all the numbers smaller, then ApproxFun doesn't seem to have a problem finding all the roots. function lhs(ε, α) A, B = sqrt((α-1)/(2ε)), sqrt((α+1)/(2ε)) return sqrt(α*α-1)*(cos(A) - 1/cosh(B)) - sin(A)*tanh(B) end u = Fun(x -> lhs(1e-3, x), [1.0,...


5

I don't know much about tracking implicit surfaces, so I'm just going to start with the optimization problem and go from there. The optimization problem is, at the core, nonlinear least squares, and can therefore be solved efficiently by the Gauss-Newton method, wherein at each step one linearizes the problem at the current guess, solves the linear least ...


5

The Netwon(-Raphson) method helps you finding the root(s) of a function $f(x)$, that is, those values of $x_0$ for which $f(x_0) = 0$. In your case, the unknown value is $c$, so we start by writing down the problem in terms of a function $f(c)$ that is zero when $c$ takes the "correct" value (which we call $c_0$). Since this looks like homework, I don't ...


5

Not sure where you get your equation for $\epsilon$, but ultimately your approximation for the Jacobian matvec operation is a finite difference approximation to the directed derivative of $F(\cdot)$. This means that you will want as accurate of an approximation to this directed derivative as possible while being sure to avoid numerical round-off issues. For ...


5

If you have a single vector equation $\vec{F}(\vec{x})=0$ then you solve it by representing that state vector $\vec{x}$ as a set of amplitudes $[x_0,x_1,...,x_{n-1}]$ after discretization by your favorite method (FD, FV, FEM, spectral); and we know how to solve it. If you also have a second equation $\vec{G}(\vec{y})=0$ then the full state vector is $[x_0,...


4

Using a line search or trust region method along with Newton's method typically does the trick of finding the zero reliably. Check any advanced numerical methods or optimization book.


4

Can you bound your zero? If you can come up with a good heuristic for the approximate location of the zero and then compute bounds for it over an interval that is not too large, the performance of many methods is improved. Here's some Matlab code I use to perform basic root bracketing for monotonic functions. In terms of interval methods, @Pedro's ...


4

They're equivalent. Both imply some variation of the root-finding method by linearization.


4

This remark is a bit too long for a comment, but regarding DanielShapero's remark with respect to changing the Jacobian matrix, for every degree of freedom $i$ corresponding to a Dirichlet boundary condition (BC), a common approach is to set $i$th row of the Jacobian matrix equal to the $i$th row of the identity matrix, meaning your (quasi-)Newton iterates ...


Only top voted, non community-wiki answers of a minimum length are eligible