We’re rewarding the question askers & reputations are being recalculated! Read more.
20

I would strongly advice against using closed form solutions since they tend to be numerically very unstable. You need to take extreme care in the way and order of your evaluations of the discriminant and other parameters. The classical example is the one for the quadratic equation $ax^2+bx+c=0$. Calculating the roots as $$x_{1,2} = \frac{-b \pm \sqrt{b^2-...


18

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

If it's shock-capturing that you're interested in, I would suggest you use the finite volume method instead of the finite element method. When applied naively, FEM is actually notoriously bad at resolving shocks -- usually there are spurious oscillations or unwanted diffusion. Provided your original PDE is a conservation law, the FVM method will preserve ...


10

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


10

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

For a single rational equation in the complex domain, the basin of attraction is fractal, the compelement of a so-called Julia set. http://en.wikipedia.org/wiki/Julia_set . For theory with some nice online figures, see, e.g., http://mathlab.mathlab.sunysb.edu/~scott/Papers/Newton/Published.pdf http://hera.ugr.es/doi/15019160.pdf Even the ''globalized'' ...


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


7

If $q:=|f'(x^*)|<1$, where $x^*$ is the solution, the fixed point iteration you talk about is locally linearly convergent with convergence rate $q$. Thus if $q$ is small or zero, the method is competitive with Newton's method. Far away from the solution, convergence is difficult to predict in the absence of global information (such as a Lipschitz ...


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

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

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


6

I'm going to answer a more general question than the one you asked: do the eigenvalues of an initial value ODE determine the stability of the solution? Here I'm referring to mathematical stability, not numerical stability. Of course, a "yes" to this question is a necessary condition for a "yes" to your question. And unfortunately, the answer is "no". In ...


6

You need to solve a Riemann problem, perhaps approximately. For a linear system of equations, the solution to the Riemann problem is just upwinding applied to the characteristics. An "exact" Riemann solver for nonlinear problems resolves the full wave structure (consisting of shocks, rarefactions, and possibly linearly degenerate contact discontinuities). An ...


6

The Banach fixed point theorem is extremely general, needing only a metric space. If the map $x \to A(x)^{-1} b$ is globally contractive, the fixed point iteration converges. For your problem (if the formulation is typical of the problem class you describe), it will likely converge linearly, with a basin of attraction that is somewhat larger than basic ...


6

For solving nonlinear partial equations, say $-\Delta y + f(y) = 0$, there are two common approaches: Fixed point iteration: Pick $y^0$ and for $k=1,\dots$ solve $$ -\Delta y^{k+1} + f(y^k) = 0.$$ (That is, you replace $y$ by $y^{k+1}$ in all linear terms and by $y^k$ in all nonlinear terms.) Newton's method: Pick $y^0$ and for $k=1,\dots$ solve $$ -\Delta \...


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

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


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

The Feigenbaum fractal is a good example of how strange fixpoint iteration can be: http://en.wikipedia.org/wiki/Feigenbaum_fractal http://en.wikipedia.org/wiki/File:Logistic_Bifurcation_map_High_Resolution.png The second link plots the behavior of fixpoint iteration applied to the logistic map as one of the parameters varies. For certain values it ...


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