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

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


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


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


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


2

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


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

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


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


1

It is useful, as a first step, to analyze method-of-lines (MOL) treating the time integration as exact. For example, suppose we are solving the 1D advection equation $ \frac{\partial{n}}{\partial{t}} =-c \frac{\partial{n}}{\partial{x}}, $ where $c$ is the advection speed. Then, for example, for central difference, we'd have $ \frac{d}{dt} n_i = - \frac{c}{2 ...


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