All Questions

Filter by
Sorted by
Tagged with
6
votes
2answers
381 views

Implementing the pressure correction method using finite elements

Ok so I am nearing the completion of my finite element Navier-Stokes solver that uses the $\theta$-method for time stepping and the pressure correction method for the pressure. I am following the ...
6
votes
2answers
1k views

matlab eigs: wrong eigenvalues for tridiagonal matrix

I try to compute eigenvalues of the tridiagonal matrix coming from finite difference scheme. For small mesh size, eigs works well. But for large size it fails. Here is an example where eigs fails. Is ...
6
votes
4answers
999 views

Need Fortran 77 compiler

Does anyone know a compiler for Fortran 77 available as a free download? I have pre-written 77 code from a source published in the early 90's that I just need to compile, build, and run. But I don't ...
6
votes
1answer
360 views

Caveats of Hessian free method

Hessian free iterative optimization techniques like Newton-CG, do not explicitly compute the Hessian but instead approximate the product of the Hessian with a vector through finite difference. The ...
6
votes
1answer
169 views

A simple PDE solution question

I need to ask a question about partial derivatives. I want to solve this equation (steady state, one dimensional continuity equation): $$\frac{\partial (\rho u)}{\partial z}=0$$ which is equivalent to:...
6
votes
1answer
167 views

Can heat distribution in an optical element irradiated by laser be oscillating?

I am modelling a heat distribution in optical element irradiated by laser. System is radially symmetric, and element is thin, i.e. heat value depends only on distance from center. Heat is received via ...
6
votes
2answers
835 views

Python package for (adaptive) function plotting

Are there any mature Python packages that can plot functions, and possibly use adaptive sampling? I am looking to pass only a function (can be a numerical black box) and a range, and expect a plot as ...
6
votes
2answers
546 views

Applying matrix square root inverse in matrix-free regime

Let $A$ be a large symmetric positive definite matrix, and suppose that we can efficiently apply $A$ and have a fast solver to apply $A^{-1}$, but we do not have access to the matrix entries for ...
6
votes
0answers
260 views

Integrating highly oscillatory functions

I have a logarithmic grid, upon which i have two functions that are similar to this one (this is only the last 100 points): These are essentially very similar to a Sin function at this point. I need ...
6
votes
1answer
205 views

Does mean removal increase accuracy of numerical differentiation?

I wish to compute the derivative of a vector through numerical differentiation. Let's say, we use a standard 2nd order central difference scheme, to arrive at a differentiation matrix, and apply it on ...
6
votes
1answer
857 views

Algorithm for solving system of quadratic equations and linear equations

Let $x \in R^N$. From a Spectral Chebyshev collocation method, I have a system of quadratic and linear equations. Denote them, $$ x^T Q_i x + L_i^T x = 0 $$ and $$ A x = 0 $$ Furthermore, I know ...
5
votes
1answer
254 views

Convergence of Jacobi's method for a semilinear elliptic PDE

I have an iterative finite difference scheme for the Poisson equation $$ \nabla^2 u=-\rho $$ It's the Jacobi method, which has the form (for 1D systems) $$ u^{n+1}_{i} = \frac{1}{2}(u^n_{i+1} + u^n_{...
5
votes
3answers
794 views

Seemingly non-unique Cholesky factor via QR rectangularisation

I am trying to implement an algorithm from a paper which makes use a QR factorization of a real matrix $A$ as a means of one of forming the Cholesky factor of $A^T A$ without explicitly forming $A^T A$...
5
votes
2answers
196 views

Is there a special algorithm for computing the convex hull ordering when the candidate points are on the hull?

I'm dealing with a set of points which are already placed on the 2D hull boundary: a convex polygon. I know this for sure. However, the point set is not ordered, and I need the polygon points to be ...
5
votes
1answer
187 views

What is the correct way of performing numerical experiments on desktops?

Suppose I want to set up an experiment to measure the performance of some numerical code on a desktop machine running Linux/Windows/MacOS. What kind of environment should I arrange in order to get ...
5
votes
0answers
85 views

Interface Formulation at Finite Volume Boundaries when using the Dual Mesh

When using the dual mesh (vertex-centered) for finite volume methods, you end up with a cell center at the boundaries between materials. It is possible that the equations being solved in each ...
5
votes
1answer
1k views

FEM: which is the correct way to impose Dirichlet B.C

I know the Neumann B.C. is implicit in FEM language. However, I have seen at least two ways to impose Dirichlet B.C. e.g. for the following problem 1D, $$\nabla^2 u + \nabla u= 0, u_{left}= 1, u_{...
5
votes
3answers
6k views

Algorithm for Principal Eigenvector of a Real Symmetric 3x3 Matrix

I have a 3x3 covariance matrix (so, real, symmetric, dense, 3x3), I would like it's principal eigenvector, and speed is a concern. Is there a fast algorithm for this specific problem? I've seen ...
5
votes
2answers
3k views

Fast algorithms to find the eigenvalues of some matrix on intervals of interest

I am curious how to quickly compute the eigenvalues for arbitrary matrices, sparse or dense, restricted on some given interval of interest. Suppose we have an arbitrary $n\times n$ matrix $A$, ...
5
votes
1answer
2k views

Stability of numerical method for 1D Burger's equation

I am trying to solve 1D viscous Burger's equation numerically and I cannot apply von Neumann analysis because the equation is non-linear. How do I predict the stability criteria for my system? I also ...
5
votes
3answers
351 views

PRIMA gives an unstable result?

I am working with Modified Nodal Admittance representation of circuits. I am doing Model Order Reduction using PRIMA on MATLAB. I am considering these circuits as Descriptor State-Space systems. I ...
5
votes
3answers
215 views

Is there a way to reduce aberration in computations of planets' trajectories?

I don't think the title is very accurate , sorry for that. I simulate bodies in space using two timestep: the TIMESTEP is the Δt wich I use to make the calculation and XTIME is the number of times ...
5
votes
0answers
170 views

Symmetric sparse direct solvers in scipy

scipy.linalg.solve, in its newer versions, has a parameter assume_a that can be used to specify that the matrix $A$ is symmetric ...
5
votes
3answers
178 views

Algorithms for radiation treatment planning

I have a medical physics problem - I want to maximise the dose absorbed by a brain tumour whilst minimising the dose in the rest of the brain, especially certain organs, such as the pituitary gland, ...
5
votes
2answers
924 views

minimization problem: sum of Rayleigh quotients

I would like to find $x$ which minimizes the following equation: $\frac{x^HAx}{x^HBx} + \frac{x^HCx}{x^HDx}$ where A, B, C, D are positive-definite. $x$ is not a very large vector (<1000 elements ...
5
votes
1answer
91 views

Is the bandwidth of indefinite A equal to its factor L in LDL^T?

In George, Liu, and Ng's book Computer Solutions of Sparse Linear Systems, it has been shown that bandwidth of $A$ is equal to bandwidth of its factors in $LL^T$.(section 4.3) However, I guess this is ...
5
votes
1answer
568 views

How is rigid bodies implemented in finite element codes

I am writing a finite element code for structural analysis, and I want to implement rigid bodies. How is this usually done? Say that I have a square mesh, with one half of the mesh being defined rigid ...
5
votes
2answers
928 views

Matlab: Solving $u_t = f(u) u_{xx} + g(u) u_x + h(u) u$

I am trying to solve the following equation numerically: $$\dfrac{\partial u}{\partial t}=\dfrac 3 x\dfrac \partial{\partial x}\left(\sqrt x \frac \partial {\partial x}(\nu(u,u^2,\cdots,u^k)u\sqrt x)\...
4
votes
4answers
338 views

Role of weight function in Galerkin methods

I have difficulties in understanding the role of the weight function $w(x)$ that occurs in the solution of PDEs via the Galerkin approach. Consider a linear differential equation of the form $$ \...
4
votes
1answer
60 views

How to deal with pseudo-compressibility of lattice Boltzmann method when you are calculating mass flux?

In lattice Boltzmann method, we have a concept, which is called pseudo-compressibility and it is defined based on the fact that LBM simulates incompressible flows by having small Mach number to ensure ...
4
votes
3answers
425 views

How to implement this trigonometric polynomial maximum finding semidefinite program

Hi All, I posted this on the math.se site, but this may be a better location. I need a method of finding the maximum of a real valued trigonometric polynomial where I can trade accuracy for speed. ...
4
votes
1answer
339 views

Why my linear congruential generator generate low quality random numbers?

I am implementing an arbitrary bits random number generator (in [$0,2^n$)) and also want to ensure that the generator always generate unique numbers until all possible numbers generated, so I ...
4
votes
2answers
251 views

How “sparse” should a sparse matrix be to see benefits?

I have a matrix, whose size scales as $2^N$ (assume even $N$). In each row of the matrix, only about $2^{N/2}$ of the entries are filled ($N$ can be somewhere between 10 and 40, depending on what's ...
4
votes
1answer
1k views

How to find conditional Lyapunov exponents

For a synchronized ODE system,I wanted to know if there is a program code in MATLAB available for plotting the conditional Lyapunov exponent. For example, synchronization of identical Rossler system ...
4
votes
2answers
142 views

Maximization variant of semidefinite programming (SDP)

Consider the following program: $$\max_{\pmb a} \sum_i z_i\\ u.c. \pmb a \pmb P_i\pmb a^\top\geq z_i$$ where $\pmb a \in\mathbb{R}^p$ and the $\pmb P_i$ are all symmetric positive semidefinite ...
4
votes
4answers
146 views

Can the mesh generation methods in FVM and FEM be totally based on the knowledge of the mesh generation theory in computer graphics?

The main references of mesh generation methods in computer graphics (CG) I found are discrete Differential Geometry [1] and a famous book "Polygon Mesh Processing" [2], while the "...
4
votes
1answer
351 views

Efficient and stable computation of inverse CDF

What is the most efficient and numerically stable algorithm for computing the inverse CDF $F^{-1}(y)$ of a probability function, assuming that both the PDF $f(x)$ and the CDF $F(x)$ are known ...
4
votes
2answers
3k views

Testing 1D Poisson Solver

I'm trying to test a simple 1D Poisson solver to show that a finite difference method converges with $\mathcal{O}(h^2)$ and that using a deferred correction for the input function yields a convergence ...
4
votes
0answers
354 views

Iteratively finding both left and right eigenvectors for non-symmetric complex matrix

I have a complex, non-Hermitian matrix $\mathbf{A}$, for which I need to find a few eigenvalues and eigenvectors in the generalised eigenvalue problem: $$\mathbf{A}\cdot \mathbf{x} = \lambda \mathbf{...
4
votes
0answers
78 views

Large residual when integrating 2nd order ode close to singularity with SciPy ode / ODEPACK

I am trying to integrate a 2nd order ODE with a singularity at close to the initial condition. Why do I get large residuals when I plug-in the result of my integration back into the ODE? The equation ...
4
votes
2answers
369 views

LAPACK - singular matrices - what does the positive integer info mean?

please can you help me with my code - I use Lapack to solve complex matrix (quite biq) and do it in two steps: I call zgetrf (LU factorization) and then ...
4
votes
1answer
251 views

Multigrid stops converging when more grid levels are used

I'm having a problem with multigrid code I wrote. If I solve Laplace's equation in 2D and use more than 5 grid levels, the V-cycles stop converging after a few cycles (see below, convergence factor > ...
4
votes
2answers
570 views

Maxwell vs Kepler for GPU computation

I am looking at the Nvidia GT-860M which comes in both the old Kepler and new Maxwell architecture. The old one (1152) seems to have almost twice the cores as the new one (640). The new one has a ...
4
votes
1answer
1k views

Using Gmsh to create a mesh with zero thickness (quad) interface elements

I acknowledge the following post where a similar question is posed and a very nice answer has been provided: Is there a mesh generator that will generate zero thickness elements for interfaces? ...
4
votes
0answers
474 views

full rank update to cholesky decomposition for multivariate normal distribution

This question is a specialization of full rank update to cholesky decomposition, to which I hope to get a more positive answer. When calculating the minus log of the multivariate normal distribution, ...
4
votes
1answer
130 views

Constructing the origin position by transforming distance information

Suppose a set of $n$ points, $n\in M$, is given in some $d-$dimensional space, $X\in\mathbb{R}^{n\times d}$. Among these $n$ points, some $k\in K$ are selected, so $k<n$, and the distances from ...
4
votes
3answers
2k views

4th-order Runge-Kutta method for coupled harmonic oscillator

I’m attempting to write a C program to gather values from a coupled spring system: There is a wall, connected to a mass $m_1$ by a spring, then this mass is connected to a second mass $m_2$ by another ...
4
votes
1answer
367 views

Best platform for complex SDPs with n and m around 5-15K?

I am looking to solve a class of SDPs with complex entries, with the semi-definite cone $S^n$, $n$ around 5000 to 15000. Also, $m$, the number of equality/inequality constraints is close to $n$. I ...
4
votes
1answer
160 views

Factorization for reweighted least squares

I am solving a problem using an iteratively-reweighted least squares method: http://en.wikipedia.org/wiki/Iteratively_reweighted_least_squares Essentially this requires solving a number of least-...
4
votes
2answers
185 views

Intervals where the sign of a polynomial can be computed reliably

This is a follow-up of a previous question. Let $p$ be a polynomial with floating-point coefficients. Is there a method for finding intervals where evaluating $p$ in floating-point arithmetic always ...

15 30 50 per page
1
11 12
13
14 15
17