On Fri, Aug 14, 2020, 12:05 AM Stephen J. Turnbull <turnbull.stephen.fw@u.tsukuba.ac.jp> wrote:
I cannot speak to engineering uses of matrix computations.  If someone
produces use cases that fit into the list of operations proposed
(addition, negation, multiplication, inverse, transposition, scalar
multiplication, solving linear equation systems, and determinants, I
will concede it's useful and fall back to +/- 0.

I'll try to speak to that below.  

On Fri, Aug 14, 2020, 1:48 AM Marco Sulla <Marco.Sulla.Python@gmail.com> wrote:
On Fri, 14 Aug 2020 at 03:16, Steven D'Aprano <steve@pearwood.info> wrote:
> Funny you mention this, I have been working on a Matrix object for
> precisely the use-case you discuss (secondary school maths), where
> performance is not critical and the dimensions of the matrix is
> typically single digits.

This is really good. I think that numpy is a practical project, but I
feel it hard to understand for people that come from Python and not
from Matlab.
Maybe the API for such a module can be released initially as a
provisional API, as at the pathlib and asyncio beginnings?

I agree this could be great. I would characterize the steps to get to matrix arithmetic in python as currently something like:

1. Learn python, with typical experience of:
- I learned it last night! Everything is so simple!
2. Learn numpy, with the typical experience of:
- Wait, is this array thing really what I need? I am looking for matrix. 
- This is kind of intimidating and after a few hours I barely feel confident I can add matrices together correctly... also, why did this tutorial cause me to spend half a day figuring out virtual environments? And this other tutorial doesn't even mention them? And this other one says never to install things without using them? Did I hurt my computer by putting numpy on it?

It's tough.

However, engineers aren't usually doing a lot of linear algebra. Linear algebra isn't even a required course in ABET accredited university engineer programs. I walked out of my linear algebra course the first day, saying to myself, "nope, I have written my last proof."

What is used is mostly NOT even matrix arithmetic of the high school variety (although it's occasionally useful to solve a system of equations). I'll give a longish example below.

The other thing I would say is there are some nice third party tools out there that assume numpy is being used (I use pint a lot for units, and of course there's notebooks like jupyter). Might be a good thing to see if the API for basic matrix objects can be made to work well with some of these existing tools.  

Here's my example code-- I utilize a lot of this sort of thing. Skip down this large block of code to see the lines that exemplify the reason I'm using a numpy array in the first place.

#  Nominal skin friction capacity, Pn, of drilled shafts/piers

## depths below grade to 20 ft in 1 foot increments
z = np.linspace(1, 20, num=20)

## create array containing unit weight of soil at each depth

### unit weight log
### for looking up unit weights (pcf) based on depth (ft)
### of form {depth: unit_weight}
log = {3: 120, 9: 118, 12: 122, 15: 106, 18: 102, 21: 110}

### fill in gaps in unit weight log
#### assume first layer extends to ground surface
log[0] = log[next(iter(log))] 
#### assume unit weight change occurs after mid point between each layer
log_depths = sorted(list(log.keys()))
log_filled = {}
for depth, next_depth in zip(log_depths[:-1], log_depths[1:]):
    # layer boundary
    log_filled[depth] = log[depth]
    # layer change occurs after midpoint of boundary
    log_filled[ (depth+next_depth)/2 ] = log[depth]
#### don't forget bottom of last layer
log_filled[next_depth] = log[next_depth]
#### final per foot unit weights, γ
γ_per_ft = np.array([])
for depth in z:
    log_depth = next(d for d in log_filled if d >= depth)
    γ_per_ft = np.append(uwt_per_ft, log_filled[log_depth])

## calculate per foot stress

### unit weight of water (pcf)
γ_w = 62.4

### effective unit weight given high water table
γ_prime = γ_per_foot - γ_w 

q_s = γ_prime * z

## Nominal capacity

### shaft diameter (ft)
D_shaft = 6

### shaft circumference
C_shaft = D_shaft * np.pi

# friction coefficient
friction_angle = 17  # degrees
μ_s = np.tan( friction_angle * np.pi/180 )

# nominal vertical capacity
P_n_per_foot = μ_s * q_s * C_shaft
P_n = np.sum(P_n_per_foot)

This operation:

q_s = γ_prime * z

And these operations:

P_n_per_foot = μ_s * q_s * C_shaft
P_n = np.sum(P_n_per_foot)

...which are just scalar or element-wise multiplication and summation, are pretty much the extent of the kinds of mathematical things an engineer is doing with numpy arrays- at least the majority of the time. We just generally are not doing linear algebra.

That being said, I do think it would be nice to have a "matrix" type in the std lib that supported both element-wise operations like this and basic matrix arithmetic.