22

This is a really interesting question, and there are a lot of possible explanations. If we are attempting to use a polynomial interpolation, then note that polynomial satisfy the following annoying inequality Given a polynomial $P$ of degree not exceeding $N$ we have $$ |P^{\prime}(x)| \leq \frac{N}{\sqrt{1-x^2}}\max _x |P(x) | $$ for every $x \in (-1,1)$...


20

It's a long-running joke that CFD stands for "colorful fluid dynamics". Nevertheless, it is used -- and useful -- in a wide range of applications. I believe your discontent stems from not sufficiently distinguishing between two interconnected but different steps: creating a mathematical model of a physical process and solving it numerically. Let me ...


20

You should definitely check out Julia. Julia is a programming language which is similar to Python or MATLAB but utilizes a strong type-inference algorithm + JIT in order to optimize code. If types can be fully inferred in a function (which it usually can), then the code compiles in a fully static manner that matches C or Fortran performance. Translating ...


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

What is it you want to achieve? If you want to develop methods/algorithms you might prefer a language that is flexible, and that you are familiar with. As you stated in your question, the Fortran code of your professor was hard to grasp for you, so you re-implemented it in MATLAB. This is very natural way of doing method development: first you implement ...


11

Start simple. Learn Python. I have been paid to write programms for over forty years and I have used all the languages mentioned in other answers (except Julia - I had never heard of it before now). Each language has its strengths, and most have their weaknesses. Like human languages, code is a way of expressing and framing ideas and when you know how to ...


11

Given code that computes a function $f(x)$, automatic differentiation tools produce a code that can compute $f(x)$ and its derivatives at the same time. Solving a differential equation is an entirely different problem and AD doesn't solve differential equations (although AD tools are sometimes useful in connection with PDE constrained optimization.) AD ...


10

OOP has great utility in scientific programming, but the natural way one would write an "OOP" scientific program after education in non-scientific software engineering may be a bad way to tackle the problem. Anecdote: My first semi-serious CFD code as an undergrad, I wrote a finite-volume 2D Euler equation solver with excessive use of OOP. Intuitively, I ...


10

Overview Good question. There is a paper entitled "Improving the accuracy of the matrix differentiation method for arbitrary collocation points" by R. Baltensperger. It's no big deal in my opinion, but it has a point (that already was known before the appearance in 2000): it stresses the importance of an accurate representation of the fact that the ...


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

I'm going to assume since you mention electrodynamics that you're interested in PDEs. You've already mentioned FEniCS in your comment. FEniCS offers a domain-specific language (DSL) called UFL. This language makes it simple to express the weak form of a PDE for discretization via the finite element method. There's another package called Firedrake*, which ...


9

In Mathematica you can do this on multiple levels. First, let's just try to solve it. We need a geometry, an equation and boundary conditions: reg = ImplicitRegion[0 <= x <= 2 && 0 <= y <= 0.5 && ! (x >= 1 && y <= 0.1) && ! (x >= 1 && y >= 0.4), {x,y}]; rp = RegionPlot[reg, AspectRatio -> ...


9

You should specify the eigenvalues you want with which="SM", for example. Check the following snippet. I also changed the solver, since your system is symmetric. import numpy as np from scipy.sparse.linalg import eigsh import matplotlib.pyplot as plt n = 200 h = 2/(n-1) # domain for x and y is [-1, 1] L = np.diag(np.ones(n-1), k=-1) - np.diag(2*np....


9

If your discrete equations represent your problem properly you should have the right eigenvalues, although numerically you can get a small imaginary part when solving the eigenvalue problem. As an example, let us consider the axisymmetric case of a cylindrical well. The Schrödinger equation would be written as $$\frac{\partial^2 \psi}{\partial r^2} + \frac{...


9

The preconditioner for the FDM method that corresponds to the one you outline for the FEM (i.e., the Sylvester-Wathen approach) will still contain the Schur complement of the FDM matrix. The Schur complement of the FDM matrix will, in general, have the same kind of structure it has for the FEM, i.e., it will be spectrally equivalent to the identity operator ...


8

As far as I can tell, this scheme just consists in replacing the uniform finite difference stencil near the boundary with a non-uniform stencil (with at least one point shifted to lie on the boundary). Basically, you take your arbitrarily shaped domain, put it in a box, discretize the box with a uniform grid, throw away all grid points that do not have at ...


8

If you substitute, at least for your analysis, $\frac{\partial u}{\partial x}$ by $u_x$, you can write your system as $$ \begin{bmatrix} 0 & 0 \\ I & I \end{bmatrix} \frac{d}{dt} \begin{bmatrix} p_h(t) \\ u_{x,h}(t) \end{bmatrix}+\begin{bmatrix} -\partial_h & \partial_h \\ -\Delta_h & 0 \end{bmatrix}\begin{bmatrix} p_h(t) \\ u_{x,h}(t) \end{...


8

Answering your last question first, do people actually use FDM for curved boundary nowadays I'd say the answer is no. In the commercial CFD world, 2nd order accurate finite volume schemes are the de-facto industry standard. One of the advantages of FV (and finite element/discontinuous galerkin approaches Jed mentioned) over FD is the much more natural ...


8

Notice that $\hat C_*=(1-r F(h))\hat C_n$, but the sign in front of $r$ is lost when you use $\hat C_*$ inside $\hat C_{n+1}$. I'm looking at p.149 of Numerical Solution of Time-Dependent Advection-Diffusion-Reaction equations by Hundsdorfer and Verwer on Google books. (I'm going to use their signs.) The stability region of Heun's method is $$|g(z)| = \...


8

The problem you have is what is called "stiff": it has two time scales, namely one for the oscillation (which is ${\cal O}(1)$) and one for the damping (which is ${\cal O}(\eta^{-1})$). If $\eta$ is very small, then these time scales are very different. Stiff ODEs are difficult to solve because in order to solve them accurately, you have to resolve the ...


8

When looking at the solution of your system, you will find that almost all entries of $x$ are nonzero although the right-hand side is "sparse". Hence, whatever algorithm you use, it'll have to visit each and every entry at least once, so one wouldn't expect that you can save a lot of time using the sparsity of $b$. Right-hand sides where you can save a lot ...


8

Numerical solution of the advection equation with centered differences in space and forward Euler in time is unconditionally unstable. So the behavior you are seeing is expected. Here is a nice lecture by Gil Strang ( MIT 18.086) where he discusses the instability of this method. He also shows how the simple centered difference method can be modified to ...


8

I think you're probably seeing artifacts that are due to numerical dispersion. In brief, in the discrete case different (spatial) frequencies of a wave function will propagate at different phase/group velocities. This stands in contrast to the continuous case, where the all frequency components travel with exactly the same speed (c). In the first simulation, ...


8

It's easy to derive that equation from Fick's law. You have this diffusion equation as: $$\frac{\partial C}{\partial t} = D \nabla^{2} C$$ The mean square displacement weighted by the concentration profile is defined as: $$\langle r^{2}(t) \rangle = \int_{\Omega} |\vec{r}|^{2} C(\vec{r},t) d^{3} \vec{r}$$ The time-evolution of this mean square ...


8

The difference between finite volumes and finite differences is really more about the form of the equations solved. In typical FV methods, the conservative form is discretised in terms of integrals and fluxes, whereas FD methods generally approximate derivatives in the non-conservative form directly. Maintaining conservation and preserving physical ...


7

What you describe as your time discretization is called the Crank-Nicolson scheme. For nonlinear differential equations it leads, as you have observed, to a nonlinear algebraic system that needs to be solved at each time step. The typical approach is to solve it with a Newton method -- in your case, that requires to solve a nonlinear system in 5 variables. ...


7

Absolutely nothing is "natural" about using procedural programming for finite differences or finite elements. In fact, all widely used software packages for finite element methods that have come up in the past 15 years are object-oriented (e.g., libmesh, openfoam, my own library deal.II, fenics, ...) and the same is true for linear algebra packages (e.g., ...


7

You can convert $\mathbf{b}-\mathbf{x}$ into polar coordinates, and do the dot product in this system. This changes $((\mathbf{b}-\mathbf{x})\cdot\nabla)^n\frac{1}{r}$ to $$\left((\mathbf{b}-\mathbf{x})_r\frac{\partial}{\partial r} + (\mathbf{b}-\mathbf{x})_{\theta}\frac{1}{r}\frac{\partial}{\partial \theta}\right)^n\frac{1}{r}$$ Here, I am using the ...


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