Hot answers tagged

8

The feature that you demand is called event location in Matlab ODE solvers pack, or rootfinding in SUNDIALS solvers suite terminology. Essentially this feature allows to stop integration exactly at the point where some vector function of free and dependent variables has a root. Namely, for system $$ \frac{d\boldsymbol y}{dt} = f(t, \boldsymbol y), \quad \...


5

Hashing floating-point numbers can indeed lead to weird results, especially if the node positions can be perturbed by some small amount or if there are denormalized values. You included the Python tag, so I assume that's what you're using and that you have scipy. A quick and dirty solution would be to construct a kd-tree of the fine grid points, then for ...


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


4

You achieve this with equidistribution to a mesh density function $\rho(x)$. If you consider $x$ as a continuous map from $\xi \in [0,1]$ into your domain $[0,L]$, then the statement $x$ equidistributes $\rho$ is equivalent to $$\int_0^{x(\xi)} \rho(x') \mathrm{d}x' = \xi \theta\,,$$ where $$\theta = \int_0^L \rho(x) \mathrm{d}x\,.$$ Not sure if this has a ...


4

In this situation, you need an integrator with dense output, which means it has a capability to calculate the solution to the ODE at any point in time between the time steps taken by the integrator. dop853 (from Hairer and Wanner's eighth-order Dormand-Prince integrator) has this capability, in addition to adaptive time-stepping. Likely what happens is that ...


4

Your numerical solution is probably just getting more accurate as you increase the number of grid points. Do you know or have you tried to derive the analytic (exact) solution for this problem? By looking at your plots, it seems like the exact solution has a shock (discontinuity) occurring at around x = 4.5, and the numerical method is resolving it with ...


3

Let me formulate my remark as an answer (and make it more precise): My experience and knowledge is that to use collocated grids with any standard (lower or higher order) methods for this problem, you must modify (stabilize) your Navier-Stokes equation. it means with some additional effort. This needs not to be necessary when using staggered grids. The ...


3

This is not a good way to do modern programming for many reasons. First of all, as you pointed out, this kind of code is hard to read and maintain. Secondly, this tends to be done in old versions of Fortran for reasons which are no longer a problem. Namely, old versions of Fortran didn't have ways to "bundle together" objects. If you made different arrays ...


3

In essence, you are asking whether you can enumerate the integer lattice sites within your domain from $1$ to $N$ in such a way that accessing the east/west/north/south neighbors of a location $n$ requires accessing positions $n_e,n_w,n_n,n_s$ so that the distance between index $n$ and these four indices is minimal. This problem is equivalent to shuffling ...


3

The solution is https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html From the documentation : ‘RK45’ or ‘RK23’ method for non-stiff problems and ‘Radau’ or ‘BDF’ for stiff problems The documentation taken from scipy: scipy.integrate.solve_ivp(fun, t_span, y0, method='RK45', t_eval=None, dense_output=False, events=None, ...


2

I ended up raising an exception in the deriv() function (and catching it in the loop) because it turned out to be the only reliable way to end integration in every single case and for every ODE. (While I considered this option already in my original question, it seems there is no other way. If anybody presents a better option I will gladly mark her/his ...


2

SRM is a standard protocol, so it is suppose to work for all storage implementations. The only things that differ are the protocol implementations. In addition to bestman2 and dcache srm tools, there is also another popular implementation, lcg-util (lcg-cp, lcg-ls...). I don't know why srmping and srm-ping do not work the same. To be honest, I have ...


2

Writing C++ code from the ground up for adaptive mesh refinement (as part of a PDE solver) is a relatively complicated endeavor and can easily involve thousands of lines of code for even simple problems. Based on the comments it sounds like you may be fairly new to programming (for example unsure how to implement a binary tree). Of course there is nothing ...


2

If your data are regularly spaced, you do not need an interpolation procedure, you can directly plot a contour through a contour or a pcolor command. That why I don't understand why you need to initially interpolate. Contouring algorithms often use marching cube method (and its variations) to render a contour plot. This kind of algorithm only work with a ...


2

A simpler answer, IMHO, is to solve your system using PyDSTool. It has built-in support for accurate event finding. In fact, it's easier to set up and more efficient (including faster) than Matlab's approach if you use the C-based solvers. I am the author of this package. A tutorial using events is given here.


2

It seems like you are inventing Barnes-Hut-type algorithm, which is a fundamental accelerated algorithm for n-body simulations. It follows a similar logic: you combine masses on a grid. But, it is done in a more elaborate fashion: when particles are close, you still use direct calculations (no combination), otherwise, you would lose accuracy. multipole ...


2

In the section about adaptive Methods Chapter 16. in "Chebyshev and Fourier Spectral Methods" from John P. Boyd several different coordinate transformations together with their application in different publications are presented. I have not used any of those transformation myself but they can serve as a starting point for your problem. y is the physical (...


2

Use the chain rule to get the derivative on the non-uniform grid, $\frac{dy}{d\zeta}$: $$ \begin{align} \frac{dy}{d\zeta} &= \frac{dy}{dx} \frac{dx}{d\zeta}\\ &= \frac{dy}{dx} \cdot 2 \zeta \\ &= \frac{dy}{dx} \cdot 2 \sqrt{x} \\ \end{align} $$ So if you have a way of computing $\frac{dy}{dx}$ on the uniform grid, you can simply scale it by $2\...


2

Let's assume the following heat transient heat transfer equation in 1D : $$ \frac{\partial T} {\partial t} = \alpha \frac{\partial^2 T}{\partial x^2} $$ If we take its finite difference approximation using an explicit time scheme we obtain : $$ \frac{T_i^{t+\Delta t} - T_i^{t}}{\Delta t} = \frac{\alpha}{(\Delta x)^2} (T_{i-1}^{t}-2T_i^{t}+T_{i+1}^{t}) $...


2

You can probably speed things up a little bit by storing the array in 4x4 square subarrays so that each of them fit in cache line (64 bytes = 4x4 32-bit integers). This changes the probability distribution of number of memory accesses from [0,0,1] (mean number of accesses = 3) (for 1,2,or 3 accesses) to [1/4,1/2,1/4] (mean number of accesses = 2). Maybe ...


2

You've almost certainly got to reduce the timestep to maintain stability due to the CFL condition imposed by your explicit timestepping method choice. That being said, for benchmarking purposes, I'd be inclined to pose the problem as computing to a particular final time and want to compare that to the quality of the answers and times to solution on other ...


2

The best choice for a numerical grid is the one that will most accurately approximate the solution to your problem (without being too computationally expensive). But beyond that the specific features will depend heavily on the type of problem you are trying to solve. A grid might be aesthetically pleasing because it cleverly exploits some symmetry of the ...


2

The easiest thing is to ensure a consistent (although arbitrary) tie-breaking scheme. If your nodes/vertices indexed, this usually means preferring the split edge with the lowest index of its lowest vertex. If the edges share the same lowest index vertex, then check the other vertex and choose the one with the lowest index. So in your situation, imaging the ...


2

Assuming that the arrays are passed already sorted (which is a reasonable assumption since you are starting from two 1D grids), I have a solution which is $\mathcal{O}(n+m)$ where $n$ and $m$ are the respective sizes of the arrays. I wrote the code in MATLAB/Octave, but should be easy to translate to Python. % Compare two arrays for overlaps % Test 1 - Easy %...


1

As pointed out in the comments, the CFL condition dictates how big timestep can be. Therefore as we progressively refine in space, we would be required to take smaller time steps if an explicit temporal discretization is employed. Regardless of the smaller time steps, it would still be acceptable to get a time-averaged quantity to compare across the various ...


1

I think the problem is that you've lost the topology upon the first rasterization: | 0 | 1 | 1 | 0.5 | 0 | could be =============== | 0 | 1 | 1 | 0.5 | 0 | or it could be ============ === | 0 | 1 | 1 | 0.5 | 0 | So the answer is no. Unless you demand that two objects cannot end in the same grid point, then you ...


1

I would second the opinion expressed in the context. For that small problem and limited usage, you don't need a library. Generation of a structured grid and retrieving points can be coded up in at most several screens of code. However, I would point out for you ViennaGrid library, which can be used and actually provides STL iterators. In addition, the ...


1

For scattered data at points with no structure, try inverse-distance-weighted-idw-interpolation-with-python. This combines scipy's fast K-d_trees with inverse-distance aka radial kernels: $\qquad interpol( x ) = \sum w_i \, f_i$, sum over say 10 points nearest $x$ $\qquad \qquad \qquad \qquad w_i = {1 \over |x - x_i|^p}$ normalized so $\sum w_i = 1$ . It's ...


1

Your data is already in the right shape, then you don't need to create a meshgrid. See the code below import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D import numpy as np I = np.array([ [10.55, 0., 0.], [0., 0., 0.01], [0.2, -0.1, 3.33], [0., 2.14, 0.], [0., 3.80, 0.], [9.02, 0., 0.]]) t = np.array([ [0., ...


1

The dataset is huge. So I would like to be more memory-efficient. The main point to consider is the following: […] most of the cells will be “unactivated”, i.e most of the elements of the matrix will store an empty list (no signal). This sounds as if an array of lists (or variably sized arrays) is just fine. In this case, the crucial factor for memory ...


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