New answers tagged python
2
As I mentioned in my comment, due to that you are searching for a method based on Python and ideally available in NumPy or SciPy, my suggestion is to use scipy.linalg.expm. It's really easy to use basically you have the $X$ matrix and you just pass it to the scipy.linalg.expm class and it would give you the exponential of $X$ (i.e. $e^{X}$).
Due to that you ...
1
Learning Python allows you so solve many text-file parsing/processing and manipulation tasks (it is universal language).
But you should write the code in a way it will be useful to the team even in the case you will leave. So if nobody else uses Python with the libraries, it would be a bad taste to become the rogue Python coder. Ask the group leadership. ...
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 ...
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 ...
1
It turns out the issue was very simple (and not related to a limitation of numpy's matrix power function as such). I initially thought there was some numerical floating point error being propagated - but the example matrix I was testing on contains only integers. The problem was that at $k=20$ the value of the matrix entries exceeded numpy's maximum possible ...
1
The Frobenius norm is not an operator norm, it is a norm on the vector space of linear operators/matrices, which is not the same thing. Just change it to any other preset norm and it should work.
It is also the case that your method of computing matrix powers is not stable. The algorithm used in Numpy is basic repeated squaring, which has no normalization ...
5
I see at least one important problem. On the right hand side you have a term that looks like
$$P_o \left( \frac{\dot{R}}{R} \right)^{3 \kappa}$$
This term is dimensionally inconsistent with the other terms in the brackets, which have dimensions of a pressure. This term should actually be
$$P_o \left( \frac{R_o}{R} \right)^{3 \kappa}$$
The paper you link ...
Top 50 recent answers are included
