7

Some of the confusion comes about because mathematicians often consider artificial functions that happen to have size $O(1)$ whereas practitioners use functions and variables that just happen to be whatever they are, often vastly different from one. One of the considerations that is important in practice is to make sure that the choice of physical units does ...


6

Looks like everything is working relatively OK? Your matrix is of order 1e10, so residuals of 1e-4 are actually close to machine precision. The convergence criterion is indeed violated, but not by much; not sure what's going on there, are you sure Arpack really guarantees it or is it a best effort kind of thing? It surprises me that you get different results ...


4

$h$ is a measure of the mesh size. In the example, they are using rectangular elements. For which a commonly used measure of the mesh size is the length of the largest diagonal. Looking at the table above the sentence you quoted: $(M_5-M_6)/(M_4-M_5)=(0.14052586 −0.14056422)/(0.14037251−0.14052586)\approx 1/4 = (h_6/h_5)^2$ where $M_i$ is the mean at $i$-th ...


4

This is not a definite answer, but I highly recommend the work of Deuflhard on Newton methods (e.g. this report, this book, and another report I cannot find just yet). The basic idea is the following: a criterion based on the norm of the Newton step is scaling independent, as the sequence of Newton iterates is actually independent of the scalings used for ...


4

These humps occur in the context of optimization when the solution is approached but we then go past it. This frequently occurs in gradient descent with momentum, where they appear regularly on plots because the gradient and momentum both decay in magnitude together. Usually this suggests a little too much momentum, since we're retaining speed as we reach ...


3

The choice of $k$ is restricted also by the discretization of the source term. To see it, rewrite your scheme to \begin{equation} u_m^{n+1} = \left(1 - \frac{k(1-x_m)}{h} - k(1-x_m)\right) u_m^n + \frac{k (1-x_m)}{h} u_{m+1}^n \,. \end{equation} You need $$ 1 - \frac{k(1-x_m)}{h} - k(1-x_m) \ge 0 $$ for all $x_m$. Taking $x_m=0$ (the worst case scenario) you ...


2

OK, so with new version in edit: it's not your fault, it's https://github.com/JuliaLinearAlgebra/Arpack.jl/issues/87. You can either call Arpack manually yourself, or even better use a pure-julia solver (KrylovKit.jl and IterativeSolvers.jl are good choices)


2

The series that converges to $\ln(2)$ appears to be suitable for Cohen-Villegas-Zagier acceleration [PDF]. This is an acceleration technique for alternating series, but continued fractions with positive partial numerators and denominators are equivalent to alternating series. In particular, if $S_m$ is the $m$-th continued fraction approximant, then let $$...


2

This is one of those scenarios where assessing the convergence order of your scheme is difficult, because, as you explained, your time-step is "linked" to your spatial discretization. What you could do is manufacture an analytical problem without any spatial error. For example, in finite elements, if you are using P2 polynomials, a second order ...


2

Based on the references you provided, I assume your initial conditions are random numbers (see second reference, Sect. IV)? If so, it's probably no surprise that you do not see higher-order convergence in time, at least not from the beginning. Your initial conditions are not smooth and it may take some time to get to a smooth solution/state. There, however, ...


1

I assume you are solving a two dimensional problem: Your refinement in your first example is not $\frac{\Delta x}{2}$. Generally this does not make problems if you refine uniformly. If you look a the number of elements, the error and the EOC something does not match. I am no expert for deal.ii, but i think the EOC routine is not able to handle different ...


1

@cfdlab's answer gives a general way to construct solutions for discontinuous coefficients, but here is a slightly more theoretical perspective to it as well. All of this can be understood in 1d, so imagine all of the functions below to be functions of just one argument $x$ for a moment. First, you can't just choose any $u$ and $\alpha$ to obtain a function $...


1

If you want to use the coefficient of the form $$ \alpha = \begin{cases} \alpha_1 & r < 0.5 \\ \alpha_2 & r > 0.5 \end{cases} $$ it is useful to work in polar coordinates. You build two pieces of the solution $$ u_1(r,\theta), \qquad r < 0.5 $$ and $$ u_2(r,\theta), \qquad r > 0.5 $$ Then at the interface you ensure solution and flux ...


1

The way I computed the solution in the linked answer is the classical one: I just took a reference solution (assuming the code was correct, i.e. the numerical solution I found is the right one) with a small enough step size $h$, say $h=10^{-9}$. Then I computed the solution with smaller $h$s and for each one of those $h$ I computed the $|| e_h ||$ in a ...


1

To add to Wolfgang's answer, some problems have specific structure that you can use to design physically-based convergence criteria. Suppose you're solving the generalized Poisson equation or some other elliptic partial differential equation that can be derived through minimization of some convex action functional $J(u)$. Very often, the action can be split ...


1

The finite difference scheme also gives exact solution at the nodes to the problem $-u''=1$ because $$ \frac{u_{j-1} - 2 u_j + u_{j+1}}{h^2} $$ is exact for a quadratic function. For the more general case $$ -u''(x) = f(x) $$ if you compute the integral $(f,\phi_i)$ exactly in $P_1$ FEM, it gives exact solution at the nodes. The FD scheme will not be exact ...


1

eigsh( ... sigma=0 ) looks like a bug in python-scipy arpack too. Here's a comparison of arpack with numpy eigh (Lapack syevd, gold standard), run on your matrix / 1e6: sigma=0 -- evals numpy eigh: [-8.23e-13 370 515 2.82e3 3.64e3 5.21e3 5.9e3 7.22e3] evals arpack: [-146 -5.68e-14 96.1 844 3.43e3 5.06e3 5.18e3] ❓ eigcheck AV - VΛ numpy: max 2.5e-12 ...


1

First some comments on why such matrices are hard for solvers. Notice that if $U$ is upper diagonal like [[. 2 . . .] [. . 2 . .] [. . . 2 .] [. . . . 2] [. . . . .]] then $(I - U)^{-1} = I + U + U^2 + ... U^{n-1} \ (U^n = 0)$. For this example, $(I - U)^{-1} =$ [[ 1 2 4 8 16] [ . 1 2 4 8] [ . . 1 2 4] [ . . . 1 2] [ . . . . 1]] ...


1

Have you plotted the solution using a visualization package or just a cut through the middle of your domain? Sometimes just looking at the solution you're generating will give you an insight into where it's arising from. Based on your description, I'd say there's a mistake in your implementation of your boundary condition. ETA: The last sentence is probably ...


1

(Not exactly an answer to my question, but some empirical experiments which might give guidance to what a worst-case analysis of the L-BFGS iteration of a quadratic function should result in). I tested directly translated algorithm pseudocode implementations of GD and L-BFGS. See the code and all my experiments here. Even in the nicest setting you can think ...


1

Your step two is to solve the original ODE, which doesn't make sense. I'll write out the steps for applying Forward Euler to your second order ODE. Forward Euler solves the first order ODE $$ M \dot{y} = f(y) $$ with the steps $$ \begin{align} M k_1 &= f(y_n) \\ y_{n+1} &= y_n + dt \, k_1 \end{align} $$ Let $v = \dot{u}$. You finite element ...


1

This is speculation, as I do not know the Lipschitz constant or the derivative scales of your simulation. Also, there might be some resonance effect in the interplay of the discrete and continuous parts. But what I would first draw attention to is that the error of RK4 (and any other method) has a V shape in a loglog plot vs. the step size. This is the ...


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