3
This problem is too small to actually be sparse. Sparse handling has a big overhead because the indexing is not "direct", i.e. you don't necessarily know where the next value will be without branch checking. So you need it to be "sparse enough" that the O(n^3) dense LU-factorization cost shrinking to the purely non-zero terms overcomes ...
3
The effect is due to quantization noise and/or aliasing, you are not computing with the true minima of the radius along the orbit, but with the closest sample point. This means that the sample point for the half-step integration can lie in-between. Close to the minimum this results in an $O(e·dt^2)$ ($e$ the eccentricity) difference between the computed ...
2
One the most useful things I found myself when converting from matlab to numpy/scipy was this page NumPy for MATLAB users. I even look at this every once in awhile to remind myself even to this day. I also highly recommend taking a look at the matplotlib gallery as what you probably want to do visually in either MATLAB or Mathematica is hopefully ...
2
You have not completely vectorized your code. z(1:10) is the first 10 numbers in the data of z, not the first 10 rows. For that you need to write z(1:10,:) like you have already correctly done on the left side.
Why the flip in the matrix orientation occurs between the two calls, that is why z(1:10) inherits the column-ness in the first case and the row-ness ...
2
odenumjac calls your function in a vectorized manner it seems, and your function is not vectorized. You can easily change that by changing the second index of f in your function to : instead of 1, for instance:
f(10,:) = 2*(x(end-1,:) - x(end,:));
I thought the setting joptions.vectvars=1 would not allow the vectorised call (see one of your other questions). ...
1
Since you already know Matlab and Mathematica, I would go with SciPy Lecture notes. It is a course on using Python for Scientific purposes.
Gael Varoquaux, Valentin Haenel, Pierre de Buyl, Gert-Ludwig Ingold, Emmanuelle Gouillart, Michael Hartmann, … João Felipe Santos. (2017, October 4). scipy-lectures/scipy-lecture-notes: Release 2017.1 (Version 2017.1). ...
1
The simple example as posted wasn't stiff, but I put it together in Julia anyways for show. I modified it to be a system of 2 PDEs with N=1200 and get in the 10's of ms
using ModelingToolkit, LinearAlgebra, BenchmarkTools, DifferentialEquations
# Setup matrices
N = 1200
mat1 = hcat(zeros(N),Tridiagonal(ones(N-1),-2*ones(N),ones(N-1)),zeros(N))
mat1[1,1] = 1
...
Only top voted, non community-wiki answers of a minimum length are eligible
