Questions tagged [numerics]

The tag has no usage guidance.

Filter by
Sorted by
Tagged with
2 votes
2 answers
2k views

Why does scipy's odeint function give a non-monotonic solution for a problem whose solution should be monotone?

The solution to the ode below looks like it is monotonically increasing: However on closer inspection we see that it is not: How can I ensure that the numerical solution is monotonically increasing?...
  • 873
1 vote
1 answer
348 views

How to deal with indeterminate function limit?

How do I ensure that my function below is well conditioned as $s$ approaches $\infty$? The problem I get is that for large $s$ the function returns an indeterminate form $\frac{0}{0}$. I would ...
  • 873
7 votes
1 answer
5k views

What is numerical damping in the context of time-dependent FEM solvers?

Comsol Multiphysics (a popular FEM package) includes two time-stepping algorithms (IDA aka BDF, and Generalized-alpha), described in their documentation as follows (quoted here under Fair Use; ...
6 votes
2 answers
593 views

Evaluating 6D Gaussian Integral in Matlab

I have to compute the accuracy of a new Gaussian mixture fitting algorithm. One of the tests include computing the probabilities in certain intervals in a 6D hyperspace. Also, the integral of the ...
3 votes
1 answer
2k views

Partial differential equations with octave [closed]

I need to find a numerical solution for $-\Delta U = f$, on the $\Omega = [0,1]^2$, with $ U|_{\partial \Omega} = 0$. I found a method: POISSONFD in ...
3 votes
1 answer
328 views

Books and references on implementing finite difference codes for PDEs

Are there any good books or references on implementing finite difference methods for PDEs? Specifically, I'm looking for something comparable to Gockenbach's book Understanding and Implementing the ...
  • 707
7 votes
2 answers
2k views

Choosing epsilons

Most numerical algorithms require an epsilon to be chosen in order to be robust and provide meaningful results. Choosing machine epsilon is usually too aggressive. Barring any special knowledge ...
3 votes
1 answer
647 views

Solving a system of nonlinear equations with an ODE solver is faster than with the Newton method?

This is somehow unexpected, but my recent experience with solving a system of nonlinear equations is that treating them as the right hand side of a system of ordinary equations and then evolve the ...
  • 133
6 votes
2 answers
240 views

What are the numerical methods for huge polynomial systems?

Let a system of $n$ polynomial equations of degree $d$ with $m$ variables. I'm interested in a sparse system with $d = 3$, $n \sim 2000000$, $m \sim 50000$ and integer coefficients. What techniques ...
1 vote
1 answer
73 views

eigenvalue of small symmetric matrices

If I am to solve a symmetric eigenvalue system $A=QDQ^T$, where $A\in\mathcal{R}^{n\times n}$ and $n$ is small (in the range 4 - 64); I want all the eigenvectors and eigenvalues; There are two major ...
4 votes
1 answer
150 views

State-of-the-art for active set optimization algorithms?

Given a problem like this: $$ \text{min } ||Ex-f|| \text{ s.t.}$$ $$ Gx \ge 0$$ $$ Cx = d $$ And assuming that the matrices are medium sized (dimensions in the low thousands) and dense, what's the ...
1 vote
1 answer
41 views

replace non-smooth discrete values with analytical function

I do have a Diffusion coefficient in a convection diffusion PDE which is discontinuous and looks like (concentration on the x-Axis): For numerical reasons i use the integrated form: I calculate the ...
  • 333
10 votes
1 answer
3k views

Fortran 90/95: Deallocating variables

I understand the crucial importance of freeing memory when certain variables or arrays need to be reused later in the program, or may not be in use for a while. However, in my experience with ...
  • 165
2 votes
1 answer
254 views

Energy of damped harmonic oscillator begins increasing with very large Q

I have numerically integrated the (reduced) homogeneous equation of a damped harmonic oscillator in order to see how the error propagates. $$\frac{d^2 X}{d\phi^2} + \frac{1}{Q}\frac{dX}{d\phi}+X(\phi)...
  • 263
1 vote
1 answer
141 views

Importance of change the Time Step Value for the Accuracy of a Transient CFD Simulation

I have a transient simulation for a case from 0 to 3 sec , actually i interest the solutions on time range from 2 to 3 sec , as my velocity and mesh size , my time step should be 0.0000000001 sec , my ...
0 votes
3 answers
2k views

Looking for ways to speed up the numeric evaluation of a symbolic expression in Matlab

{Summary: I have a symbolic expression DCritnF expressed in terms of two variables x1 and x2. I need to find it's numeric value and I used combination of double and subs as given below. ...
  • 35
9 votes
1 answer
527 views

Algorithm to calculate the exponential of an Hessenberg matrix

I am interested in computing the solution of a lage system of ODEs using a krylov method as in [1]. Such method involve functions related to the exponential (the so-called $\varphi$-functions). It ...
9 votes
2 answers
2k views

Help deciding between cubic and quadratic interpolation in line search

I'm performing a line search as part of a quasi-Newton BFGS algorithm. In one step of the line search I use a cubic interpolation to move closer to the local minimizer. Let $f : R \rightarrow R, f \...
  • 193
0 votes
2 answers
521 views

Numerical evaluation of the first and second complete elliptic integrals

To get a numerical evaluation of the first (K) and second (E) complete elliptic integrals: $$K(k)=\int_0^1\frac{dt}{(1-t^2)^{1/2}(1-k^2t^2)^{1/2}}, \ \ \ \ \ E(k)=\int_0^1\frac{(1-k^2t^2)^{1/2}}{(1-t^...
  • 103
3 votes
1 answer
800 views

Accurate implementation of the logarithm of the incomplete Beta function in C++?

I need an accurate implementation (for use in C/C++) of the logarithm of the incomplete Beta function: $$\log \mathrm{B}(x,y;\alpha,\beta) = \log \int_x^y t^{\alpha-1}(1-t)^{\beta-1}\mathrm{d}t$$ ...
  • 1,629
4 votes
3 answers
6k views

How to get ODE solution at specified time points?

The code below basically illustrates my problem. It is a test code for a pendulum. I solve it using a method suggested on https://stackoverflow.com/questions/12926393/using-adaptive-step-sizes-with-...
  • 873
2 votes
0 answers
199 views

boundary conditions with non-constant coefficients in cell centred finite volume method

Suppose am solving the heat conduction equation in 1d with Dirichlet boundary conditions. The thermal conductivity $k$ is a non constant function. So $-(k(x)u'(x))' = f(x)$ The value of $k$ enters the ...
  • 435
3 votes
1 answer
315 views

Decompositions of sparse symmetric matrices and methods for solving large linear equations

I asked the same question on mo.se and it was suggested that scicomp would be a better forum for it. So here it is: I am writing code for solving linear equations of the form $$A_{n\times n}\cdot x=...
  • 417
5 votes
1 answer
3k 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_{...
  • 603
3 votes
1 answer
155 views

Imposing symmetry plane boundary condtition

I want to impose symmetry plane boundary condition for a solid mechanics problem. I googled around and found out that in many places people say to "forbid displacemnts out of symmentry plane and ...
  • 325
8 votes
3 answers
16k views

plotting discontinuous functions

I need help plotting the Heaviside function: Real analysis often involves constructing bizarre functions which are intuitively correct, but ultimately wrong. See the great book Counterexamples in ...
1 vote
0 answers
571 views

Convolution of two radially symmetric functions for double logarithmic plot

Question: I have problems with sidelobes when computing the convolution of two radially symmetric functions. I wonder whether I should just switch from second order to fourth order and refine the grid ...
5 votes
1 answer
156 views

Scheme to alleviate (numerical?) instability in system of coupled nonlinear ODEs

I'm solving a system of nonlinear ODEs that take the form $Q_{nm} \ddot{y}_m + S_{nkl}\dot{y}_k\dot{y}_l +V_n = 0$ where Einstein summation is assumed, $y_i$ are the dependent (complex) variables, ...
  • 385
1 vote
0 answers
43 views

Error analysis and the Model Problem [closed]

In numerical methods for ODE's, the model problem y' = cy where c is complex is regarded as sufficient in performing error analysis for different methods in ...
  • 121
2 votes
1 answer
438 views

LSA, SVD and the Frobenius norm

In Latent Semantic Analysis one uses the SVD to perform a dimensional reduction of the term-document matrix, via the Eckart-Young theorem. Now, the rank $k$ approximation obtained by E-Y is proven to ...
  • 141
2 votes
1 answer
144 views

Lanczos algorithm with thick restart on a dynamic matrix

According to a recommendation, this is a re-post of that. currently, I'm working on a way to compute the 2 biggest eigenvalues of a real, symmetric, huge and sparse matrix that changes a few entries ...
  • 21
4 votes
1 answer
191 views

Numerical spherical integration

In a high-dimensional setting, say $d \gg 5$, what is a recommended way of evaluating a spherical integral of a smooth (non-symmetric) function $f(\mathbf{x})$? $ \int_\mathcal{S_r} f(\mathbf{x}) \...
  • 870
1 vote
1 answer
570 views

RATTLE numerical integrator example

I want to understand how the RATTLE algorithm works. Can somebody give me an example (in pseudocode or using any programming language like python or matlab) of how would I implement a numerical ...
8 votes
1 answer
197 views

Radial integration of expensive function with Bessel weights

I need to calculate the integral $$I = \int_0^R f(r)J_n\left(\frac{z_{nm}r}{R}\right)rdr$$ where $J_n$ is the $n^{\mathrm{th}}$ order Bessel functions of the first kind, $z_{nm}$ is its $m^{\mathrm{...
  • 81
3 votes
1 answer
956 views

Von Neumann stability analysis in 3d

I need to get a stability criterion for the numerical scheme for equation $$\frac{\partial u}{\partial t}-\frac{\partial^2 u}{\partial x^2}-\frac{\partial^2 u}{\partial y^2}-\frac{\partial^2 u}{\...
  • 501
3 votes
1 answer
564 views

Condition number of an algorithm

I am stuck with a problem about finding the condition number of an algorithm. I tried to find an example, but i couldn't. Can anybody help me, please? Given is $f(x)=\ln(x)$. We have the algorithm $...
3 votes
1 answer
1k views

Courant Friedrichs Lewy condition - how to get it?

I am interested, how can we get CFL condition for every type of PDE? It's known that for 1st order linear equation $$\frac{\partial u}{\partial t}+a\frac{\partial u}{\partial x}=0 $$ CFL is get from ...
  • 501
0 votes
0 answers
135 views

Algorithms to compute largest gap between smallest nonzero eigenvalues of sparse symmetric matrix

I am looking mainly for c/c++ implementations but also for theoretical algorithms to compute gaps between smallest positive eigenvalues of symmetric, singular matrix or real numbers. To be precise, I ...
1 vote
0 answers
275 views

Numerical integral with a weakly singular kernel with a satisfactory precision

I am working on a numerical method for time fractional PDE. One problem is that I must compute a numerical integral of the following form: $$ \begin{equation} \int_0^{t_m} (t_m-s)^{-\beta}f(s)ds \end{...
  • 121
3 votes
2 answers
909 views

C# implementation of the gamma function that produces correct answers at positive integer inputs?

I need a C# implementation of the gamma function that produces correct exact answers at positive integer inputs. I took a look at MathNet.Numerics Meta.Numerics. In both cases, if you calculate ...
2 votes
1 answer
967 views

Finite difference discretization on a circle

I am trying to discretize the differential operator $\frac{d^2}{dx^2}$ acting on $S^1 = [0,1]$ using finitely many points around a circle at $0, \frac{1}{N}, \frac{2}{N}, \dots, \frac{N-1}{N}$. Here ...
4 votes
2 answers
303 views

Solving system of differential equations with interconnected boundary conditions

I am trying to solve the following system of differential equations numerically over the domain $x=0$ to $x=D$. The main difficulty is that the boundary conditions are interconnected and depend on the ...
  • 5,299
5 votes
1 answer
2k views

Numerical gradient in spherical coordinates

Assume that we have a function $u$ defined in a ball in a discrete way: we know only the values of $u$ in the nodes $(i,j,k)$ of spherical grid, where $i$ is a radius coordinate, $j$ is a coordinate ...
  • 501
1 vote
1 answer
1k views

"boundary" vs "interface"?

I am working with biofilm and there are many documents talking about boudary conditions while others talks about interface or both of boundary and interface. So, boundary and interface are the same (...
2 votes
1 answer
169 views

Extended finite element method vs $P_k$-bubble element

Can you show me the main differences between 2 methods? I find out 2 reasons but I don't know they are right or not. XFEM is constructed base on enrichment functions whereas P1-bubble is constructed ...
9 votes
3 answers
3k views

Finite Element Method vs Extended Finite Element Method (FEM vs XFEM)

What are main differences between FEM and XFEM? When should we (not) use XFEM intead of FEM and vice versa? In other words, when I meet a new problem, how I can know to use which one of them?
4 votes
1 answer
371 views

Methods for solving BVP for DAE

I look for a numerical method to solve boundary value problems for systems of differential and algebraic equations of the form F(x,y,y') = 0, G(x,y) = 0, y(a) = ya, y(b) = yb, where y = (y1, y2, ... ...
  • 43
7 votes
2 answers
175 views

For a non-linear PDEs should the source term be discretised at $u_j$ or averaged over $(u_{j+1} + u_{j-1})/2$?

The non-linear Poisson equation in one-dimension, $$ 0 = \frac{\partial^2u}{\partial x^2} - f(u) $$ can be discretised as to give, $$ u_{j-1} -2u_{j} + u_{j+1} = h^2 f(u_j) $$ where $h$ is the ...
  • 5,299
2 votes
1 answer
98 views

Implicit Finite difference scheme for a PDE with only one boundary

I am looking at a few reaction-diffusion equations of the form $\frac{dP}{dt} = D\left(\frac{d^2P}{dr^2} + \frac{2}{r}\frac{dP}{dr}\right) - a(P)$ I know the initial conditions and the boundary ...
  • 121