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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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


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

You're using a backward Euler finite difference time stepping method. This is stable, but only first order accurate, so I suspect explains the discrepancy in the reduction factor. More explicitly, assuming a perfect spatial representation and substituting your initial condition into your difference equation gives $$(1+2\pi^2\Delta t)u^{n+1} = u^n.$$ For a ...


2

It seems that the Jacobian ratio is defined as the ratio between the maximum and minimum Jacobian determinant in an element [1, 2]. And, that a value between 0.33333 and 1 is good-enough [3]. Nevertheless, for linear elements, the Jacobian is constant and thus the same over each element. As mentioned by @GustavoCosta, 3 descriptors commonly used for element ...


2

Taking $x_0 = 0$, we have that $x_1 \in <M^{-1}b>$. For the next iteration, we get $x_2 \in <M^{-1}b, M^{-1}AM^{-1}b>$. For the next iteration, we get $x_3 \in <M^{-1}b, (M^{-1}A)^2M^{-1}b>$. Continuing this argument, you will eventually find that $x_k \in \mathcal{K}_k(M^{-1}A,M^{-1}b)$. This space can also be described as $Mx_k \in \...


2

You can arrive at the Jacobian analytically it just takes a few steps So assuming we have our typical FE field values: $$ u_i = \sum_j \phi^j u_i^j $$ Where $i$ represents coordinate direction, $\phi$ our shape function, and $j$ index to DOF. Start with our stabilization parameter, assuming you are using Euclidian norm: $$ \tau = \sqrt{\alpha \lVert u \...


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