Hot answers tagged

28

The equation you're solving does not permit right-going solutions, so there is no such thing as a reflecting boundary condition for this equation. If you consider the characteristics, you'll realize that you can only impose a boundary condition at the right boundary. You are trying to impose a homogeneous Dirichlet boundary condition at the left boundary, ...


23

Let's look at the simplest case of Poisson's equation $$-\Delta u = f \tag{1}$$ on a domain $\Omega\subset \mathbb{R}^n$ together with homogeneous Dirichlet conditions $$ u|_{\partial\Omega} = 0 \tag{2}$$ on the boundary $\partial\Omega$ of $\Omega$. We assume for now that $\partial\Omega$ is as smooth as we want (e.g., can be parametrized by a $C^\infty$ ...


19

Pseudo time-stepping, probably better known as pseudo-transient continuation, is the technique of solving for the steady-state solution of time-evolving partial differential equations by setting an initial guess and using a time-stepper to evolve the solution forward. It tends to succeed where standard globalization strategies fail by taking advantage of ...


19

Starting with the advection equation is conservative form, $$ \frac{\partial u}{\partial t} = -\frac{\partial (\boldsymbol{v} u)}{\partial x} + s(x,t) $$ The Crank-Nicolson method consists of a time averaged centered difference. $$\frac{u_{j}^{n+1} - u_{j}^{n}}{\Delta t} = -\boldsymbol{v} \left[ \frac{1-\beta}{2\Delta x} \left( u_{j+1}^{n} - u_{j-1}^{n} \...


18

Yes, it is much more difficult to do so. For the $N$ body problem, all you need to compute are the trajectories $\mathbf x_i(t), i=1\ldots N$ which are just $N$ functions of a single variable. On the other hand, even for a single electron, the solution of the Schroedinger equation is a function $\Psi(x,y,z,t)$, i.e., a function of four variables. For two ...


18

Nothing stops you from doing that technically, but when you integrate by parts you get more flexibility with the solution space in that they need not have $H^2$ regularity (required for the non I.B.P formulation). The linear elements you suggest generally have enforced continuity between elements, and so could not be in $H^2$. The I.B.P formulation ...


18

Short answer: No, you don't have to do integration for certain FEMs. But in your case, you have to do that. Long answer: Let's say $u_h$ is the finite element solution. If you choose piecewise linear polynomial as your basis, then taking $\Delta$ on it will give you an order 1 distribution (think taking derivative on a Heaviside step function), and the ...


16

I think that one of your problems is that (as you observed in your comments) Neumann conditions are not the conditions you are looking for, in the sense that they do not imply the conservation of your quantity. To find the correct condition, rewrite your PDE as $$ \frac{\partial \phi}{\partial t} = \frac{\partial}{\partial x}\left( D\frac{\partial \phi}{\...


15

This is a well-framed question and a very useful thing to understand. Korrok is correct to refer you to von Neumann analysis and LeVeque's book. I can add a bit more to that. I'd like to write a detailed answer, but at the moment I only have time for a short one: With $\alpha=\beta=1/2$, you get a method that is absolutely stable for arbitrarily large ...


15

Excellent answers already on this page, but there is still a (small) missing point. The OP asked: Now, let's say that I have a PDE with higher order derivatives, does that mean that there are many possible variational forms, depending on how I use Green's formula? And they all lead to (different) FEM approximations? Integrating by parts (in the ...


15

I'll begin with a general remark: first-order information (i.e., using only gradients, which encode slope) can only give you directional information: It can tell you that the function value decreases in the search direction, but not for how long. To decide how far to go along the search direction, you need extra information (gradient descent with constant ...


13

Finding the eigenvalues for the Schrödinger equation is really similar to finding the eigenvalues for the wave equation. You start with your differential equation $$\left[-\frac{1}{2}\nabla'^2 + V(r)\right]\psi(\mathbf{r}) = E' \psi(\mathbf{r})$$ where we did the change of variable $(x,y,z) \rightarrow (a_0 x, a_0 y, a_0 z)$, with $a_0 \equiv 1$ Bohr, $E'= ...


12

Adams-Moulton method is significantly more stable. The analogy used when I was taught the difference is the same as extrapolation and interpolation. Interpolation is relatively safe numerically. Extrapolation can blow up if you happen to have an asymptote or some other odd feature. For instance, solving the ode $y'(t) = -y(t)$ with $y(0) = 1$ using ...


12

Yes, this is the standard Aubin-Nitsche (or duality) trick. The idea is to use the fact that $L^2$ is its own dual space to write the $L^2$-norm as an operator norm $$\|u\|_{L^2} = \sup_{\phi\in L^2\setminus\{0\}} \frac{(u,\phi)}{\|\phi\|_{L^2}}.$$ We thus have to estimate $(u-u_h,\phi)$ for arbitrary $\phi\in L^2$. To do that, we "lift" $u-u_h$ to $H^1_0$ ...


12

Suppose that the solution $u$ of the PDE lives in some function space $X$. We'll write the PDE as a bilinear form $A(u, v) = f(v)$ for all $v$ in $X$, where $f$ is some element of the dual space $X^*$. To approximate the solution $u$ of this PDE, we can instead look for some field $u_N$ that lives in a finite-dimensional subspace $V_N$ of $X$. Typically, ...


11

Already one good answer is available here, I just want to highlight some things, dual time stepping scheme uses pseudo time in addition to real time (so in your equations two time parameters will come one real and one pseudo). Real time act as a new dimension to the equations and generally discretised implicitly (refer paper). Real time step size should ...


11

You define a sequence in, say, $C^{\infty}(\Omega)\times C^{\infty}(\Omega)$ by $$ \frac{\text{d}^2u^k}{\text{d}x^2} + \frac{\text{d} v^{k-1}}{\text{d} x} =f\\ \frac{\text{d}^2v^k}{\text{d}x^2} + \frac{\text{d} u^{k-1}}{\text{d} x} =g\\ $$ (plus boundary conditions). It is clear that if this sequence converges, it will be a solution of your original set of ...


11

Curved boundaries are covered in most CFD books, e.g., Chapter 11 of Wesseling or Chapter 8 of Ferziger and Peric. While not a fundamental theoretical problem, the practical complexity of implementing boundary conditions for high-order methods on curved boundaries is a significant reason for interest in more geometrically-flexible methods such as the finite ...


11

In general, you cannot just transfer the same polynomial basis from tetrahedral to quadrilateral elements.1 In particular, the whole point of quadrilateral elements is to work with tensor products of one-dimensional polynomials, which is not possible for tetrahedral elements. There are in fact quadrilateral Raviart-Thomas elements, but their definition is ...


11

Yes, lots of people have. Automatic Jacobian sparsity handling shows up in the second tutorial of DifferentialEquations.jl, where it's able to run sparsity detection on normal Julia code to get the sparse form and perform coloring to then specialize the matrix computations. Then the tutorial ends by showing how you can swap out linear solvers for Newton-...


10

Generally speaking, you'll want to use an implicit method for parabolic equations (the diffusion part) -- explicit schemes for parabolic PDE need to have a very short timestep to be stable. Conversely, for the hyperbolic part (advection) you'll want an explicit method as it's cheaper and doesn't disrupt the symmetry of the linear system you have to solve by ...


10

From a numerical perspective, it's perhaps easiest to discuss the discretizations directly. For the Poisson equation with homogeneous Dirichlet boundary conditions, there is a unique solution for any right-hand side. Once discretized, the equation can be written in the form $Ax = b$, where $A$ is the standard discretization of the 3D Laplacian operator with ...


10

Let me first answer all the questions: What is the theoretical convergence rate for an FFT Poison solver? The theoretical convergence is exponential as long as the solution is sufficiently smooth. How fast should this energy converge? The Hartree energy $E_H$ should converge exponentially for a sufficiently smooth solution. If the solution is less ...


10

To implement your problem in FEniCS, you have to replace the integrals in terms of boundaries by integrals in terms of edges. This introduces jumps/averages in the test functions, which you entirely miss in your implementation. Hence, the system is not invertible and your solution does not look right. Equation (3.3) in Arnold et. al. 2002 gives you a tool to ...


10

The reason is that with the exception of linear problems, if you do a Fourier (or other) decomposition in time, you end up with a significant number of problems that are coupled globally in time. In other words, you have to solve lots of problems on the entire time interval concurrently. That will typically bust your computational or memory budget. The ...


10

The biharmonic equation is the Euler-Lagrange equation of the Laplacian energy $\frac{1}{2} \langle \Delta u,\Delta u \rangle$. A systematic approach to discretize higher order problems is to convert the unconstrained problem to a constrained problem: Minimize $\frac{1}{2} \langle v,v\rangle$ s.t. $\Delta u=v$; that is, \begin{equation} \frac{1}{2} \langle v,...


10

How does this adjoint 'trick' improve the cost of the optimization per iteration in the case where the number of design variables is large? I think about the cost from a linear algebra perspective. (See these notes by Stephen G. Johnson, which I find more intuitive than the Lagrange multiplier approach). The forward approach amounts to solving for ...


10

Here I have an example: x = linspace(-5,5,100); y = linspace(-5,5,100); z = linspace(-5,5,100); [X, Y, Z] = meshgrid(x, y, z); Ex = sin(2*pi/5*Z); Ey = 0*X; Ez = 0*X; [Bx, By, Bz, V] = curl(X, Y, Z, Ex, Ey, Ez); Eplot = 0*x; Bplot = 0*x; for i=1:100 %% Integration-like procedure Eplot(i) = mean(mean(Ex(:,:,i),1),2); Bplot(i) = mean(mean(By(:,:,...


10

This phenomenon is often called "ringing" and plagues methods that are not $L$-stable. This can be seen in this motivating example from Hairer & Wanner (1999) "Stiff differential equations solved by Radau methods". Consider the equation $$ \dot y = -50 (y - \cos t) $$ and apply explicit Euler with time step near the stability limit, implicit midpoint ...


10

The best way to do this is (as you said) to just use the definition of periodic boundary conditions and set up your equations correctly from the start using the fact that $u(0)=u(1)$. In fact, even more strongly, periodic boundary conditions identify $x=0$ with $x=1$. For this reason, you should only have one of these points in your solution domain. An open ...


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