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


8

You haven't told us exactly what optimization routine you're using, so it's difficult to provide a very specific answer to your question. However, if you don't supply your own Jacobian function then the optimization routine that you're using is presumably approximating the derivatives using a finite difference approximation scheme. Computing finite ...


7

Like with any other thing you want to test, you need to come up with a list of situations where you know that something is true, and then check this. In many cases, knowing that something is true does not require being able to actually compute the numerically correct answer, which simplifies things. For example: As @BlaisB suggestes, if you rotate the ...


7

Julia's DifferentialEquations.jl has a lot of tooling for automatically deriving (sparse) matrices. For more information, see the JuliaCon 2020 video on Auto-Optimization and Parallelism in DifferentialEquations.jl. Combined with the orders of magnitude acceleration commonly seen over the MATLAB solvers, this might be a good option for you and is a quick ...


6

One word: Modularity. There are a lot of repeated expressions in your Jacobian that could be written as their own function. There's no reason for you to write the same operation more than once, and that will make debugging easier; if you only write it once there's only one place for an error (in theory). Modular code will also make testing easier; ...


6

Given that $u \frac{\tau^2}{2} \ll 1$, one way of tackling the numerical oscillations, even before the actual term emerges, is a Taylor approximation in $u$ of the sine term (Thanks to Kirill for the full series): $$ \sin(\underbrace{x_5 + x_6 \tau}_{u_0} + \underbrace{u \frac{\tau^2}{2}}_{\Delta u}) = \sin(x_5 + x_6 \tau) + \sum_{n=1}^N \frac{(-1)^{\lfloor ...


6

Finite difference approximations of the Jacobian are really only good if the step lengths are chosen appropriately for each coordinate. But a black-box solver like CVODE has no way of knowing what these step lengths should be, and so has to use heuristics to choose them. This may or may not work. You are almost always better off if you provide an ...


5

You've started with a singular linear system of equations $Ax=b$. As a practical matter, it's unlikely that $b$ lies exactly in the range of $A$, so at best you can find a least squares solution that minimizes $\min \| Ax - b \|_{2}$ Because the system is singular, the null space of $A$ is non-empty, and there will be an infinite number of solutions to ...


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

Well, the first obvious test is to check if the determinant of the Jacobian matrix is positive. If it not, it means that your element has inverted and your answer is going to be invalid. Otherwise, I am not sure of any other test that would be good... For a unit test you could rotate a triangle 360 degrees by increment of 1 degree and measure that the ...


5

There are a couple of options here: Don't compute elements of the Jacobian corresponding to the fixed parameters and don't include these parameters in your vector of parameters to be fitted. In your code this would require a different version of the function and Jacobian for each collection of parameters that might be fixed in a particular call to the LM ...


5

There are several strategies to consider: Find the derivatives in symbolic form using a CAS, then export code for computing the derivatives. Use an automatic differentiation (AD) tool to produce code that computes the derivatives from code to compute the functions. Use finite difference approximations to approximate the Jacobian. Automatic ...


5

My first question: Will split ODE solvers likely work for my problem? From your description, this sounds like a textbook use-case for a split ODE solver. Neither an implicit method nor an explicit method is practical across the entire ODE. IMEX methods sound like a reasonable choice as long as the atmospheric temperature dynamics that you want to treat ...


4

You could consider $d=1$ and $$ f(x) = e^x, \qquad y=0. $$ I don't think a function's derivative being positive everywhere implies that the function's range is all of $\mathbb{R}$. It does imply, though, that solutions are unique. If $(x_1,y)$, $(x_2,y)$, $x_1\neq x_2$ are two solutions, then $$ 0 = (x_2-x_1)\cdot(f(x_1) - f(x_2)) = \int_0^1 (x_2-x_1)\cdot ...


4

I am not much into the thing but when we use l-m for training neural nets there is a trick we use: The original equation is $$\Delta w = - (J^T J)^{-1} J^T e$$ by subtracting the term $uI$ in order to make the matrix invertable where u is a constant and I is the identity matrix so the equation will become $$\Delta w = - (J^T J+uI)^{-1} J^T e$$ so when $u$ ...


4

It looks like you have a advection diffusion PDE discretized with finite differences. This gives an ODE of the form $$ y' = f(y) = A y + D y, $$ where $A$ is the discretized advection operator and $D$ is the discretized diffusion operator. In your case, it seems you have zero Dirchlet and Neumann boundary conditions which can be encoded in the $A$ and $D$ ...


3

You need to understand how to actually compute derivatives when you want to take the derivative with respect to a function. I've recorded a lengthy example in lecture 31.55 here: http://www.math.colostate.edu/~bangerth/videos.html In short, you need to consider the derivative as the limit of a variation.


3

From what you describe, you have an ill-posed problem: The solution is not unique (not even locally). A standard way of dealing with this is the following: Instead of trying to solve $F(x)=y$, where $x$ is your vector of geometrical parameters, $y$ the vector of optical parameters, and $F$ the mapping that computes the latter via the former, you minimize $$ ...


3

So, if I framed this question properly, is it possible at some point during the RKF calculation to construct a Jacobian to pass out to Scipy along with the objective function evaluation? The Runge-Kutta-Fehlberg 4(5) method you describe is an explicit numerical method for solving an ODE. As such, it does not require the solution of a nonlinear algebraic ...


3

In your case (line 12 in your code) $$A = D-R.$$ Then $$ \begin{aligned} A\mathbf{x} {}&{}= D\mathbf{x}-R\mathbf{x}\\ {}&{}= \mathbf{b}\,, \end{aligned}$$ which rearranges to $$\mathbf{x} = D^{-1}\left( \mathbf{b} + R\mathbf{x} \right)\,,$$ which leads to the iteration which is the the first line in your while loop. This only differ's ...


3

Here is an example of where we have used automatic differentiation using Sacado in one code: http://www.dealii.org/developer/doxygen/deal.II/step_33.html


3

If your advection problem had Dirichlet of Neumann boundary conditions, the linear system would be tridiagonal and you could apply the Thomas algorithm. With periodic boundary conditions, however, we lose this. If c(x) is a constant independent of x, the matrix would be circulant and linear systems could be solved efficiently using FFTs. An even better ...


3

This problem is too small to actually be sparse. Sparse handling has a big overhead because the indexing is not "direct", i.e. you don't necessarily know where the next value will be without branch checking. So you need it to be "sparse enough" that the O(n^3) dense LU-factorization cost shrinking to the purely non-zero terms overcomes ...


3

There are some unknowns in what you are doing but for simplicity, suppose we want to find $u(t)$ as discrete times $t_1, t_2, \cdots, t_n$. Let $\textbf{F} = [F(t_1), F(t_2), \cdots, F(t_n)]^T$ and $\textbf{u} = [u(t_1), u(t_2), \cdots, u(t_n)]^T$ be column vectors representing $F$ and $u$ evaluated at the desired times. From your problem statement, you wish ...


2

For the past few months, I've been working on a similar problem. I only solve the steady-state potential equation for the moment, assuming homogeneous, steady-state ion concentrations. Scaling your equations is going to be a rudimentary preconditioner. I found that I got good results using algebraic multigrid (AMG) for preconditioning instead, since my ...


2

The answer depends very much on what optimization function in MATLAB you're talking about. The base MATLAB system includes a function fminsearch that uses the simplex method (not the simplex method for LP but rather the Nelder-Mead-Lagarias algorithm of the same name.) This method does not make use of derivatives. The functions of the optimization ...


2

I'm thinking pseudo-spectral method might solve the problem? A priori, I don't see why it would. As a continuous expression, the Butler-Volmer equation varies tremendously with small changes in potential. Any discretization inherits this behavior, so I think you're stuck with ill-conditioning no matter what. Maybe you could nondimensionalize the expressions,...


2

If you look at the full output of your script, sol has fields ipvt and qtf. Both of these fields come from MINPACK (http://www.netlib.org/minpack/lmdif.f), which is the backend here (https://github.com/scipy/scipy/blob/2526df72e5d4ca8bad6e2f4b3cbdfbc33e805865/scipy/optimize/_root.py#L246; https://docs.scipy.org/doc/scipy-0.19.0/reference/generated/scipy....


2

If you don't care too much about which $b$ you're working with (say just for linear algebra), then you can use the Method of Manufactured Solutions (MMS) in Linear Algebra much like we differential equations geeks would for our problems to start with a known exact solution, $x^*$, derive a $b$, and then apply your solver to see how decent an $x$ you can ...


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