Stefano M
• Member for 9 years, 6 months
• Last seen more than a month ago

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

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

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

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

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

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

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

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

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

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

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

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 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 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 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 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: u^h(x,t) = \mathbf{\Phi}(x)^T \, \... View answer 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 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 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 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 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 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 \...

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

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

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

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

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

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

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