Skip to main content
Share Your Experience: Take the 2024 Developer Survey
25 votes

More stable algorithm to calculate `sqrt(a^2 + b^2) - abs(a)` in MatLab

A manipulation that may help is the following. Assume for simplicity $a>0$. We have the identity $$b^2 = (a^2+b^2)-a^2 = (\sqrt{a^2+b^2}-a)(\sqrt{a^2+b^2}+a),$$ hence $$ \sqrt{a^2+b^2}-a = \frac{b^...
Federico Poloni's user avatar
18 votes
Accepted

More stable algorithm to calculate `sqrt(a^2 + b^2) - abs(a)` in MatLab

You can break down the domain of your function into three distinct cases: $|a|\gg |b|$: In this case, $\sqrt{a^2+b^2} \approx |a|$ and a naive application of the formula will likely result in poor ...
Wolfgang Bangerth's user avatar
7 votes

'eigs()' in Matlab gives inaccurate eigenvector when only several eigenvalues are calculated

Eigenvectors are not unique, especially when you have repeat eigenvalues. Let's define each unique eigenvalue as $\lambda_1$, $\lambda_2$, ..., $\lambda_i$, where $\lambda_i$ has a multiplicity of $...
helloworld922's user avatar
5 votes

enough conditions to check that a matrix doesn't have Cholesky factorization while factorizing it

You are trying to figure out what is known as Modified Cholesky. There are a lot of resources on it, but McSweeney's thesis is a good starting point "Modified Cholesky Decomposition and ...
Abdullah Ali Sivas's user avatar
5 votes

More stable algorithm to calculate `sqrt(a^2 + b^2) - abs(a)` in MatLab

This is basically a variant of Federico Poloni's excellent answer, based on the principle that one might want to try proven building blocks first. I do not know if Matlab provides this functionality, ...
njuffa's user avatar
  • 1,875
4 votes
Accepted

Generalized eigenvalue problem for large, potentially ill-conditioned systems

For large systems, any direct solver methods tend to be a dead end as what often starts as a sparse system ends up becoming dense. In fact, just storing all eigenvectors is itself typically impossibly ...
Mikael Öhman's user avatar
4 votes

Matlab eigs function with function handle

This is for algorithms like Lanzcos-Arnoldi or Krylov-Schur, as mentioned in the reference part of the documentation page. Generally these are Krylov-space generalizations of the power iteration, ...
Lutz Lehmann's user avatar
  • 6,109
4 votes

Reordering eigenvalues in Schur factorization - MATLAB ordschur and LAPACK dtrsen not producing the same results

You should use COMPQ = 'V'; instead of COMPQ = 'N'; accroding to the official documentation of LAPACK ...
138 Aspen's user avatar
  • 153
4 votes

Reordering eigenvalues in Schur factorization - MATLAB ordschur and LAPACK dtrsen not producing the same results

First of all, please include a definition of Z in your minimal example so that we can reproduce your output. The bug clearly seems to be in the Fortran code. Have you verified that $Z \approx UTU^*$ ...
Federico Poloni's user avatar
4 votes

Step size constraint in Euler backward

This answer is about illustrating the comments on the structure of the Jacobian, the resulting Lipschitz constant and its consequences on the step size. There are only two non-zero entries per row ...
Lutz Lehmann's user avatar
  • 6,109
3 votes

What would be the right approach to numerically integrate in MATLAB an improper integral with Bessel function?

I solved this by observing the behaviour of the rapidly decreasing exponential function. It dropped to zero at about $\xi = 40$. So instead of integrating it from 0 to $\infty$, I integrate it from 0 ...
ishan_ae's user avatar
3 votes
Accepted

How can I find the current for a nonlinear electrical circuit using the Shockley equation, in Octave?

The easiest way is to create a lambda function at call time which lets you set "extra" arguments. You can then call this in a for loop and vary the extra argument. ...
helloworld922's user avatar
3 votes

Matlab eigs function with function handle

The non-function handle version of the eigs function assumes that the A-matrix (and optionally the B-matrix) are defined as matlab sparse matrices and operated on using matlab sparse matrix functions. ...
Bill Greene's user avatar
  • 6,074
3 votes
Accepted

How to remove triangles in a hollow hemisphere shape?

Rather than relying on a delaunay triangulation, you could consider making a structured mesh directly, as you know how you are constructing the points and the ...
Mikael Öhman's user avatar
2 votes

How does the MATLAB backslash operator solve $Ax=b$ for square matrices?

(This answers the question in the original title, but only kind-of, sort-of the question in the original body.) From https://www.mathworks.com/help/releases/R2023b/matlab/ref/mldivide.html, where ...
Kenn Sebesta's user avatar
2 votes

Educational Purpose: How to synchronize chaotic systems

The exact manner in which you couple the systems seems to be the issue here. This can be done for the Lorenz system. Consider the coupled Lorenz systems given by $$ \begin{aligned} \dot{x}_1 &= \...
whpowell96's user avatar
  • 2,546
2 votes
Accepted

Prof Lawrence Shampine's hpde matlab code

There is a zip archive on the Wayback Machine from Shampine’s homepage in 2010: https://web.archive.org/web/20100624204352/http://faculty.smu.edu/shampine/hpde.zip
jdgleeson's user avatar
  • 376
2 votes

ode23, 45, 15s, 15i in matlab for conservative ODEs

None of those ODE solvers are conservative*. For a conservative integrator you typically need a symplectic integrator. There are implementations of various symplectic integrators on Matlab file ...
helloworld922's user avatar
2 votes

How to remove triangles in a hollow hemisphere shape?

I agree with Mikael's suggestion to use a structured mesh if you can get away with it. But there is a more fundamental reason why the Delaunay function from matlab isn't doing what you want and I ...
Daniel Shapero's user avatar
1 vote

How to plot the power spectrum

I pick a column, the sixth, take the FFT, then take the absolute-value squared I guess this is not one PROPER method to estimate power spectrum. The code ...
138 Aspen's user avatar
  • 153
1 vote

Multivariable Newton's method for-loop

It's not entirely clear from the way the question is worder but I guess the idea is given 3 circle equations to find the intersection point of the 3 circles (if it exists). That is: $$\begin{cases} \|...
lightxbulb's user avatar
  • 2,162
1 vote

shooting method to compute the interface shape

Given your differential equation is nonlinear, it is not clear that using the bisection method to find the desired value for $\theta_q$ so that $\theta|_{s = l} = \pi/2$ is reasonable. Let $f(\theta_q)...
spektr's user avatar
  • 4,248
1 vote

'Integral2' error in MATLAB for invalid integrand

If you wish to evaluate this numerically there is no need to involve symbolic math, and we can just write functions; ...
Mikael Öhman's user avatar
1 vote

Is there a Python version of the ODE tool pplane?

I would like to provide with my codes with a result resembling that of matlab, although not 100% the same. Codes are maintained and extended to more general cases via pplane for python. ...
Wei Shan Lee's user avatar
1 vote

How to curve-fit the lower envelope of random sequence?

Well, let our data be $(x_i, y_i)_{i=1}^{n}$ and suppose you have a parametric model $f(x|\beta)$ with parameters $\beta$, that you want to fit to your data. A simple formulation that will help give ...
spektr's user avatar
  • 4,248
1 vote

How to curve-fit the lower envelope of random sequence?

Assuming you have a model $M(x)$ that you can evaluate at $x$ and compare to your data at each point $Y(x)$. Try minimizing $\sum D^2$ where $D = Y(x) – M(x)$ if $Y(x) > M(x)$ and $D = P$ if $Y(x) &...
Ed R's user avatar
  • 11

Only top scored, non community-wiki answers of a minimum length are eligible