People who code: we want your input. Take the Survey
19

There are two issues that you are likely to be encountering. Ill-conditioning First, the problem is ill-conditioned, but if you only provide a residual, Newton-Krylov is throwing away half your significant digits by finite differencing the residual to get the action of the Jacobian: $$ J[x] y \approx \frac{F(x+\epsilon y) - F(x)}{\epsilon}$$ If you ...


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

This is the stochastic root-finding problem, as in The stochastic root-finding problem: Overview, solutions, and open questions.


11

Is the shooting method the only general numerical method for solving BVP of nonlinear ODE(s)? No. Most other methods consist of three parts: Discretization. This may be done with finite differences, finite volumes, finite elements (Galerkin or collocation), spectral methods, and so forth. This reduces the problem from an infinite-dimensional one to a ...


9

If you can stand to use complex arithmetic, simultaneous iteration methods might be preferable for computing all the roots of your polynomial. The simplest simultaneous iteration method, the (Weierstrass-)Durand-Kerner method, is effectively equivalent to applying Newton-Raphson to the Vieta relations relating the coefficients and roots of a polynomial, ...


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

You can solve this numerically in Python without symbolic computation. from __future__ import print_function, division import numpy as np from numpy import exp from scipy.integrate import quad from scipy.optimize import root def f1(a1, a2, x): return exp(a1 * x + a2 * x * x * x) / (1 + x * x) def f2(a1, a2, x): return exp(a1 * x + a2 * x * x * x) *...


8

PDEs are a form of dynamical system where there is another continuous variable. Usually this is space, so you're looking at how things over time and space instead of just over time. Here's an example of generalizing an ODE to a PDE. Take your ODE model of chemical reactions which models the concentration of certain chemicals over time. Now generalize the ...


8

Julia has a whole ecosystem for generating sparsity patterns and doing sparse automatic differentiation in a way that mixes with scientific computing and machine learning (or scientific machine learning). Tools like SparseDiffTools.jl, ModelingToolkit.jl, and SparsityDetection.jl will do things like: Automatically find sparsity patterns from code Generate ...


7

To answer your questions in order: Any implicit method for solving an ordinary differential equation involves solving a system of nonlinear equations. You can do this through variants of Newton's method, successive substitution, full approximation scheme, or any other approach that solves systems of nonlinear equations. (The caveat is, of course, depending ...


7

What you describe as your time discretization is called the Crank-Nicolson scheme. For nonlinear differential equations it leads, as you have observed, to a nonlinear algebraic system that needs to be solved at each time step. The typical approach is to solve it with a Newton method -- in your case, that requires to solve a nonlinear system in 5 variables. ...


7

As pointed out by David Ketcheson, one method is to use the companion matrix and find its eigenvalues (this is what Matlab does for the roots function). However, if you want to code everything by yourself, you can try to use the Sturm sequences, which you use as a first step to find an interval with only one zero. Then, you can apply one standard methods ...


7

For your example equation, taking the average approach, the local consistency error $$ \frac{1}{h^2}[u(x-h) - 2 u(x) + u(x+h)]-f(\frac{1}{2}[u(x-h)+u(x+h)]) = \frac{1}{2}f_uu_{xx}h^2 + hot. $$ will be of order $2$ (instead of order $3$). ($hot.$ means higher order terms) Therefore, if your overall approximation is of order $1$, e.g. if you use upwind ...


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

$s_k$ is the "approximate Newton" search direction. So in essence, when they say Choose $s_k$ such that $\|F(x_k)+F'(x_k)s_k\| \le \eta_k \|F(x_k)\|$ they are saying: Solve the Newton system $F'(x_k)s_k = - F(x_k)$ inexactly for $s_k$ until the norm of the residual $F(x_k)+F'(x_k)s_k$ is smaller than the norm of the right hand side $-F(x_k)$ by a ...


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

You should rather think of what you will need in the following step, which is probably the numerical time integration of the semi-discrete equations. If you are going to use a (semi) explicit time stepping scheme, all you need is a function that for a given $\phi_0$ assembles the vector $\langle (\nabla \phi_0)^2, v \rangle$, where $v$ are your test ...


6

$K u$ equals the internal forces only in the linear case. The tangent stiffness matrix, $K$, in a nonlinear problem is normally used in a Newton-Raphson algorithm to calculate updates to the displacement vector as follows: $$K \Delta u = f - f_{internal}$$ $$ u_{i+1} = u_i + \Delta u$$ The vector of internal forces, $f_{internal}$ must be calculated from ...


6

If you change variables to optimize for the residual of the linear part, then the Hessian will be a low-rank update to the identity. Then L-BFGS would work very well. Specifically, your problem takes the form $$ \min_x \frac{1}{2}\|Ax-b\|^2 + \frac{\mu}{2}\|g(x)\|^2 $$ where $Ax=b$ is the linear PDE and $g$ is the nonlinear part, and $\mu$ is a tradeoff ...


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^{...


6

Two systematic ways of smoothing a function $h$ would be: 1. Join the piecewise smooth parts of your function using Hermite interpolation so that the derivatives are matched to your satisfaction. 2. Convolve your function $h(x)$ with a heat kernel of the form $f(x) = \frac{\exp\left\{-\frac{x^2}{2 \sigma^2}\right\}}{\sqrt{2 \pi \sigma^2}}$ so that instead ...


5

That's actually very easy to do in Dolfin: from dolfin import * # define mesh, function space (piecewise linear) mesh = UnitSquareMesh(64,64) V = FunctionSpace(mesh,'CG',1) # inhomogeneous boundary conditions (otherwise the solution is trivial) bc = DirichletBC(V, Constant(1.0), lambda x,on_boundary: on_boundary) # define bi(non)linear form # note that ...


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

As it looks to me your formulation of the discrete system is incomplete. You need to state your equation at each (grid?) point associated with $n$ put back the summation over $l$ (grid) points in your last equation assign a time index $i+1$ or $i$ to your $c_l$s In this way, you get a nonlinear (if you take $c_l^{n+1}$) or linear (if you take $c_l^{n}$) ...


5

A common strategy is to employ a damping strategy, i.e., to compute $\vec{w}^{\ast} = F \left( \vec{w}^{(n)} \right)$ and then set $\vec{w}^{(n+1)} = \alpha \vec{w}^{(\ast)} + (1-\alpha)\vec{w}^{(n)}$ where $\alpha\in[0,1]$. You typically choose $\alpha$ in such a way that it minimizes $\|[\alpha \vec{w}^{(\ast)} + (1-\alpha)\vec{w}^{(n)}] - F(\alpha \vec{w}^...


5

Certified homotopy continuation methods are used both for finding roots and for proving that they indeed exist (inside a certain interval). A quick web search turned out this paper: Reliable homotopy continuation by Joris van der Hoeven.


5

No it is not. There is also multiple shooting collocation finite differences fixed point iterations and probably some more.


5

Since $\phi$ is a scalar between $0$ and $1$, the easiest method for finding a root is bisection. If you cannot calculate the expectation of the nonlinear function $$f_\phi(\theta) = \left(\phi r_z +(1-\phi)(r_k+\theta)\right)^{1-\gamma}(r_k+\theta-r_z)$$ in terms of $\phi$ analytically, you can use quadrature to approximate. For the first variant, this ...


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