5

So, you want to invert your matrix $A=\Phi^T\Phi$. For $A$ to be invertible it must not have zero eigenvalues. We can show that $A$ is positive semi-definite as follows. Positive semi-definite means that the eigenvalues of $A$ are $\geq 0$. This is equivalent to showing $y^TAy \geq 0, \forall y \neq 0 $. $$ y^TAy = y^T\Phi^T\Phi{y}=(\Phi{y})^T(\Phi{y}) \geq ...


4

Stability does not necessarily imply accuracy. I'll demonstrate this with a simple scalar ODE $y' = \lambda y$ (known as Dahlquist's test equation). This simple ODE is generally interpreted as a linearisation of a generic ODE $y'=f(y)$ around a given solution point, i.e. it reproduces, at least locally, the behaviour of the more complex ODE you might be ...


4

While $a^{b}$ is usually computed via $\exp(b \log (a))$, in real-life high-quality implementations of functions like pow() the logarithm is computed with extra precision to counter the error magnification effect of $\exp$, as demonstrated in this answer for example. This suggests that we should favor computation via a built-in general exponentiation ...


4

It's not overcomplicated: That's how $x^y$ is actually defined -- via $$ x^y = \exp(y \log x). $$ The point is that for non-integer $y$, it's not at all obvious what $x^y$ actually means: It is not just a product of $x$ with itself, repeated $y$ times. The definition of $x^y$ then is done by pulling the expression back to functions we have previously ...


4

Look up something on Tikhonov regularization, also known as ridge regression in machine learning. This is a standard technique (but I agree that the explanation in that notebook is somewhat poor). Technically speaking, it does not affect the numerical stability of that algorithm, but it modifies the problem to a more well-conditioned one, from $\min \|\Phi \...


4

Think of the simplest case when $\Phi$ is a scalar value. Not well defined: $$ \boldsymbol \theta^\text{ML} = (0^T 0)^{-1}0^T ~ y = \frac{1}{0} 0~y= \frac{0}{0} $$ Well defined: $$ \boldsymbol \theta^\text{ML} = (0^T 0 + \kappa)^{-1}0^T~y =\frac{1}{\kappa} 0 ~y= 0 $$


4

The explanation in the book does not use von Neumann analysis at all but the absolute stability regions and the eigenvalues of the discrete Laplacian operator. For the result you specifically mentioned we use the fact that the maximum eigenvalues is $$ \lambda_m \approx -\frac{4}{h^2} $$ from the expression given. We then want this eigenvalue to lie inside ...


3

I believe Von Neumann's stability analysis would give you the answer here. Consider the heat transfer equation: $$\frac{\partial \mathcal{T}}{\partial t} = \alpha \frac{\partial^{2} \mathcal{T}}{\partial x^{2}}$$ By using Forward Euler time integration and central difference in space discretization: $$\mathcal{T}^{t+\Delta t}_{x} = \mathcal{T}^{t}_{x} + \...


3

There is an entire research field on system theory. The problem is well known but that doesn't mean the answer is simple. Some basics First, let me write your problem in a more common notation. In general, a discrete linear system is described as: $x_{k+1} = Ax_k + Bu_k $ $y_{k} = Cx_k + Du_k $ We call $x$ the internal states and $y$ the observable output. I ...


3

A few things jump out at me from your code as potential problems. You seem to be using Newton's method separately for the $y_1$ and $y_2$ variables. This is not the same as using Newton's method for a nonlinear system involving both $y_1$ and $y_2$. For the full Newton method, you'll be solving a 2x2 linear system at every step. You'll have to calculate not ...


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

Disclaimer, I just look at the problem with $r_1=r_2=r_3=r=1$. But I expect, that one can generalize this approach for different $r_i$. I suggest the following mapping: Project the surfaces of an interior cube onto the surface of your superellipsoide. This divides the surface into 6 parts. Because of the symmetry, I will restrict this now to the mapping of ...


2

Yes. That's all there is to the stability condition. Taking the material properties - shear modulus ($\mu$), bulk modulus ($\kappa$) and density ($\rho$) - into account, the global critical time step is evaluated as the minimum of the critical time step for each element ($\Delta t^e$) $\Delta t^e = CFL * h^e / c_{\kappa}$ where CFL is the Courant-Friedrichs-...


2

It turns out that stable summation of numbers is a topic that is still being researched today -- but you can get the basics by looking up "Kahan's summation algorithm". That said, there really is only a stability issue if you have numbers of widely varying size. In that case, you need to do the summation in a particular order -- intuitively from ...


2

The stability analysis is made considering two close trajectories, if the difference between the state is guaranteed to reduce after one iteration, it means that the method is stable, the numerical errors introduced in one iteration (rounding errros) will not accumulate. The numerical solutions satisfy the equation $$ u_i^{k+1} = r\,(u_{i+1}^k + u_{i+1}^k) + ...


1

Here $u$ is complex so the energy is $u^* u$ where $u^*$ is complex conjugate. Then you must compute $$ (u^* u)_t = u^* u_t + u^*_t u= \frac{i}{2} u^* u_{xx} - \frac{i}{2} u^*_{xx} u = \frac{i}{2}( (u^* u_x)_x - (u^*_x u)_x ) $$ Integrating over $x$ and using zero boundary conditions on $u$ $$ \frac{d}{dt}\int_0^1 u^* u dx = \frac{d}{dt}\int_0^1 |u|^2 dx = 0 ...


1

You are solving the wrong equations with the Newton method, even in the sub-optimal version presented in the code. The BDF method requires the solution of an implicit equation, removing denominators it can be written for a system $\dot U=F(U)$ as $$ 25·U[I]-12·H·F(U[I]) = R = 48·U[I-1] - 36·U[I-2] + 16·U[I-3] - 3·U[I-4] $$ If this were for a scalar equation,...


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