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


6

The problem with ABC in the form of derivative operators at the boundary is that the exact ABCs are non-local in space and often in time, which makes them hard to use in numerical simulations. You can derive approximate ABCS based on series expansions, but this can be tedious and it may be difficult to get stable schemes. This approach was made famous by a ...


5

It sounds as if you're running a time-dependent linear elasticity simulation, right? Most likely, you're running an "Explicit" time-stepping scheme, which means that all of your information at time $t_{i+1}$ is computed entirely using information from time $t_i$ (or $t_{i-1}$...). One consequence of such a scheme is that in order to be a stable time ...


5

I think that you have two problems: your timestep is too big (the CFL condition is not satisfied); and you are not updating the values for $\phi=2\pi$. I am not sure about the CFL condition for polar coordinates, but I believed that it should be $$2c\frac{\Delta t}{\Delta r \Delta \phi} \leq C_\max \enspace ,$$ what implies that $\Delta t < \frac{\...


4

This may not be your definitive solution, but what you need in any case is to change the discretization method that replace the continuous Eikonal equation with some discrete numerical scheme on the grid. The classical methods like fast marching method etc. are based mostly on the Rouy-Tourin scheme. In the work of Gibou and Min you can find the ...


4

Higher order methods often have a smaller radius of convergence, i.e., they require smaller time steps. In your context, this means that they require a smaller CFL number, often significantly smaller than 1. Did you take this into account? As for the question of whether there are methods that can deal with CFL numbers much larger than one: This really doesn'...


4

The most commonly used explicit ODE solver in structural analysis is the central difference method. Because it is explicit, the solution becomes unstable if the time step is larger than a so-called critical time step. The calculation of this critical time step is straightforward and can be found in many references (e.g. page 808 in this text by Bathe ) and ...


3

I was going to write a comment, but the equation seems to view better in answers.. I assume Von Neumann analysis is the proper approach to derive this equation, but a coordinate transformation from the cartesian CFL condition (I took from wikipedia) is not somehow equivalent? Specifically: \begin{equation} \Delta t \sum_{i=1}^3 \frac{u_i}{\Delta x_i} = ...


3

High order RK methods work fine for a large number of wave propagation (an elastic example using a DG discretization). The procedure is usually to discretize the spatial derivatives in the first order form of the wave equation. This is called semi-discretization and it results in a system of ODEs involving a vector of solution coefficients $U$ (for both $u$...


3

The stress-velocity formulation has been used extensively in DG context on account of the fact that it can lead to a first order system form of the elastic equation. The latter proved to lend itself more easily to a DG formulation. However, recently a DG formulation of the elastic equation in the second-order form has been proposed, which you may find ...


3

For the heat equation imagine a rod that is heated on the left and cooled on the right. Now imagine that instead of a constant prescribed temperature on the left what we want is the heat to steadily rise and fall. This would be an example of a time-dependent boundary condition. This would be written as: \begin{align} \frac{\partial{u}}{\partial{t}} = c\...


3

The problem is that odeint solves first-order ODEs of the form $$\dot u(t) = Au(t),$$ but the wave equation is (after discretization in space) a second-order ODE of the form $$\ddot{u}(t) = Au(t).$$ (And in fact, your plot looks quite reasonable for the solution to a parabolic equation.) To apply a black-box ODE solver (which is really not such a great idea,...


3

There exist Absorbing Boundary Conditions for the wave equation that are stable and that go up to any order of accuracy (limited only by the accuracy of discretization of your model), so that they are good competition for PMLs. A commonly used ABC of the Enquist-Majda type mentioned above, which is of the second order, is: $u_{tt}-u_{xt}-\frac{1}{2}u_{yy}=...


3

I have not studied this particular system before, and I'm sure that someone who has could say much more than I will. I don't think there is any reason to expect that discretizing the equations in the form you have them will lead to highly accurate conservation of total energy. In this form, conservation of total energy depends on cancellation of certain ...


3

You're looking for waveguide port boundary conditions. I think the most accessible treatment is within Jin & Riley's Finite Element Analysis of Antennas and Arrays, Chapter 5. It's available on Amazon, see https://www.amazon.com/dp/0470401281/. A lot of these formulations were first introduced by Jin-Fa Lee, you can find his works in IEEE Microwave ...


3

There is no reason to believe that two random fields with the same arithmetic mean would yield solutions that have anything to do with each other. In fact, for the case you consider, one might imagine that maybe the harmonic mean is actually a better indicator, but even that is unclear -- it could also be the geometric mean. Apart from this, you have to ...


3

Regarding the boundary conditions: Don't be fooled by Wikipedia. Yes, the scenario in the picture suggests an absorption at the boundaries, and yes, one could use absorbing boundary conditions in order to reproduce that numerically. In the simple case of a wavepacket these are readily available, because in the end, there exists an analytical solution for the ...


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


2

My hunch is that you may be prescribing boundary data at both ends of the domain (which leads to oscillations) or incorrectly specifying the boundary data. For a hyperbolic problem, you need to respect the flow of information in your numerical method. For example, consider the linear hyperbolic equation $$ \mathbf u_t = \mathbf A \mathbf u_x $$ with ...


2

A very simple example #include <iostream> #include <cmath> // sin #include <fstream> //ofstream #include <sstream> //stringstream #include <string> //string #include <cstring> //c_str #define PI 3.14159265 // Modifies the elements in an array to five an evenly spaced number of points // Along the specified interval (...


2

Since the wave equation involves a second derivative in time, two initial conditions are necessary: initial displacement $U(x,0)$ and initial velocity $U_t(x,0)$. I don't see the initial velocity condition specified in your question. This leads me to suspect that there is something wrong with the tutorial you linked to... As I suspected, the code in the ...


2

I agree with the suggestion of starting with a simple problem and with the elastic solution. Probably the simplest wave problem is the 1D, infinite bar/string. The analytical solution to this problem is well-known, e.g. http://mathworld.wolfram.com/dAlembertsSolution.html You can model this with a strip of elements in the x-direction. You want zero-stress ...


2

This is not something I've worked on, so I do not know if there is a recommended approach to this, only how I would approach this problem. I'm assuming you can evaluate $C(x) = (\ln(S(x)))_x$ wherever you want. Since $\phi$ is prescribed at $0$ and $L$, I would discretise on $x_i = L\frac{i}{N+1}$ for $i$ from $1$ to $N$. Then the finite difference ...


2

I understand that you want to solve the differential equation $$\nabla^2 u - \frac{1}{c^2}\frac{\partial^2 u}{\partial t^2} = f(t)\delta(\mathbf{r}) \enspace ,$$ for an unbounded domain, $\mathbf{r}$ is the position vector, and $$f(t) = -\frac{2(t-t_0)}{T^2} e^{-\frac{\left( t-t_0\right)^2}{T^2}} \enspace ,$$ although, we don't need to write this ...


2

You can make up as many solutions to this equation as you want using the Method of Manufactured Solutions (PDF). Is there a reason you're focused on this particular case only, or could you use any solution/forcing-function pair?


2

I believe it is OK to use the PML the first way you described. It's not necessary to put the source at the center, otherwise PML is not better than many other absorbing boundary conditions. The PML even doesn't need to be a spherical shell, which is just a choice of convenience.


2

You are asking very complicated questions for which there are likely no answers that can be rigorously proven. If you go back for a second and ask the same question for the solution of the Laplace equation, then you are asking whether, given that $u\in H^{1+\delta}$, you can find an algorithm that constructs a mesh with $N$ cells so that $\|u-u_h\|_{L_2} = {...


2

The answer is very simple: you provide the code with geometric information, i.e. nodal coordinates (which in your case are expressed in metre), and not only topological information, i.e. how the mesh is laid out and connected to the nodes. So depending on the nodal coordinates you provide, you can have different domain "sizes". E.g. if you have four nodes ...


2

There are two problems with your code: 1) You are mixing $i$ and $n$ in the following line: SUM = SUM+Ai*cos((i-0.5)*pi*xi/L)*sin((n-0.5)*pi*ti/(L/c)) It should be: SUM = SUM+Ai*cos((i-0.5)*pi*xi/L)*sin((i-0.5)*pi*ti/(L/c)) 2) You have to normalize the rod displacement. The coefficients are then: Ai = 2*sin((i-0.5)*a*pi/L)/(pi*(i-0.5)**2*pi*a/L) ...


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