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


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


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

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

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

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


4

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


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


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

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

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


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

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

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


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


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

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

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

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


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

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


1

So it looks to me like there is a misunderstanding. your residual, or right hand side, presumably depends on multiple cells, so when building your jacobian, you should ask yourself how each entry in your residual depends on each cell. Now this should give you an nxn sparse matrix. If you in fact linearize this by hand, you should be able to store it in a ...


1

I finally found the solution to that problem: In order to be able to reuse as much code as possible, I had one single function for calculating the right hand side of my problem, i.e. $-F(u)$, and for calculating the residual value. This function returned $-F(u)$ when calling, but I forgot to remove the minus-sign when using the same function for the ...


1

I known that this is not a definitive answer, but it can give an idea how to move with CUDA. At this stage is difficult to give advice because for this is necessary a detailed know about the actually code As already write in comments, in general, is better not invert the matrix, but solve the linear system (for detail about why see this question). There is ...


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