Multiplying 4x4 matrices with complicated elements - taking long time

Questions
, , ,
1

I’m very new to Mathematica and to programming. Complete computer novice. I might be coming at my problem from too mathematical a viewpoint

Basically what I’m trying to do is find the transmission of a generic wave through a series of potential barriers with spin orbit interaction and a magnetic field. Now, I can do all the maths, but at a certain point you have to compute these hideous 4x4 matrices so I hoped mathematica could do that part for me and give me graphs.

I want to keep the functions as general as possible until right at the end so I can input parameters right at the graphing stage. But the evaluation gets so big that graphs take hours to plot even over small ranges. Also, the calculation is complicated enough that even if I define 5 out of 6 variables and only keep the variable I want to us as x axis, the evaluation uses up all of my memory and quits before it even finishes.

I don’t understand parallelization - and the little bit of experimenting I’ve done returns <> messages.

The crucial line is

Anz[n_, a_, d_, kx_, l_, g_] = 
 Inverse[L[((n - 1)/2) (a + d), kx]].Smz[a, d, kx, l, 
   g].(MatrixPower[Sz[a, d, kx, l, g], (n - 3)/2]).L[a, kx]

Where L, Smz and Sz are all 4x4 matrices with huge elements. Like so huge I cannot post even one here.

Please help,
My computation as is isn’t even as complicated as it need to be and it takes hours to plot a single graph!

2

Without posting more details I’m afraid there is little that can be helped. We don’t know SMz and Sz and L functions…

Doing MatrixPowers on symbolic expressions creates crazy big and complicated matrices very quickly, i think ‘powering’ numbers (not symbols) makes way more sense…

Please post the notebook!

3

I kinda figured you’d need more info. But I didn’t want to spam everyone right away!
Here is everything I have with annotations. In desperation I have also attached the file. Any tips at all would be appreciated.

P.S. I should say that the end goal of this will be to plot that value of T[n,a,d,kx,l,g] {k,-10,10} without it taking infinite time.

Remove["Global`*"];

$Assumptions = a > 0 && d > 0 && l > 0 && g > 0 && n \[Element] Integers;

Hamiltonian for region with SOI and Zeeman field:

Hsoz[qx_, l_, g_] = {{-g, 0, px + I py, -2 I l}, {0, g, 0, px + I py}, {px - I py, 0, -g, 
  0}, {2 I l, px - I py, 0, g}}

Evals[qx_, l_, g_] = Simplify[Eigenvalues[Hsoz[qx, l, g]]];

Ez1[qx_, l_, g_] = Evals[qx, l, g][[1]];
Ez2[qx_, l_, g_] = Evals[qx, l, g][[2]];
Ez3[qx_, l_, g_] = Evals[qx, l, g][[3]];
Ez4[qx_, l_, g_] = Evals[qx, l, g][[4]];

Evecs[qx_, l_, g_] = Simplify[Eigenvectors[Hsoz[qx, l, g]]];

Normalization below (I suspect there is a better way to do this)

n1[qx_, l_, g_] = g^2 + l^2 + Sqrt[l^4 + g^2 qx^2 + l^2 qx^2];
n2[qx_, l_, g_] = g^2 + l^2 - Sqrt[l^4 + g^2 qx^2 + l^2 qx^2];

Nz1[qx_, l_, g_] = 
  Simplify[1. + (Abs[l]^2*Abs[qx]^2 + Abs[l*(Ez1[qx, l, g] + g)]^2)/
    Abs[Ez1[qx, l, g]*g + n2[qx, l, g]]^2 + 
    Abs[Ez1[qx, l, g]*(n1[qx, l, g] - g^2) + 
       g*(n1[qx, l, g] - g^2 - qx^2)]^2/(
    Abs[Ez1[qx, l, g]*g + n2[qx, l, g]]^2*Abs[qx]^2)];
Nz2[qx_, l_, g_] = 
  Simplify[1. + (Abs[l]^2*Abs[qx]^2 + Abs[l*(Ez2[qx, l, g] + g)]^2)/
    Abs[Ez2[qx, l, g]*g + n2[qx, l, g]]^2 + 
    Abs[Ez2[qx, l, g]*(n1[qx, l, g] - g^2) + 
       g*(n1[qx, l, g] - g^2 - qx^2)]^2/(
    Abs[Ez2[qx, l, g]*g + n2[qx, l, g]]^2*Abs[qx]^2)];
Nz3[qx_, l_, g_] = 
  Simplify[1. + (Abs[l]^2*Abs[qx]^2 + Abs[l*(Ez3[qx, l, g] + g)]^2)/
    Abs[Ez3[qx, l, g]*g + n1[qx, l, g]]^2 + 
    Abs[Ez3[qx, l, g]*(n2[qx, l, g] - g^2) + 
       g*(n2[qx, l, g] - g^2 - qx^2)]^2/(
    Abs[Ez3[qx, l, g]*g + n1[qx, l, g]]^2*Abs[qx]^2)];
Nz4[qx_, l_, g_] = 
  Simplify[1. + (Abs[l]^2*Abs[qx]^2 + Abs[l*(Ez4[qx, l, g] + g)]^2)/
    Abs[Ez4[qx, l, g]*g + n1[qx, l, g]]^2 + 
    Abs[Ez4[qx, l, g]*(n2[qx, l, g] - g^2) + 
       g*(n2[qx, l, g] - g^2 - qx^2)]^2/(
    Abs[Ez4[qx, l, g]*g + n1[qx, l, g]]^2*Abs[qx]^2)];

wz1[qx_, l_, g_] = 
  Simplify[(1./Sqrt[Nz1[qx, l, g]]) Evecs[qx, l, g][[1]]];
wz2[qx_, l_, g_] = 
  Simplify[(1./Sqrt[Nz2[qx, l, g]]) Evecs[qx, l, g][[2]]];
wz3[qx_, l_, g_] = 
  Simplify[(1./Sqrt[Nz3[qx, l, g]]) Evecs[qx, l, g][[3]]];
wz4[qx_, l_, g_] = 
  Simplify[(1./Sqrt[Nz4[qx, l, g]]) Evecs[qx, l, g][[4]]];

Conservation of energy:

qz3x[kx_, l_, g_] = Sqrt[
  g^2 + kx^2 - 2 Sqrt[g^2 kx^2 - g^2 l^2 + kx^2 l^2]];
qz4x[kx_, l_, g_] = Sqrt[
  g^2 + kx^2 + 2 Sqrt[g^2 kx^2 - g^2 l^2 + kx^2 l^2]];

Generic wave equation constructed:

M1[x_, kx_, l_, g_] = 
  Simplify[wz3[qz3x[kx, l, g], l, g]*E^(I qz3x[kx, l, g] x)];
M2[x_, kx_, l_, g_] = 
  Simplify[wz4[qz4x[kx, l, g], l, g]*E^(I qz4x[kx, l, g] x)];
M3[x_, kx_, l_, g_] = 
  Simplify[wz3[-qz3x[kx, l, g], l, g]*E^(-I qz3x[kx, l, g] x)];
M4[x_, kx_, l_, g_] = 
  Simplify[wz4[-qz4x[kx, l, g], l, g]*E^(-I qz4x[kx, l, g] x)];

And in matrix form:

Mz[x_, kx_, l_, g_] = ({
    {M1[x, kx, l, g][[1]], M2[x, kx, l, g][[1]], M3[x, kx, l, g][[1]],
      M4[x, kx, l, g][[1]]},
    {M1[x, kx, l, g][[2]], M2[x, kx, l, g][[2]], M3[x, kx, l, g][[2]],
      M4[x, kx, l, g][[2]]},
    {M1[x, kx, l, g][[3]], M2[x, kx, l, g][[3]], M3[x, kx, l, g][[3]],
      M4[x, kx, l, g][[3]]},
    {M1[x, kx, l, g][[4]], M2[x, kx, l, g][[4]], M3[x, kx, l, g][[4]],
      M4[x, kx, l, g][[4]]}
   });

Hamiltonian for region without SOI and Zeeman field:

L[x_, kx_] = 1./Sqrt[2.] ({{E^(I kx x), 0, -E^(-I kx x), 0},
{0, E^(I kx x), 0, -E^(-I kx x)},
{E^(I kx x), 0, E^(-I kx x), 0},
{0, E^(I kx x), 0, E^(-I kx x)}});

Generalising for potential barriers of width “d” every “a” units of distance

Sl[a_, kx_] = Simplify[L[2 a + d, kx].Inverse[L[a + d, kx]]];

Where it starts to really slow down

Smz[a_, d_, kx_, l_, g_] = Mz[a + d, kx, l, g].Inverse[Mz[a, kx, l, g]];

Sz[a_, d_, kx_, l_, g_] = Sl[a, kx].Smz[a, d, kx, l, g];

Impossible to compute beyond here. Even takes forever if I define every variable:

Anz[n_, a_, d_, kx_, l_, g_] := Inverse[L[((n - 1)/2) (a + d), kx]].Smz[a, d, kx, l, 
   g].(MatrixPower[Sz[a, d, kx, l, g], (n - 3)/2]).L[a, kx]

Tferz[n_, a_, d_, kx_, l_, g_] := ({
   {1, 0, -Anz[n, a, d, kx, l, g][[1, 3]], -Anz[n, a, d, kx, l, g][[1,4]]},
   {0, 1, -Anz[n, a, d, kx, l, g][[2, 3]], -Anz[n, a, d, kx, l, g][[2,4]]},
   {0, 0, -Anz[n, a, d, kx, l, g][[3, 3]], -Anz[n, a, d, kx, l, g][[3,4]]},
   {0, 0, -Anz[n, a, d, kx, l, g][[4, 3]], -Anz[n, a, d, kx, l, g][[4,4]]}})

Auz[n_, a_, d_, kx_, l_, g_] := Inverse[Tferz[n, a, d, kx, l, g]].{{Anz[n, a, d, kx, l, g][[1,1]]}, {Anz[n, a, d, kx, l, g][[2,1]]}, {Anz[n, a, d, kx, l, g][[3,1]]}, {Anz[n, a, d, kx, l, g][[4,1]]}}

Tuz[n_, a_, d_, kx_, l_, g_] := Abs[Auz[n, a, d, kx, l, g][[1]]]^2 + Abs[Auz[n, a, d, kx, l, g][[2]]]^2

Ruz[n_, a_, d_, kx_, l_, g_] := 1 - Tuz[n, a, d, kx, l, g]

Adz[n_, a_, d_, kx_, l_, g_] := Inverse[Tferz[n, a, d, kx, l, g]].{{Anz[n, a, d, kx, l, g][[1, 2]]}, {Anz[n, a, d, kx, l, g][[2, 2]]}, {Anz[n, a, d, kx, l, g][[3, 2]]}, {Anz[n, a, d, kx, l, g][[4, 2]]}}

Tdz[n_, a_, d_, kx_, l_, g_] := Abs[Adz[n, a, d, kx, l, g][[1]]]^2 + Abs[Adz[n, a, d, kx, l, g][[2]]]^2
Rdz[n_, a_, d_, kx_, l_, g_] := 1 - Tdz[n, a, d, kx, l, g]

Tz[n_, a_, d_, kx_, l_, g_] := (Tuz[n, a, d, kx, l, g] + Tdz[n, a, d, kx, l, g])/2.
Rz[n_, a_, d_, kx_, l_, g_] := 1 - Tz[n, a, d, kx, l, g]
4

If you define everything with := rather than =, and just evaluate it numerically (machine numbers), it should be quite fast right? Doing this symbolically will probably be a nightmare anyhow, I quickly inspected the code, if you substitute everything inside each other it will be a HUGE expression. And evaluating this huge expression numerically probably takes longer than just evaluating it numerically directly.
Have you tried?

I would not expect any sensical answer to appear with all these nested definitions.

I see only function definitions, what do you want to do with it in the end?

5

In the end I want to make plots of T[n,a,d,kx,l,g] {k,-10,10} for various values of the other variables (and eventually want to include an angular dependence but I’m a ways from that yet)

I have played around with := for the parts that take a long time to evaluate, but when I try to make a graph the kernel eventually times out and mathematica shuts down.

I don’t quite know what you mean by “numerically” and “numerically directly”

Thank you for responding.

6

When you make the definitions, like you do know, it will evaluate it already and thus ‘copying/replacing’ the definition that occur in it. Since your definitions are very much nested, it will nest all these functions inside each other symbolically:

Compare the two ways of defining functions (immediate and delayed assignment):

ClearAll[a,b,n]
b[n_]:=n^2
a[m_]:=m b[m]
DownValues[a]
DownValues[b]

ClearAll[a,b,n]
b[n_]=n^2
a[m_]=m b[m]
DownValues[a]
DownValues[b]

Have you tried to evaluate T for some suitable numbers when you have all your definitions done using := (this way the kernel will not die of memory requirements because it does not substitute all the sub-definitions inside each other)

7

To add, and please put the numbers with a period in it: 3.0 (not 3) et cetera. So you will trigger the usage of machine numbers rather then exact arithmetic.

8

Yes, I have tried to evaluate T for suitable numbers. It takes, I would say, over 30 minutes the way the code is written now. And that’s for a single value of kx.

I know about putting periods in numbers. That sped things up a lot when I was trying this for a much simpler case.

Are you saying I should hold off the evaluation of everything in the book until the end?

9

Yes, do all the definitions using := to stop Mathematica from substituting everything inside each other… Could you give me some ‘typical’ values for the function T you want to calculate?

10

Typical values I’m trying now

n=5. (this really should be exactly 5 but if this runs faster 5.0 will do)

a=0.23

d=0.12

kx= 1. (this is what I want to plot over so 1. is fine for now)

l=0.01

g=0.09

11

and the function T itself would also be useful :wink:

12

Oh sure. Sorry. In this workbook it’s Tz[n,a,d,kx,l,g]

I was thinking about an older workbook. I want to plot Tz over kx

13

In your function L[x_,kx_]… your exponent are k*x not kx.
and
SL should be:

Sl[a_, kx_, d_] := Simplify[L[2 a + d, kx].Inverse[L[a + d, kx]]]

I presume?

If I remove all the simplify commands, and I run:

Table[Tz[5.0, 0.23, 0.12, k, 0.01, 0.09], {k, DeleteCases[Range[-10.0, 10.0], 0.0]}]

it takes roughly takes 15 seconds. But there might be still an error as it return a 1D vector, not a number. Is that correct. I can further optimize it a lot, because you recalculate Evals and Evecs many times.

14

kx! Yes. I caught that mistake as I was sending you the notebook but I didn’t change it.

I suppose it should be a 1D vector cos that is for 20 values of k (?)

What is the DeleteCases command?

Other optimizations would be great if you have time. I am such a novice at this. I’m literally transposing a handwritten work-through of this problem.

15

DeleteCases does what is kind-a says: it delete certain cases, in this case 0.0 is removed from the list, as Tz blows up for k=0.
No i meant that ‘Tz’ returns a vector (of 1 value) not just a value. is that correct?

16

Hmmm. I guess it should just be a value. It’s the transmission probability. It has to be between 0 and 1. But in principle it should just be a number.

17

I sped it up another factor 10 just now, but there is still lots to go. I’m sure.

18

Ooh. That will plot Tz.

Unfortunately I just noticed an error in the notebook I sent you

Anz[n_, a_, d_, kx_, l_, g_] := 
 Inverse[L[((n - 1)/2) (a + d), kx]].Smz[a, d, kx, l, 
   g].(MatrixPower[Sz[a, d, kx, l, g], (n - 3)/2]).L[a, kx]

should read

Anz[n_, a_, d_, kx_, l_, g_] := 
 Inverse[L[((n - 1)/2) (a + d), kx]].Smz[a, d, kx, l, 
   g].(MatrixPower[Sz[a, d, kx, l, g], (n - 3)/2]).L[a, kx]

(I had tried to make a simplification that didn’t work.)

Changing that the Table[] still evaluates very quickly. And the graph of Tz comes out in 60 seconds.

It is, however not between 1 and 0 which is troubling. I may have more errors than I thought.

19

I don’t see a difference!

20

That is because I copy and pasted the same thing twice!

Should have read:

Anz[a_, d_, kx_, l_, g_] := 
 Inverse[L[(a + d), kx]].Smz[a, d, kx, l, g].Sz[a, d, kx, l, g].L[a, 
   kx]

should read:

Anz[n_, a_, d_, kx_, l_, g_] := 
 Inverse[L[((n - 1)/2) (a + d), kx]].Smz[a, d, kx, l, 
   g].(MatrixPower[Sz[a, d, kx, l, g], (n - 3)/2]).L[a, kx]

Sorry about that