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

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


7

You might want to check out DifferentialEquations.jl. It supports ODEs, PDEs, stochastic equations, delay equations, and basically everything else. It also has really good automatic sparsity detection and ability to work on GPU for big systems. DiffEqOperators.jl is the submodule which has the automated finite difference operators (with lazy stencil ...


6

Finite difference approximations of the Jacobian are really only good if the step lengths are chosen appropriately for each coordinate. But a black-box solver like CVODE has no way of knowing what these step lengths should be, and so has to use heuristics to choose them. This may or may not work. You are almost always better off if you provide an ...


5

The short answer is that you need $$\phi_{-1} = \phi_0$$ $$\phi_N = \phi_{N-1}$$ to impose $\nabla\phi=0$. A quick check by making the following change if idx == -1: idx = 0 elif idx == N: idx = N-1 in the code, you have posted shows that the average $\phi$ remains constant up to 14 decimal places. To see why this is the correct boundary condition ...


5

I'm afraid the method only works to compute derivatives of real-valued functions (of which you happen to have an implementation that also works on complex values).


5

I think that the main problem might be with the solver you are using. The Hamiltonian (matrix) in this case is Hermitian, it is even symmetric since it is purely real. You could use eigh instead of eig to take advantage of this. Furthermore, you are not removing only the first and last points but intervals of size 1 at each end. Following, I show you a ...


4

You simply have a bug in your code. The flux is $\frac{1}{2} u^2$ and not $\frac{1}{4} u^2$.


4

The two quoted approximations to the mixed (second) derivative are simply different formulas. There is no direct way to arrive from one at the other. The first formula evaluates the derivative on the points $\Big((x-h,y-k), (x-h,y+k), (x+h,y-k),(x+h,y+k)\Big)$, the second formula uses another stencil that involves points that are also used by the first ...


4

Yes, that finite difference is correct. You can obtain it using a finite difference in each direction for the Laplace operator in each coordinate. \begin{align} \nabla^2 u =& \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} + \frac{\partial^2 u}{\partial z^2}\\ \approx& \frac{1}{h^2}[u(x + h, y, z) - 2u(x, y, z) + u(x -h, y, z)]\\...


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

Given that you are using center difference formula to to get second order derivatives in your domain. The common practice at the ends is to use forward and backward difference formula ( start and end ) . Again this strictly depends on the nature of your problem and what you consider as reasonable approximation


3

Crank-Nicholson is usually used for linear ODE but can also be extended to a nonlinear ODE $$\dot z = f(z)$$ like so: $$\frac{z_{n + 1} - z_n}{\delta t} = \frac{f(z_{n + 1}) + f(z_n)}{2}.$$ Another similar alternative is the implicit midpoint rule: $$\frac{z_{n + 1} - z_n}{\delta t} = f\left(\frac{z_{n + 1} + z_n}{2}\right).$$ The Crank-Nicholson and ...


3

You did nothing wrong, the numbers just are what they are. Plotting the error for the backwards difference gives the plot which is as expected for a method of order one with step size about $h=0.01$. For the second order methods one expects an error of the magnitude $h^2=10^{-4}$, the plot below confirms this


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

The notation $$\frac{\partial U}{\partial \eta}$$ means usually $$\eta \cdot \nabla U$$. This is correct even if the domain is the interval $[a,b]$. The normal vector on the interval $[a,b]$ @a is $\eta=-1$ and @b $\eta= 1$ both pointing outwards of the domain. Hence in 1D $\frac{\partial U}{\partial \eta}$ means $$\eta\cdot\nabla U=\eta \frac{dU}{dx}$$.


3

If your advection problem had Dirichlet of Neumann boundary conditions, the linear system would be tridiagonal and you could apply the Thomas algorithm. With periodic boundary conditions, however, we lose this. If c(x) is a constant independent of x, the matrix would be circulant and linear systems could be solved efficiently using FFTs. An even better ...


3

$c$ depending on time is not the issue. You will use an RK scheme which takes care of this. The issue is $c$ is discontinuous in $x$. I recommend SBP-SAT schemes for this. (1) Derive an energy equation at PDE level. (2) Search literature for SBP-SAT schemes which enforce interface conditions via SAT penalty terms, which are designed to mimic the energy ...


3

I think you can understand this using the concept of a modified equation. As you have shown, your discretization $$ \frac{u_i - u_{i-1}}{h} + f(x) = 0 $$ is an approximation of your differential equation $u_x + f(x) = 0$ with accuracy $\mathcal{O}(h)$. Now consider the Taylor expansion of your exact solution $$ u(x_{i-1}) = u(x_i) - h u_x(x_i) + \frac{1}{2} ...


3

You lose convergence order, or in the worst case convergence altogether. You can try this out: Take $$ f(x) = \begin{cases} 0 & \text{if $x<0$} \\ x^2 & \text{if $x\ge 0$}.\end{cases} $$ The function is differentiable, but not twice differentiable at $x=0$. It's exact derivative is $f'(0)=0$. Now compute the (second-order) symmetric finite ...


3

Generally the step from compressible Euler equations to the Navier-Stokes equations is not that hard, at least the coding part. If you want to implement it with an explicit scheme you have to consider the severe time step restriction of the parabolic contributions. One tricky part, at least for a consistent FV implementation, is the calculation of the ...


3

Regarding performance, Python is definitely the bottleneck. I have experienced the same issue with a 2D Euler code I had developed, even with vectorised operations everywhere possible. It was actually even worse, as I was using solve_ivp time schemes which reallocated memory at every step... You can try and profile your code to see where the bottlenecks are. ...


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

Indeed, the roundoff error is much lower with the complex step approach, because all the derivative terms are handled in the imaginary part, without direct interaction with the real part, as opposed to a finite difference calculation. The advantage of complex step is that you can take a very low perturbation size $h$ (I usually use $10^{-50}$), whereas you ...


3

Probably not what OP was waiting for, but I think it could be pretty instructive and useful. FEM codes use a much different approach to build the so-called stiffness matrix. In practice, they loop over elements and compute for each element small matrices (in your case if you use linear elements 3by3) which are distributed to the right entries of the global ...


3

There are some unknowns in what you are doing but for simplicity, suppose we want to find $u(t)$ as discrete times $t_1, t_2, \cdots, t_n$. Let $\textbf{F} = [F(t_1), F(t_2), \cdots, F(t_n)]^T$ and $\textbf{u} = [u(t_1), u(t_2), \cdots, u(t_n)]^T$ be column vectors representing $F$ and $u$ evaluated at the desired times. From your problem statement, you wish ...


3

The choice of finite-difference scheme depends on several factors, such as the smoothness of your data, how uniformly-spaced the data actually is, etc. You may also want to consider just how accurate your velocity estimate actually needs to be. For example, if there are large error bars in the trajectory then it probably makes little sense to use a high-...


2

Looking briefly to your implementation, it seems to me that you may not discretize correctly the (independent) variable $S$. The coefficient of your PDE are depending on $S$, so you should evaluate them using correct (discretizied) values. Choosing smaller values od $\Delta S$ requires more discrete values of $S$, you can not choose a fixed number of "...


2

It depends on the equation you are solving. For example, consider the following heat equation with a source term: $$ \frac{dT}{dt} - \nabla^2 T = S$$ If $S \neq 0$, then you have a source term that changes the value of the temperature as time evolves. Consequently, even if your initial condition is $T(\textbf{x},t)=0$, your temperature will increase because ...


2

The BOUT++ project http://boutproject.github.io offers a set of tools for finite-difference solution of systems of PDEs, primarily targeting fluid dynamics and plasma physics but not limited to those application areas. Equations in human-readable form are automatically discretized in space and integrated in time, using efficient parallel solvers.


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