Stefano M
  • Member for 9 years, 6 months
  • Last seen more than a month ago
When should log1p and expm1 be used?
Accepted answer
29 votes

We all know that \begin{equation} \exp(x) = \sum_{n=0}^\infty \frac{x^n}{n!} = 1 + x + \frac12 x^2 + \dots \end{equation} implies that for $|x| \ll 1$, we have $\exp(x) \approx 1 + x$. This means ...

View answer
What is the purpose of using integration by parts in deriving a weak form for FEM discretization?
15 votes

Excellent answers already on this page, but there is still a (small) missing point. The OP asked: Now, let's say that I have a PDE with higher order derivatives, does that mean that there are ...

View answer
What are the symptoms of ill-conditioning when using direct methods?
Accepted answer
13 votes

When is a matrix ill conditioned? It depends on the accuracy of the solution you are looking for, as much as "beauty is in the eye of the beholder"... May be your question should better rephrased as ...

View answer
Permute a matrix in-place in numpy
12 votes

According to the docs, there is no in-place permutation method in numpy, something like ndarray.sort. So your options are (assuming that M is a $N\times N$ matrix and p the permutation vector) ...

View answer
solve $xA=b$ for $x$ using LAPACK and BLAS
Accepted answer
10 votes

Trivial answer for square $A$: use dgesvx which solves also for $A^T x = b$ when TRANS = 'T'. Please note that with BLAS or LAPACK you hardly have to transpose (swapping elements in memory) a matrix:...

View answer
Understanding wall time jitter in MATLAB computations
Accepted answer
10 votes

There are several issues in your code. There is no semi-colon after the Y1 = Y0+Y0*(I-A*Y0) statement. In fact you are timing the screen output of Y1, not the computational time The iteration should ...

View answer
Loop optimization with f2py, Cython and Numba
Accepted answer
9 votes

I think that the problem is linked to the way in which f2py generates the fortran interface: the argument to fortranrun.f2py should be stored as a F_CONTIGUOUS array, otherwise the interface will ...

View answer
Using multiple languages in scientific codes
9 votes

I think that the "two language approach" is sound and I feel very comfortable in using it. When you start a new project from scratch you never know beforehand which will be the critical code sections ...

View answer
Data structures for finite volume code: Arrays vs Classes
Accepted answer
9 votes

Simple answer: in modern python every data type is a class, so formally there is no difference between the two solutions you proposed. (Please remember to use new-style classes: classic classes are ...

View answer
Symmetric boundary condition
Accepted answer
8 votes

For the sake of simplicity let me slightly change your notation. Let $u, v, w$ be the components in the $x, y, z$ directions of a true (polar) vector, and $y=0$ the symmetry plane. For a true vector, ...

View answer
LAPACK - singular matrices - what does the positive integer info mean?
Accepted answer
8 votes

Short answer: nothing more than $U_{ii} = 0$, i.e. that your computed $U$ factorization is exactly singular. xGETRF is not safe as a rank revealing factorization, so I would not draw any conclusion, ...

View answer
Positive semi-definiteness of a (symmetric) matrix
8 votes

According to MathWorld a matrix $A \in \mathbb{R}^{n \times n}$ is positive definite iff $$ (x^T A x) > 0 $$ for all non zero vectors $x\in\mathbb{R}^n$. It is trivial to obtain that $$ x^T\,A\,x =...

View answer
Condition number of A'A and AA' formulations
Accepted answer
8 votes

If $A\in\mathbb{R}^{N\times M}$ with $N<M$, then $$ \mathop{\mathrm{rank}}(A^TA) = \mathop{\mathrm{rank}}(AA^T) = \mathop{\mathrm{rank}}(A) \leq N < M $$ so that $A^TA \in \mathbb{R}^{M\times M}...

View answer
Is there an efficient way to form this block matrix with numpy or scipy?
7 votes

The code proposed by the OP can indeed made be more efficient, mainly by noting the fact that to form the sequence $A^i B$, with $i=0\,\dots,N$ you do not have to compute $A^i$ at each step, but you ...

View answer
Finite Element Method vs Extended Finite Element Method (FEM vs XFEM)
7 votes

Both Mike's answer and Jed's one describe well the XFEM/FEM dichotomy and correctly point out that the most important area of application is 3D Fracture Mechanics, where you have a crack, i.e. a ...

View answer
What is the difference between implicit FEM and explicit FEM?
Accepted answer
7 votes

The FEM method for transient problems typically uses the method of lines, i.e. the spatial discretization is decoupled from the time discretization: \begin{equation} u^h(x,t) = \mathbf{\Phi}(x)^T \, \...

View answer
Stress due to the mismatch of thermal expansion coefficients of two different attached materials in COMSOL
Accepted answer
6 votes

Thermal stresses are self stresses that arises in two main cases. If one imposes displacement continuity at the interface between two materials with different thermal expansion subjected to a uniform ...

View answer
debugging a rotation matrix for elastic constants
Accepted answer
6 votes

You are facing one of the most tedious and error prone aspect of elasticity theory: change of reference frame in engineering (or Voigt) notation. Recap theory If $\boldsymbol{\epsilon}$ and $\...

View answer
How should I report profiling/timing information about my code?
Accepted answer
6 votes

Total running time (wall clock) is the only metric that matters in industry or real life applications: this figure should never be omitted, even if embarrassing. Of course this metric is very ...

View answer
What are the fastest available implementations of BLAS/LAPACK or other linear algebra routines on GPU systems?
5 votes

Let me focus only on CUDA and BLAS. Speedup over an host BLAS implementation is not a good metric to assess throughput, since it depends on too many factors, although I agree that speedup is usually ...

View answer
Manipulating Matrices in Matlab
5 votes

If I got the question right, the problem is given $a_{ijk}$ and $b_{ijkl}$ to form $$ c_{ijkl} = a_{ijb_{ijlk}} $$ or in matlab notation C(i,j,k,l) = A(i,j,B(i,j,l,k)); 4 nested loops Without colon ...

View answer
Choice of basis in FEM
4 votes

In Engineering nodal bases are a good starting point for solid mechanics problems because the principle of virtual work for the discretised system $\mathbf{K} \mathbf{u} = \mathbf{f}$ reads $$ \delta \...

View answer
Parallelizing a for-loop in Python
4 votes

Before looking for a "black box" tool, that can be used to execute in parallel "generic" python functions, I would suggest to analyse how my_function() can be parallelised by hand. First, compare ...

View answer
Efficient computation of tangent of fraction of angle
Accepted answer
4 votes

I cannot think of any useful trigonometric identity that could help evaluating \begin{equation} a = \tan ( f \tan^{-1} g) \end{equation} so I would try a series expansion \begin{equation} fg + \frac{1}...

View answer
Generalized Singular Value Decomposition: only compute the r largest singular values
4 votes

I'm not an expert in this field, but being your $A$ and $B$ sparse, matlab and LAPACK are not a good choice. For sparse matrices a quick literature search confirmed that algorithms for the ...

View answer
Repeatedly solving $\mathbf{A} \mathbf{x} = \mathbf{b}$ with same $\mathbf{A}$, different $\mathbf{b}$
4 votes

Suppose $A$ is a $n\times n$ dense matrix and you have to solve $Ax_i = b_i$, $i=1\dots m$. If $m$ is big enough then there is nothing wrong in V = inv(A); ... x = V*b; Flops are $O(n^3)$ for inv(A) ...

View answer
variational formulation of linear elasticity
4 votes

Your derivation is correct up to minor imprecisions: \begin{multline} \int_\Omega \bigg((\lambda +2\mu) \nabla^2 \mathbf u \cdot \mathbf v +(\mu+\lambda)\nabla \times (\nabla \times\mathbf u)\cdot \...

View answer
What does fundamental solutions stand for in boundary element method?
Accepted answer
4 votes

In the framework of mathematical physics, the fundamental solution is the response of an infinite domain to a point source. E.g. in electrostatics the electric potential field $\varphi$ satisfies ...

View answer
FEM simulation of a material being stretched
4 votes

Modelling the stretching of objects made of rubber-like material is quite a complex task: you need to take into account hyperelastic constitutive equations, large strains, and write momentum ...

View answer
derivative of linsolve
4 votes

X = linsolve(A,B) is a convenience function of MATLAB that tries to find a sensible numerical solution to the linear system $$AX = B$$ in the general case, $A\in\mathbb{C}^{m\times n}$ with $m \...

View answer