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

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


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


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


9

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


8

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


8

A phenomenological model is based on observations of a system rather than on physical theory. Other physically based models are based on fundamental physical principles such as Newton's laws of motion. Both kinds of models might end up being expressed in the form of mathematical equations and called mathematical models. In practice, models used in many ...


8

Up until a couple of decades ago, science was based on two large pillars. Those were theory and actual physical experiments. It is an exciting time to see a third pillar arise with numerical simulations. In between pure theory and expensive real-world experiments, we can now run simulations! When it comes to these simulations, you may observe two types of ...


7

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


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


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


6

A first step, if you "have never been up in computing", is to read the literature and see what others are doing and have done. The second step is that you will likely learn that what you want to do is not possible today -- at least unless you have access to supercomputers. I suspect that 3 billion particles is possible today, but only if you have access to ...


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


5

The short answer is: it depends. Some journals have a requirement or a guideline to have experimental results that corroborate the numerical simulations. For example, IEEE Transactions on Microwave Theory and Techniques rarely publishes papers that do not have experimental results in one form or another. However, this is very different for IEEE Transactions ...


5

I agree with the spirit of the comments from @user14717. You can burn the candle from both ends here. I suggest splitting your time between small, fundamental problems that can be solved in < 100 lines of code project-based problems that foster knowledge of production code architecture and implementation For 1., your textbook problems are good. I also ...


5

There are basically two methods: You can disrectize the angular part via grid points, or you can discretize it via basis expansion. I will focus on spherical symmetry here, the cylindrical case is quite similar. In the basis expansion approach, one applies the ansatz $$ v(x,t) = \sum_{klm} a_{klm}(t)\,R_{klm}(r)\, Y_{lm}(\theta,\phi) $$ This is inserted ...


5

From a performance view, you are always interested in preserving as much 'structure' in your grid as possible. Computations on a simplex- or a hexaedral mesh, where every cell looks like the next will be more performant, as you do not have to transform from local to global coordinates differently for each cell. Also, you do not have to save the cell ...


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