Tag Info

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

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

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

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

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

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

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

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

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

2

odenumjac calls your function in a vectorized manner it seems, and your function is not vectorized. You can easily change that by changing the second index of f in your function to : instead of 1, for instance: f(10,:) = 2*(x(end-1,:) - x(end,:)); I thought the setting joptions.vectvars=1 would not allow the vectorised call (see one of your other questions). ...

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