57
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 ...
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 ...
18
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 ...
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 ...
5
The particular set of constraints you have chosen does not prevent a rigid body rotation about node 1. Thus the stiffness matrix is singular, as you have noted.
One way to prevent this rigid body rotation is to set the y-displacement at node
2 to zero. You could also constrain the x-displacement at either node 3 or node 4 to prevent the rotation.
One way ...
4
Crank-Nicolson is a very good classical approach for parabolic PDE like
the heat transfer PDE to which it was originally applied. It is relatively
easy to understand and implement so it is often presented in basic courses
on numerical methods for PDE. pdepe is also very well-suited to this class of PDE
(the second "p" in pdepe stands for parabolic). It has ...
4
First of all, MATLAB's gmres assumes that the preconditioner you use is linear. This is important! Actually it is the main difference between FGMRES and GMRES. Right preconditioned GMRES and FGMRES are exactly the same if you use a linear preconditioner, however, FGMRES allows the use of non-linear preconditioners.
What do I mean by a non-linear ...
4
Standard examples of PDE to solve with the typically taught basic discretization methods (Crank-Nicolson et al.) are
Transport equations, and other first order equations like Burger's, have often explicit solutions and conservation laws that the numerical methods more-or-less satisfy
The heat equation with different boundary conditions and source terms is ...
4
As mentioned in a comment by @AloneProgrammer, it's unlikely that your system has a solution. But to make any kind of progress, it's also useful to rewrite the system in a way that makes its structure simpler. To this end, notice that in the first equation, you always have terms of the form $A_{ij}=B_{ij}/V_i$. So the first of the two equations might be ...
2
You should use a modeling language so that your code is independent of the underlying solver.
cvxpy is a good choice.
When I rewrite your model in cvxpy:
#!/usr/bin/env python3
import cvxpy as cp
import numpy as np
from numpy.random import normal as randn
Sample = 10
H = randn(size=(4,2,Sample))+1j*randn(size=(4,2,Sample))
h = randn(size=(4,1))+1j*randn(...
2
If you're interested in modelling any type of PDE within MATLAB, the Partial Differential Equation Toolbox should be able to handle anything you're interested in.
The complete documentation for the toolbox can be found here. .
A suggested workflow for some simple examples can be found here.
The solvable equations via the toolbox are described in detail here....
2
In essence, what you want to do is interpolate data: You want to find a function $U(x,y)$ so that $U(x_i,y_i)=u_i$. This is then the function for which you'd like to do draw isocontours.
You should take a look at books on numerical methods and search for "two-dimensional interpolation". Among the methods you will find is to triangulate your points $(x_i,y_i)...
1
The only way I know of to visualise such data is to first triangulate it. Matlab can do the triangulation. See the Matlab functions 'delaunay', 'trisurf', 'pdecont'. This works if the domain is convex. If not then you need to construct the triangulation by accounting for the domain boundaries. Maybe Matlab has methods to do this, I am not very much familiar ...
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. ...
Only top voted, non community-wiki answers of a minimum length are eligible
