We changed our privacy policy. Read more.

Hot answers tagged

10

I don't think you can avoid using a tolerance for floating-point comparisons. Error due to round-off, discretization, and so on using floating-point numbers is unavoidable. What I typically do to test FEM code I write is: test the stiffness & mass matrices on a single element to make sure I get local element assembly right, compare against a known ...


7

Geoff has already given an excellent overview, but I wanted to provide another real world look on it. In the deal.II project (http://www.dealii.org/) we run some 7,000 tests with every change in the code base, on multiple platforms. The tests Geoff describes are mostly "integration tests", i.e., they run through a significant part of the code base. You need ...


7

You will want to read up on operator splitting methods. In essence, in every "macro time step" you would treat fast processes by doing many "micro time steps" in one half of the algorithm, and then do a single macro time step for the slow processes in the other half. For higher order, you will want to use what's known as "Strang splitting".


6

Off the top of my head, there are a number of things that may be going wrong. You might want to verify the following Are the function values or derivatives you are computing numerically stable? Is the rounding error in their evaluation smooth? You can verify this visually by plotting an extremely small interval around a known root such that the function ...


6

The term that you want to search for is multiple timestepping (see, for instance, [1-3]). [1] http://www.cs.unc.edu/Research/nbody/pubs/external/Berne/tuckerman-berne-rossi91.pdf [2] http://www3.nd.edu/~izaguirr/papers/newM3paper.pdf [3] http://arxiv.org/abs/1307.1167


6

There are so many Runge Kutta methods, including Dormand-Prince 45 Cash-Karp 54 Fehlberge (sic) 78 Is there any comparison between them? Well, sure. Here are some traits to compare: Is the method implicit or explicit? (All of your examples are explicit RK methods.) What is the order of convergence? Are there any embedded error ...


6

There are many different ways to do this. One of the standard is a work-precision plot where you plot the amount of time or function calls that it takes in order to achieve a certain level of accuracy. You can find tons of examples at DiffEqBenchmarks.jl. Generally you slide a timestep or adaptivity tolerances along a window and plot all of the (time,error) ...


5

Since your problem is small, you're probably best off trying fsolve or root. Both of these are interfaces to MINPACK and call HYBRD or HYBRJ. Since calculating a Jacobian matrix for your system shouldn't be hard (either do it by hand, or use your favorite computer algebra system, like SymPy, Sage, Maple, or Mathematica), you should supply a Jacobian matrix. ...


5

Your question really doesn't admit a simple answer- we need to know more specifics about your problem to provide a useful answer. In general, iterative methods can be faster than direct factorization for large sparse systems of equations if the system is reasonably well conditioned or if it is badly conditioned but you have a good preconditioner or if you ...


5

What do people suggest for the linear system of this type? I know about trilinos, petsc, and sundials, but don't know the other alternatives or have exposure to them. Generally speaking PETSc, Trilinos, and KINSOL (from SUNDIALS) are the best-of-breed when it comes to scalability. From an ecosystem standpoint, PETSc seems most flexible, since it does not ...


5

I am a PETSc developer so take my suggestion with a grain of salt, but I would use PETSc because the problem sizes are large enough that execution overhead should be minimal, you can trivially switch between various sparse and dense solvers (1/2 sparse should be treated as dense, but it might pay off to use a sparse solver for 1/9 at your sizes), a suite of ...


5

Daniel Shapero's answer is excellent, but I felt I should add the following: Properly preconditioned iterative methods will almost certainly win here, unless your systems have a very, very special structure. There's a considerable recent literature, unfortunately mostly of the "asymptotic a priori bounds on running time" variety, on solvers for graph ...


5

If your problems are small, you could try a simple solver with dense data structures and use stack allocation (old Fortran codes). IPOPT uses sparse data structures and indeed, it takes a relatively long time to initialize IPOPT if your problems are small. Old Fortran codes with stack allocation are surprisingly fast if you have small problems. If the ...


5

What is typically done for ODEs is to solve at each timestep using two methods of different order. The discrepancies between the results are then used to estimate the error and then decide whether to decrease or increase the timestep length. You can find more information on the subject here. Two examples of these types of schemes are RKF45 and Dormand-...


5

for 1-D transport problem, under implicit method, we use courant number, a dimensionless number to choose appropriate time steps: $\frac{u\Delta t}{\Delta x}\le C_{max}$ $C_{max}$ should be less than 1, if we make it as 1, we get $\Delta t\le \frac{\Delta x}{u}$ The basic ideal is that, at one time step, a particle placed in the problem domain should ...


5

I highly recommend the following read: J.R. Shewchuk, "An Introduction to the Conjugate Gradient Method Without the Agonizing Pain" In short, if the matrix is non-positive definite, there is no guarantee that CG will fail. It might be able to solve it (for some RHSs and certain tolerances), but it is just not supposed to and, probably, converges slower ...


4

Using a line search or trust region method along with Newton's method typically does the trick of finding the zero reliably. Check any advanced numerical methods or optimization book.


4

Can you bound your zero? If you can come up with a good heuristic for the approximate location of the zero and then compute bounds for it over an interval that is not too large, the performance of many methods is improved. Here's some Matlab code I use to perform basic root bracketing for monotonic functions. In terms of interval methods, @Pedro's ...


4

Which one was faster when you tried using Eigen? You may also consider Trilinos if you're using C++. By "generic Laplacian matrix", do you mean a finite difference / finite element discretization of the Laplace operator on some domain, or do you mean the Laplacian of some graph $G$? Whether or not your solver is overwhelmed as you scale up the number of ...


4

It is quite straightforward to demonstrate this for implicit Euler (aka backward Euler) on a scalar example. Consider the initial value problem $\dot{y}(t) = i \alpha y(t), \ y(0) = 1$ with solution $y(t) = \exp(i \alpha t) = \cos(\alpha t) + i \sin(\alpha t)$. The solution is a harmonic oscillation with amplitude 1 and this amplitude does not change ...


4

Many numerical tips and theoretical explanations can be found in this book from Hairer and Wanner: https://www.springer.com/gp/book/9783540566700 In this book, a strategy is described, which uses a time step such that the relative variation of the solution during the first time step is below a certain threshold if you were using explicit Euler (omitting the ...


4

So so many places you have to rewrite. The whole mesh handling (accessing faces and edges from cells, neighbors from cells, ...). Shape functions. Dealing with the question of how the normal vector of a face when seen from one cell matches that on a neighboring cell. You will also likely encounter that 3d problems are always much larger and that solvers, ...


3

For a sparse parallel solver, it's your own responsibility to provide a matrix vector product and a suitable preconditioner. The data for the vector itself should fit into main memory in any case. If the matrix has at most a fixed small number of non-zero elements (<20) per column or row, then the same is also true for the matrix itself. In this case, an ...


3

You can try using SCS, either the direct or indirect solver. SCS uses first order methods, and hence may be able to solve larger problems than second order solvers such as SDPT3, SeDuMi, MOSEK, etc. On the downside, given that it is a first order solver, it can be very slow - but very slow may be better than not at all. Paper: http://web.stanford.edu/~...


3

No, ABAQUS is not a general-purpose PDE solver. It solves the PDE for a number of pre-defined applications such as structural analysis or heat transfer in solids. Of course many physical phenomena are governed by the same equation as heat transfer so you can solve those problems by treating them as heat transfer problems in ABAQUS.


3

A few additional points I would like to add to other answers. Corner test cases should be part of the regression test suite - ill conditioned problems, ill conditioned - bad aspect ratios, orthotropic and un-isotropic material properties, improperly constrained models. Make sure reaction forces match. Sturm checks for eigenvalue problems especially when ...


3

I cant think of any truly generic ones except the Lax-Friedrichs flux, which is very dissipative. The next simplest would be the Rusanov flux which needs knowledge of minimum and maximum speeds arising in the Riemann problem, which of course depends on the particular form of the flux function. While this is also called a central flux, it can be interpreted ...


3

If you want a ten-line solution that is decently fast and stable, you can implement yourself the structured doubling algorithm: set up the coupled iteration $A_0 = A, G_0 = G = BR^{-1}B^T, H_0 = Q$ While $\frac{\|H_{k+1}-H_k\|}{\|H_{k+1}\|} \geq \varepsilon$: $\quad \quad A_{k+1} = A_k(I+G_kH_k)^{-1}A_k$ $\quad \quad G_{k+1} = G_k + A_k(I+G_kH_k)^{-1}G_kA_k^...


3

It depends on the selected level of abstraction and chosen classification. The question would be in the usability of this abstraction and the chosen classification of problems. If you are allowed to assume that an objective function can mean a boolean function that is true if and only if the found solution satisfies the given constraints (rules of the ...


3

Every iterative solver -- Jacobi, SSOR, CG, etc -- starts with an initial approximation. One often just uses the zero vector, but there is nothing wrong with using the solution of the previous time step. In fact, extrapolating from previous time steps to the current one is an even better idea -- one the authors apparently missed! For some iterative solvers, ...


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