# Tag Info

61

Except for code which does a significant number of floating-point operations on data that are held in cache, most floating-point intensive code is performance limited by memory bandwidth and cache capacity rather than by flops. $v$ and the products $Av$ and $Bv$ are all vectors of length 2000 (16K bytes in double precision), which will easily fit into a ...

37

If you want to see what a\b does for your particular matrix you can set spparms('spumoni',1) and figure exactly what algorithm you were impressed by. For example: spparms('spumoni',1); A = delsq(numgrid('B',256)); b = rand(size(A,2),1); mldivide(A,b); % another way to write A\b will output sp\: bandwidth = 254+1+254. sp\: is A diagonal? no. sp\: is band ...

37

In Matlab, the ‘\’ command invokes an algorithm which depends upon the structure of the matrix A and includes checks (small overhead) on properties of A. If A is sparse and banded, employ a banded solver. If A is an upper or lower triangular matrix, employ a backward substitution algorithm. If A is symmetric and has real positive diagonal elements, ...

24

Disclaimer: I sometimes get annoyed when somebody tries to tell me what they think I ought to do rather than answering the question I asked. But I'm going to take a risk and suggest an alternative to you. I would suggest looking at Python's scientific computing packages: numpy, matplotlib, and scipy. Together, they provide you most of the core ...

20

Going from MATLAB to Python does introduce quite a bit of syntax overhead. One way to quantify it is the nice QuantEcon cheatsheet which showcases how there's a lot of extra "stuff" going on when trying to write simple linear algebra commands in Python. The verbose NumPy syntax is really just a symptom of how it was not developed as a technical computing ...

19

Your code is limited by memory bandwidth. For trivial math, it's often better to count memory accesses rather than flops. You'll get the following table: operation memory reads/writes matrix + matrix 3n² matrix * vector 2n²+n (if vector is not cached) matrix * vector n²+2n (if vector is only read once) vector + vector ...

18

The most important aspect of interpolation and curve fitting is to understand why high order polynomial fits can be an issue and what the other options are and then you can understand when they are/are not a good choice. A few issues with high order polynomials: Polynomials are naturally oscillatory functions. As the order of the polynomial increases, the ...

18

Short answer, you want to have the leftmost index on the innermost loop. In your example, the loop indices would go k, j, i and the array indices would be i, j, k. This has to do with how MATLAB stores the different dimensions in memory. For more, see #13 of this reddit post.

17

GNU Octave is "mostly compatible with Matlab", certain subtleties means not all scripts are portable from MATLAB to Octave. It is worth reading the documentation for the language and/or compatibility notes in the FAQs or on wikibooks. There are also porting notes. Packages similar to MATLAB toolboxes exist, but you will need to check them out to work out ...

17

Matlab interprets sequences of multiplications and/or divisions from left to right. Hence $A*B*C*v$ is much more expensive than $A*(B*(C*v))$, as you have two matrix products and one matrix-vecor product in place of three matrix-vector products. On the other hand, $A*(B*(C*v))$ should be slightly faster than if you save the intermediates in separate ...

16

For sparse matrices, Matlab uses UMFPACK for the "\" operation, which, in your example, basically runs through the values of a, inverts them, and multiplies them with the values of b. For this example, though, you should use b./diag(a), which is a lot faster. For dense systems, the backslash-operator is a bit more complicated. A brief description of what is ...

16

Both of them are direct solver to solve linear systems (opposing to iterative solver). mldivide does perform the tests for $A$ in solving $Ax = b$. Please see Allan's answer in this thread for more information. Also see MATLAB's help on mldivide algorithm here. mldivide for square matrices: If A is symmetric and has real, positive diagonal elements, ...

16

Note that $\pi/2$ is represented in double precision format in a way that is not exactly equal to $\pi/2$. It's only accurate to about 15 digits. Thus you're starting every so slightly away from the equilibrium position. Since the equilibrium is unstable, it will eventually start moving.

15

A difficulty with any of these types of questions is that the answer is highly community-dependent. To answer some of your questions in haphazard order: MATLAB is used a lot both in academia and in industry. One of the reasons it's used quite a bit in industry is because it is taught in academia. I know for a fact that MATLAB is used at Lincoln Laboratory ...

15

The function error('error message'); will exit your program and print the error message to the console.

15

The question has two very different subquestions. I will address the first one only. Matlab's version runs on average 24 times faster than my python equivalent! The second one is subjective. I would say that letting know the user that there is some problem with the integral is a good thing and this SciPy behavior outperforms the Matlab`s one to keep it ...

15

I think the two main points have already been made by Brian and Ertxiem: your initial value is an unstable equilibrium and the fact that your numerical computations are never really exact provides the small perturbation that will make the instability kick in. To give a bit more detail how this plays out, consider your problem in the form of a general ...

14

For starters, I wouldn't use intermediate variables, but brackets. Unless, of course, you're interested in the intermediate results, but I'm guessing not. I tried the following in Matlab: >> N = 500; >> A = rand(N); B = rand(N); C = rand(N); v = rand(N,1); >> tic, for k=1:100, A*B*C*v; end; ...

14

Since the matrices are so small, all of the cost is going to be in call overhead. If you will do the transformation many times, it will be faster to precompute D=A*B*C once and then for each vector apply v_f=D*v_i. You could also consider bringing this out to a mex file.

14

I personally have come the way from Gnuplot to Matplotlib with PGFPlots as an intermediate step. I will try to name all aspects of Matplotlib that I like. It is very versatile. You are not limited to line or scatter plots, you can easily do bar plots, images (matrix visualization!), basic 3D plotting and even some animation. You can use Matplotlib as a GUI ...

14

The most obvious thing you can do is to precompute [L,U] = lu(A) ~ O(n^3) Then you just compute x = U \ (L \ b) ~ O(2 n^2) This would reduce the cost enormously and make it faster. Accuracy would be the same.

14

Here is R1, as computed in MATLAB: 1.0e+07 * -7.382605957465515 -9.599867106092937 -2.830412177259742 -0.000000000002830 -0.000000000002830 -1.230434326244253 -1.599977851015490 -0.471735362876624 -0.000000000000472 -0.000000000000472 3.691302978732758 4.799933553046468 1.415206088629871 0.000000000001415 0.000000000001415 -5....

14

First, see Mark L. Stone's answers, which is completely correct. Second, realize that this is the reason why people told you to use relative errors in your numerical analysis class. :) Third, the real question here is why the results do not coincide exactly, since both languages call some BLAS library functions for their computations. There are several very ...

13

Computing the term $e^x$ is definitely significantly more expensive than computing a lower-order polynomial -- say $x^4$. But it may be ten to 100 times more expensive at most, not "crazy" expensive. So I suspect that if Matlab takes forever to compute something, then that's because the character of your ODE changes significantly. For example, the presence ...

13

There are libraries that you can use in Python that will give you all (or at least nearly all) of the functionality of MATLAB. For example, scipy.integrate.solve_ivp() supports a number of methods for ODE integration that are comparable to what you can get with the various odexxx() functions in MATLAB. So no, you wouldn't have to write your own ODE ...

12

Yes you can learn MATLAB via Octave. But Octave syntax is less restrictive and more in line with modern scripting languages. MATLAB seems behind in this respect. See this wiki link MATLAB Programming/Differences between Octave and MATLAB Another major difference to me were the availability of certain libraries for MATLAB but not for Octave.

12

I've written my own integrator, quadcc, which copes substantially better than the Matlab integrators with singularities, and provides a more reliable error estimate. To use it for your problem, I did the following: >> lambda = 0.00313; kappa = 0.00825; nu = 0.33; >> x = 10; >> E = @(r) r.^4.*(lambda*sqrt(kappa^2 + r.^2)).^(-nu-5/2) .* ...

11

As Pedro points out, Levin-type methods are the best established methods for these kinds of problems. Do you have access to Mathematica? For this problem, Mathematica will detect and use them by default: In:= e[r_] := r^4 (l Sqrt[k^2 + r^2])^(-v - 5/2) BesselK[-v - 5/2, l Sqrt[k^2 + r^2]] In:= {l, k, v} = {0.00313, 0.00825, 0.33}; In:= Block[{...

11

As Misha and Geoff Oxberry pointed out, Mathematica really has a different focus (just because you can pound in a nail with a screwdriver doesn't mean you should). So I take your question as being "If I know Matlab, why should I learn Python?" [Edit: and so, apparently, did you.] For all intents and purposes, Matlab is the English of scientific computing -- ...

11

MATLAB has a couple of "exact" functions for this, cond and rcond, with the latter returning a reciprocal of the condition number. Matlab approximate function condest is more fully described below. Often estimates of the condition number are generated as by-products of the solution of a linear system for the matrix, so you might be able to piggyback the ...

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