149 questions
1
vote
1
answer
85
views
Iterative solvers with right-preconditioning: inconsistency between applying the preconditioner to the matrix, or to x
Summary
Using the Eigen library, I see significantly different numbers of BiCGStab iterations when I solve Ax = b with a right-preconditioner, i.e., [A M^(-1)] Mx = b, (a) by precomputing A M^(-1) ...
- 196
-1
votes
2
answers
190
views
Is there a faster method to compute these distance-based sums over a nonuniform grid of cells?
I'm attempting to implement an algorithm to compute an array of values. Unfortunately, I've found the literal implementation of that algorithm to be too slow for my use case. I suspect there's an ...
1
vote
0
answers
64
views
multivariate Fox-H numerical computation
I have some troubles in numerically calculating PDF below:
SNR of particular system model
This PDF is the expression of SNR of practical RIS-aided system: https://ieeexplore.ieee.org/document/9387559 (...
- 11
1
vote
2
answers
136
views
How many digits can I rely on?
Context
I am trying to improve my ability to recognize the number of correct decimal digits after performing short calculations involving floating-point numbers. In particular, I find it challenging ...
- 23
1
vote
1
answer
283
views
Associating B-Spline Coefficients to Control Points
From the seminal paper by Hughes I learnt that any B-Spline curve (or surface) can be described by a linear combination of B-Spline basis and control points, i.e.
C(x) = \sum_{i=1}^n N_{i,p} B_i
...
0
votes
0
answers
46
views
Vectorization of power law fitting
I am moving a code base to run using PyTorch tensor on GPU. For loops are terrible on GPU especially for small-sized data. I am trying to vectorize the function below, i.e. to have no for loops.
I ...
- 33
0
votes
3
answers
129
views
OpenMP parallelization is highly inefficient
I'm trying to numerically solve a 2-D Poisson equation with around a 200x200 grid. I'm trying to implement the diagonal method to enable parallelism:
#include <math.h>
#include <omp.h>
#...
- 131
0
votes
0
answers
160
views
NumPy np.linalg.solve give two different solution on different machines
I am studying computer science at the moment. We have the following exercise:
Execute LU decomposition given a matrix A and vector b using NumPy. We had to derive L and U by hand and came up the ...
- 1
0
votes
0
answers
166
views
How to calculate the Automatic differentiation of an arbitrary function in python
I'm using the mpmath library to compute a self defined function f(x) and needed to compute its higher order derivatives.
I found the Automatic differentiation on Wikipedia and found it useful.
# ...
0
votes
0
answers
26
views
Is it necessary to parallelize the nested function in python?
Computing higher order derivatives is a nested function and can be quite consuming. I tried to use the parallel programming to speed up the nested function.
from concurrent.futures import ...
1
vote
1
answer
404
views
How to tell how small `h` is the in the finite difference?
I'm doing a project required to compute the precise(up to a certain precession N) numerical derivatives of a function.
The usual approach was to use the finite difference types of algorithms, i.e.
(f(...
0
votes
1
answer
54
views
How to improve sub-matrix-vector multiplication c code
I have C code for generic-sub-matrix-vector (gesubmv) multiplication that computes
y(rinds) = A(rinds, cinds)*x(cinds) + y(rinds), where rinds and cinds are vectors of row and column (not necessarily ...
- 1
0
votes
1
answer
248
views
Numeric calculation of nullspace with Sympy is slow
I want to obtain the bases for kernel (nullspace) of a matrix using the Row Reduced Echelon Form (RREF).
Although scipy has the function to calculate nullspace, it does not give me sparse vectors ...
-1
votes
2
answers
715
views
How to solve a double integral with any code? [closed]
I have an double integral:
Because exp(−𝑥^2) is a non-integrable function, I have tried to solve this using the quadgk function in MATLAB, but I don't get a good result. Changing the integral's ...
- 1
0
votes
1
answer
161
views
Multiprocessing in python for numerical simulation with large objects
I am trying to put together a numerical simulation (specifically, Beta cell dynamics) based on Betram et al. 2007 (https://www.sciencedirect.com/science/article/pii/S0006349507709621). The model ...
0
votes
2
answers
152
views
Numerical Solver in Python is not able to find a solution
I broke my problem down as follows. I am not able to solve the following equation with Python 3.9 in a meaningful way, instead it always stops with the initial_guess for small lambda_ < 1. Is there ...
0
votes
1
answer
397
views
Sympy - speed up expressions evaluation
I compute something in python with sympy but this require too much time.
The code looks like this:
def compute_aux_lyupanov_val_k(y, A, B, C, D, k):
res = []
res.append(y)
for i in range(k-1):
...
- 115
20
votes
1
answer
940
views
why does numpy matrix multiply computation time increase by an order of magnitude at 100x100?
When computing A @ a where A is a random N by N matrix and a is a vector with N random elements using numpy the computation time jumps by an order of magnitude at N=100. Is there any particular reason ...
- 450
-2
votes
1
answer
68
views
Can I obtain different results from iterative and recurive functions?
My code is supposed to calculate the 100th element of the sequence $x_0=1 ; x_i=\dfrac{x_{i-1}+1}{x_{i-1}+2}, i=1,2, \ldots$
I wrote iterative and recursive functions, but the results are not equal. ...
- 125
-2
votes
2
answers
606
views
Solving a non linear equation using the secant method in Fortran
Sooner than expected, I have yet another problem with my Fortran code.
This time I'm trying to find the zeros of a function, which I can normally do, but now I have a function that has more than one ...
- 149
0
votes
1
answer
672
views
Create a web interface for a Python numerical simulation code (and where to host it)
I need to create a web interface for a Python numerical simulation code, and host both (interface and python code).
When a user will run a simulation, it could take hours, or even days for the ...
user14028806
-1
votes
1
answer
545
views
Wrapping around very large angles to [0,2*PI] adds too big numerical error [closed]
In below sample, I try different methods to wrap an angle between 0 and 2 x PI. While doing this to compute cosine of same angles, both methods cause a different amount of numerical error compared to ...
0
votes
1
answer
181
views
Implementation of piecewised function in python with numerical ODE solution. TypeError: float() argument must be a string or a number, not 'OdeResult'
I have a reaction system whose mechanism is defined in the form of a nonlinear first-order ODE as:
dx(t)/dt = - gaf(1 - rho)(1 - exp(-kt))(x(t)^2 + bx(t))/(cx(t) + p)
where, g, a, f, b, c, rho, k and ...
- 1
0
votes
0
answers
139
views
taking derivative using sparse matrix
I’m trying to solve the following problem for around 2 days, but have not had success.
I want to calculate the derivative using the sparse matrix, but the result doesn’t true. I think there is a ...
- 1
3
votes
1
answer
980
views
Overloading of addition for arrays in Rust
Trying to sum two arrays ([1.0, 2.0] + [3.0, 4.0]) fails in Rust as addition isn't defined for arrays. I tried to overload the operation:
use std::ops::Add;
type float = f64;
impl<const N: usize>...
- 321
0
votes
1
answer
44
views
Create B such that B @ X.ravel() == (A.T @ X + X @ A).ravel()
I'm not new to numpy, but the solution to this problem has been eluding me for a while now, so I thought that I would ask here.
Given the equation Y = A.T @ X + X @ A
where A and X are square matrices....
user9413641
1
vote
1
answer
571
views
Should one calculate QR decomposition before Least Squares to speed up the process?
I am reading the book "Introduction to linear algebra" by Gilbert Strang. The section is called "Orthonormal Bases and Gram-Schmidt". The author several times emphasised the fact ...
- 313
0
votes
3
answers
257
views
Need to find when value is == 0, but I cannot due to numerical errors
I have 2 large lists of vectors (>10,000 vectors each, say vi and wi) and I am trying to find when vi cross-product wi = 0, or, when vi x wi = 0.
The lists of vectors are previously calculated (...
0
votes
1
answer
1k
views
How to remove spikes in solution and produce smooth interpolation with scipy?
I'm writing a numerical solution to a partial integro-differential equation, and I need it to run quickly so I found that scipy.integrate.simps is best, but it's not always 100% accurate and produces ...
1
vote
2
answers
96
views
How to compute Brjuno function in Maxima CAS?
I try to compute real 1-Brjuno function for real numbers in x in [0, 1] range
G(x):= block(
[g],
if (x=0)
then g : 0
else g : float(1/x - floor(1/x)),
...
- 1,406
0
votes
1
answer
168
views
Finding f(g) from f(x) and g(x) numerically
I'm studying an article where there are two functions:
$e(x) = (e_0/8)\cdot[(2x^3 + x)\sqrt{(1 + x^2)} − \sinh^{−1}(x)]$
$p(x) = (e_0/24)\cdot[(2x^3 - 3x)\sqrt{(1 + x^2)} + 3\sinh^{−1}(x)]$
The ...
- 23
0
votes
1
answer
265
views
C floating point exception handling
The manual pages for fenv.h (feraiseexcept and its ilk) are unusually uninformative; the examples I find online (cppreference.com and others) are very minimal.
How are they supposed to be used? In ...
- 11.9k
0
votes
1
answer
98
views
numerical prediction of a pendulums motion in python not predicting properly
I'm working on a numerical prediction of a pendulums motion for a project. while running the program I've noticed the time between each peak (when the value reaches the same theta value as the ...
1
vote
1
answer
935
views
fortran : invalid memory reference
I encountered this error while solving the diffusion equation on fortran using alternating directional implicit(ADI) method.
Program received signal SIGSEGV: Segmentation fault - invalid memory ...
- 51
1
vote
1
answer
300
views
Filling pixels under or above some function
Seems like a simple problem, but I just cant wrap my head around it.
I have a config file in which I declare a few functions. It looks like this:
"bandDefinitions" : [
{
"0&...
- 337
2
votes
2
answers
2k
views
Velocity verlet algorithm - energy increasing for n-body problem
The problem
I implemented velocity verlet algorithm to compute trajectories of 2 bodies interacting with each other gravitationally (newtonian gravity only). Orbiting smaller body has a very small ...
- 89
0
votes
1
answer
595
views
Numerical solution for a pendulum
I am having a trouble in providing a graphical representation of my system, which happens to be a harmonically driven pendulum. The problem is shown below for reference.
Problem
The source code I used ...
3
votes
1
answer
168
views
More Robust Integrator cpp
I'm trying to produce a relativistic voigt distribution, the convolution of a relativistic Breit-Wigner distribution and a gaussian function.
MWE:
double relativisticBreitwigner_pdf(double energy, ...
0
votes
1
answer
190
views
Functions intersection approximation
Let f1 and f2 be two functions in the range of [a, b], and maxerr the required approximation. They both differentiable and continuous in this range. I should return an iterable of approximate ...
- 69
0
votes
0
answers
871
views
Alternatives to scipy.integrate
I'm new to numerical integration with Python and I was wandering what is the best Python package for this purpose. Of course I discovered scipy.integrate that is pretty nice, but I wanted to know if ...
- 25
0
votes
0
answers
116
views
Is there a way to manipulate function arguments by their position number?
I wish to be able to manipulate function arguments by the order in which they are given. So,
void sum(int a, int b, int c)
std::cout<< arguments[0]+arguments[1];
sum(1,1,4);
should print ...
- 101
0
votes
2
answers
95
views
C function call Numerical Recipes Chapter 10.3
I have this function. I'm kinda green in C and I don't know how to call this function to get a result:
float dbrent(float ax, float bx, float cx, float (*f)(float),float (*df)(float), float tol, float*...
- 11
0
votes
1
answer
92
views
CUDA access matrix stored in RAM and possibility of being implemented
Recently I started working with numerical computation and solving mathematical problems numerically, programing in C++ with OpenMP. But now my problem is to big and take days to solve even ...
- 15
-1
votes
1
answer
310
views
issue with converting rate equation to python code
I am trying to convert these rate equations to python code, I have made I lot of research but can't seem to get any clear path to follow to achieve this, please any help will be appreciated
This is a ...
3
votes
2
answers
108
views
How to improve the runtime of this matrix generation loop in python3?
I have code that is simulating interactions between lots of particles. Using profiling, I've worked out that the function that's causing the most slowdown is a loop which iterates over all my ...
- 301
2
votes
5
answers
939
views
Accurate sqrt(1 + (x/2)^2) + x/2
I need to compute sqrt(1 + (x/2)^2) + x/2 numerically, for positive x. Using this expression directly fails for very large values of x. How can I rewrite it to obtain a more accurate evaluation?
- 21.2k
0
votes
0
answers
30
views
Quark in producing fitted values using LM model in R
A colleague and I noticed this interesting quark in the lm function in R.
Say, I am regressing y variable on an x variable and x variable is a factor level variable with two categories (0/1).
When I ...
- 1,328
1
vote
1
answer
489
views
How to ensure my optimization algorithm has found the solution?
I am performing a numerical optimization where I try to find the parameters of a statistical model that best match certain moments of the data. I have 6 parameters in total I need to find. I have ...
- 137
0
votes
1
answer
340
views
store values of function to prevent from running again
Say I have some complicated function f(fvar1, ..., fvarN) such as:
def f(fvar1,..., fvarN):
return (complicated function of fvar1, ..., fvarN).
Now function g(gvar1, ..., gvarM) has an expression ...
- 123
1
vote
1
answer
486
views
Numerical Integration - Summations and Functions in Python
Forgive my ignorance, but the function H is actually an infinite sum (including -n terms). I'm hoping to truncate at larger values on the order of the 100s ideally. My code seems to work but I am not ...
