Group Abstract Group Abstract

Message Boards Message Boards

Multiplying 4x4 matrices with complicated elements - taking long time

Eleanor Holmes
Posted 10 years ago
POSTED BY: Eleanor Holmes
35 Replies
Sort By:
Eleanor Holmes
Eleanor Holmes
Posted 10 years ago
Attachments:
POSTED BY: Eleanor Holmes

I am far from an expert, but I know that as a rookie one can be far from optimal. On the other hand, it could also be that the problem you are trying to do is just too hard (it certainly looks that way). I will look at your notebooks some more. I was hoping to see something besides code. I am not very familiar with condensed matter.

Best,

OL.
POSTED BY: Otto Linsuain
POSTED BY: Otto Linsuain
Eleanor Holmes
Eleanor Holmes
Posted 10 years ago

Thanks for replying. I don't know how far down the thread you have read.

It's the basic square potential well problem. My tweaks are that I want a generic wave including spin, arbitrarily many identical barriers uniformly spaced ("a" apart and "d" wide) in graphene, with spin orbit interaction and a zeeman field under the potential.

My approach is to first try to replicate some simplified systems (only spin orbit, only one barrier, normal incidence etc) and then to generalise entirely so I can have any number of barriers at any incidence angle and with both SOI and zeeman.

I understand the maths of it (for the most part) but I am having huge trouble translating it to mathematica language cos I am far from fluent. I am making all of the rookie mistakes and it's costing me computing time, memory, and accuracy!
POSTED BY: Eleanor Holmes

Delayed evaluation works for me, but only if you delay everything but the function that rely on Eigensystem commands...
POSTED BY: Sander Huisman
Eleanor Holmes
Eleanor Holmes
Posted 10 years ago

You're right! Delaying everything but the eigenvalue calculations (which are always pretty quick) works great!

Now I need to generalise my equations so I can have arbitrarily many barriers. I'm getting an error I don't recognise.

(apologies if this is a duplicate reply, the thread seems to be bugging out for me)

Attachments:
POSTED BY: Eleanor Holmes
Eleanor Holmes
Eleanor Holmes
Posted 10 years ago

<<Progress Update>>

Now I have successfully expanded this for two barriers (or a square potential well) and am attempting to expand fro arbitrarily many barriers. Getting errors I don't understand though. Also still can't delay evaluation without getting the wrong result.

Attachments:
POSTED BY: Eleanor Holmes
Eleanor Holmes
Eleanor Holmes
Posted 10 years ago

The symbolic result - I think so yes.

So I got it to plot the graphs I want. (Well almost. rup and rdown are inverted for k<1 for some reason, but that isn't so important.

But, the way I have it obviously takes more time than is necessary because I make it evaluate symbolically before telling it numbers. As you have seen from my more complex example, I really can't afford to lose that optimization.

How would I have to rewrite this in order to still get those same graphs out in the end?

Attachments:
POSTED BY: Eleanor Holmes
Eleanor Holmes
Eleanor Holmes
Posted 10 years ago

Delayed evaluation works for me, but only if you delay everything but the function that rely on Eigensystem commands...
POSTED BY: Eleanor Holmes
Eleanor Holmes
Eleanor Holmes
Posted 10 years ago

Hmmm - how do I avoid this problem then?
POSTED BY: Eleanor Holmes

I'm not quite sure, but I guess generally it has to do with cases there the numbers become '0' and you have to do some division. So is it possible for your elements in your matrix to be zero? Just check what you expect, you can put some values in there and try. I put 3 random values for qx qy and l in there...

Also you can check the symbolic result, is that what you want?
POSTED BY: Sander Huisman
Eleanor Holmes
Eleanor Holmes
Posted 10 years ago

Ok. I've had to go back to a simpler example to figure out what the errors are in the more complex one.

So the (almost*) graphs I want are the first three which I get when I evaluate everything separately without delaying (no :=). The graphs I get using := are the second three which are wrong.

I am also having problems explaining to mathematica the a vector k should be k = Sqrt[kx^2 + ky^2] and that ArcTan[ky/ky] is \thetak therefore kx + i ky = k e^(i \thetak) - which is a very longhand maths way of thinking but I understand that and I don't understand computers.

I also can't figure out how to tell mathematica that at the barrier (x=0) then E3[q3] = k = E4[q4]

*these graphs are actually still wrong. But I don't know why. It should be the exact same curve but between 0 and 1 on the y axis.

Attachments:
POSTED BY: Eleanor Holmes

Have a look at the following functions:

AbsArg
AngleVector
ReIm
ArcTan[x, y]

The first gives the magnitude and the 'angle' (argument) of a complex number at once. The second can construct a vector by a radius-angle pair The third gives you real and Imaginary part of a number as a pair. And perhaps use Arctan with two arguments; so it takes care of the 4 quadrants (the signs of x and y et cetera).

The problem is in your eigensystem function:

enter image description here

First symbolic finding the general eigenvalues and eigenvectors may involve certain assumptions, not valid for all cases! Take for example the simple case of this polynomial:

enter image description here

In the general case 3x^2/x is 3x but not for x=0 ! So 'underneath' the symbolic calculation there are always some assumptions (unavoidable in some sense). I'm not sure why the general solution of the eigensystem 'works' for your input but not the non-general case!
POSTED BY: Sander Huisman
Eleanor Holmes
Eleanor Holmes
Posted 10 years ago
POSTED BY: Eleanor Holmes

Once you have some sensible answers from all your sub-functions and you need a big speed up, just continue this thread. I will be notified and others will also see it. I guess it can done in a millisecond or so, and when parallelised even faster of course...

Yes, please use Normalize as it also deals with the case {0,0,0} which (when not programmed properly) blows up if you just divide by the length of the vector...
POSTED BY: Sander Huisman
Eleanor Holmes
Eleanor Holmes
Posted 10 years ago
POSTED BY: Eleanor Holmes

I don't see a difference!
POSTED BY: Sander Huisman
Eleanor Holmes
Eleanor Holmes
Posted 10 years ago

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
POSTED BY: Eleanor Holmes
Attachments:
POSTED BY: Sander Huisman
Eleanor Holmes
Eleanor Holmes
Posted 10 years ago
POSTED BY: Eleanor Holmes

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

Attachments:
POSTED BY: Sander Huisman

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!
POSTED BY: Sander Huisman
Eleanor Holmes
Eleanor Holmes
Posted 10 years ago
Attachments:
POSTED BY: Eleanor Holmes
POSTED BY: Sander Huisman
Eleanor Holmes
Eleanor Holmes
Posted 10 years ago

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.
POSTED BY: Eleanor Holmes
POSTED BY: Sander Huisman
Eleanor Holmes
Eleanor Holmes
Posted 10 years ago

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?
POSTED BY: Eleanor Holmes

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?
POSTED BY: Sander Huisman
Eleanor Holmes
Eleanor Holmes
Posted 10 years ago

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
POSTED BY: Eleanor Holmes

and the function T itself would also be useful ;)
POSTED BY: Sander Huisman
Eleanor Holmes
Eleanor Holmes
Posted 10 years ago

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
POSTED BY: Eleanor Holmes

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.

Attachments:
POSTED BY: Sander Huisman
Eleanor Holmes
Eleanor Holmes
Posted 10 years ago

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.
POSTED BY: Eleanor Holmes

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?
POSTED BY: Sander Huisman

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.
POSTED BY: Sander Huisman
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Tag limit exceeded
Note: Only the first five people you tag will receive an email notification; the other tagged names will appear as links to their profiles.
Reply Preview
Attachments
Remove
or Discard