16

Note that $\pi/2$ is represented in double precision format in a way that is not exactly equal to $\pi/2$. It's only accurate to about 15 digits. Thus you're starting every so slightly away from the equilibrium position. Since the equilibrium is unstable, it will eventually start moving.


15

I think the two main points have already been made by Brian and Ertxiem: your initial value is an unstable equilibrium and the fact that your numerical computations are never really exact provides the small perturbation that will make the instability kick in. To give a bit more detail how this plays out, consider your problem in the form of a general ...


13

Atmosphere and ocean have highly-stratified flows in which the Coriolis force is a major source of dynamics. Maintaining geostrophic balance is extremely important and many numerical schemes are intended to be exactly compatible (at least in the absence of topography) to avoid radiating energy in gravity waves. Due to the stratification, limiting vertical ...


13

Computing the term $e^x$ is definitely significantly more expensive than computing a lower-order polynomial -- say $x^4$. But it may be ten to 100 times more expensive at most, not "crazy" expensive. So I suspect that if Matlab takes forever to compute something, then that's because the character of your ODE changes significantly. For example, the presence ...


12

At least one difference is that in a system of ODEs, all the equations are differential, e.g.: $$ \dot{x}=f(x,y)\\ \dot{y}=g(x,y) $$ whereas the definition of DAEs that I'm familiar with includes some non-differential (i.e. algebraic) equations in the set, e.g.: $$ \dot{x}=h(x,y)\\ y=l(x,y) $$ where $l$ is non-trival, and its solution can't be easily ...


10

We should keep in mind that models are just representations of a portion of reality (a narrow portion), therefore a divide by zero error or other mathematical error (negative concentration, for instance) is not necessarily a flaw in the conceptual model, only that the model is used beyond its proper range of application. For example, let's consider a model ...


8

Just simply by being consistent in all of my code? Yes this is the only way. Matlab or any other programming language does not know about units. They only know about numbers. As an example consider incompressible flow. If you set your velocity in m/sec, length in meters (how you generate the grid), pressure in Newton/m^2, kinematic viscosity in m^2/sec, ...


7

The reason the resulting geodesic curve was deviating was because the calculated Christoffel symbol of second kind was incorrect. Using the correct Christoffel symbol : C = Matrix([[(0, -tan(v)), (0,0)],[(sin(v)*cos(v),0),(0, 0)]]) results in the proper output (as displayed in the reference) : Now, I suppose I have to figure out why the calculated ...


7

Differential-algebraic equations (DAE) are equations of the form $F(t,x,x')=0$, with the unknown function being $x(t)$. So in a way are generalizations of ODEs. A nice place to start is here. On the other hand an algebraic differential equation is a totally different thing. The wikipedia page gives an overview, but basically is an equation involving ...


7

I would say that you have, mainly, two methods: Being consistent in all your code, as already suggested in another answer. For that purpose I always keep a table like this one with me, since it might prove really useful. Use nondimensional equations. That way, all my parameters and variables are consistent already. For that purpose I suggest reference 1. ...


6

You should also read about the Method of Manufactured Solutions (PDF) which will show you how to generate analytical solutions to your problem.


6

A spheroid is really just a sphere that has been squashed in the different coordinate directions. So to get a mesh for a spheroid is the same as getting a mesh for a sphere: start with the latter, then do the obvious and trivial transform. So how do you get a mesh for a sphere? If you want quadrilaterals, you start with six quads that cover the surface of ...


6

The restrictive CFL of DG schemes typically comes from the combination of high order accuracy and a compact stencil (see this reference for example). The CFL depends on bounding the variational form in terms of the $L^2$ norm of the solution, which depends on derivative and traces of polynomials. Bounds for each of these quantities (using Bernstein or ...


6

One thing that makes your non-dimensionalized ODE confusing is that you use the same symbols for dimensionalized and nondimensionalized variables, even though they are different variables. Consider writing it in a different way: let there be a new set of variables, $\tilde t, \tilde x,\ldots$, all unitless, defined by $$ \tilde t = n_{\mathrm{time}}t, \...


6

The function $q(e)$ satisfies a first order linear ODE $$ \frac{\mathrm{d}q}{\mathrm{d}e} = \frac{111 e^4+876 e^2+288}{(e^2-1) (121 e^2+304)} q(e), $$ which can be solved very easily by using an integrating factor: $$ q(e) = C_1 \exp\left( \tfrac{111}{121} e-3\operatorname{arctanh}e+\tfrac{10440}{1331 \sqrt{19}} \operatorname{arctan}\left(\tfrac{11}{4 \sqrt{...


5

The equation $\frac{\partial}{\partial z}(\rho u)=0$ only has a unique solution for the product $\rho u$ (namely a constant), but not for each factor separately. You need a different equation to tease them apart. In your context, the equation describes the conservation of momentum. You will need a separate equation for the conservation of mass, which is ...


5

Trace the arc of a ping pong ball you are shooting through the room using a slingshot. The equations that describe this are trivial (gravity acts downward, friction acts in the opposite direction of velocity and is proportional to the square of the speed). You can record the arc on a cell phone and I'm sure there is cell phone software that allows you to ...


5

The implicit Euler method is unconditionally stable alright, but what you are doing is not the implicit Euler method. Rather, what you do is compute where the particle would be at the end of the time step using only 2-particle interactions, compute the location update, and then sum these updates up for all particle interactions. But the implicit Euler method ...


5

Though the following will not necessarily solve your problem, it is highly recommended to implement the ode dimensionless such that all relevant quantities are around 1 to stabilize solvers. The result is to be scaled back accordingly. This will also give you a hint what goes wrong in terms of the size of the coefficients or step size in the discretization. ...


5

I think you are missing a very important and crucial step that lies exactly between the physics and simulation: the mathematical model. In order to model any physics, one has to formulate the mathematical description of the physical phenomenon. Depending on the goals of the simulation, different approximations and assumptions can be made resulting in ...


5

A classic paper for evaluation of the integrals commonly present in computational electromagnetics (EM) is: D. R. Wilton, S. M. Rao, A. W. Glisson, D. H. Schaubert, O. M. Al-Bundak, and C. M. Butler, "Potential integrals for uniform and linear source distributions of polygonal and polyhedral domains," IEEE Trans. Antennas Propag., vol. AP-32, no. 3, pp. 276-...


5

To calculate the reaction forces at a node, Abaqus (or any structural FE code) simply sums the internal forces for all elements attached to that node. The reaction forces are the negative of that sum. For an Abaqus user element, the internal forces for the element are returned from subroutine UEL in the RHS array. The returned stiffness matrix (Jacobian), ...


5

Yes, the problem with mixed boundary conditions is well posed. What's not clear to me is this: Why do you approximate the derivative via the two-sides approximation? Shouldn't it be enough to just the following? $$ \frac{f_6 - f_{6-\Delta x}}{\Delta x} = 0. $$ In your animations, what is the size of $\Delta x$? Your curve looks very smooth, which can be ...


4

Jed Brown described the traditional approach used in the mesoscale and larger scale models. Actually, in microscale many atmospheric models are very close to to traditional CFD codes, use similar finite volume discretizations, similar 3D grids where vertical is treated similarly as horizontal, and so on. Depending on the resolutions even features like ...


4

There are many simple cases that you could look at, e.g. 2D-flow over a backward-facing step. Lid-driven cavity flow. (Inviscid) wave resistance of thin ships or travelling pressure distributions. They should be easy to do in a month. There are also hidden depths (e.g. 3D-turbulence effects, non-linearity, etc) to the problems that could give years of ...


4

Yes! Any reentrant corner has a singularity at the corner (thus the classic L-shaped domain problem). There's lots of literature on this, but I don't know that it's a completely solved problem.


4

Your proposed discretisation appears to be consistent, but wouldn't normally be interpreted as a finite volume discretisation. Indeed, it looks a lot more like a finite difference method with a slightly unusual stencil. Further, this looks like a nonlinear hyperbolic, rather than elliptic PDE, which fits with your desire to use FVM. Generically, you can ...


4

The main advantage is that it reduces the Runge phenomenon and leads to faster convergence rates. It also presents less numerical dispersion and need less nodes per wavelength (see 1 and 2). So, I would say that you would prefer the method for wave propagation scenarios. Regarding software that includes SEM, I am aware of the following: FSELib: Matlab ...


4

The initial assumption was that the initial position was at a stable equilibrium (i.e., a minimum of the potential energy) with zero kinetic energy and the system started moving away from the equilibrium. Since physically it can't happen (if we consider classical mechanics), two things came to my mind: The first one is that maybe the initial position is: ...


4

Look at the components of the forces calculated in your functions. You will probably find they are never exactly zero, because as other answers have said, you can't represent the value of $\pi$ exactly in computer arithmetic, and the routines that calculate trig functions are not exact either. Eventually, the tiny forces (probably of order $10^{-16}$ at ...


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