# HP Forums

Full Version: Using the Prime to solve for eigenvalues
You're currently viewing a stripped down version of our content. View the full version with proper formatting.

The Prime has an eigenvalue solver as part of it's matrix utilities that will solve the common eigenvalue problem Det (A - lambda) = 0, where A is a square matrix and lambda is the eigenvalue. The problem I have, however, is in the form Det (A - lambda * B), where B is a diagonal square matrix, which does not work with the eigenvalue solver. I could convert it to the base form only if all the elements of B were the same. In the past under RPN and RPL, I used the root solver to do this by creating a function Det (A - lambda * B) = 0 to solve for lamda, but it doesn't seem possible to do this with the Prime solver. Perhaps someone can show me how to do this on the Prime.

Edited: 23 Oct 2013, 12:44 p.m.

Try this in CAS mode:

solve(DET(A - x*B)=0,x) // or cSolve

Or

zeros(DET(A-x*B),x) // or cZeros

You can store the matrices before solving. Hope this helps, Eddie

Thanks a lot, Eddie ! It works perfectly. One thing I discovered is that I had to uncheck "Return exact symbolic results" in the CAS settings for it to work. Also, if you insert the expression as an equation i.e. expression = 0, then you don't need to include the second variable parameter i.e. Solve(expr(x)=0) is the same as Solve (expr(x)=0,x)

OK, so now I have another question. I want to automatically create a matrix using a program that first inputs spring values k(i) into a vector V of size N, and then generates the square matrix K of size NxN that starts out populated with all zeros and then makes a tri-diagonal matrix based on the rules:

K(i,i-1) = -V(i) ; i-1>=1
K(i,i) = V(i)+V(i+1) ; i+1 <= N
K(i,i+1) = -V(i+1) ; i+1 <= N
K(N,N) = V(N)

Under RPL this was very easy to do with the GET and PUT commands, but I don't see anything in the Matrix programming commands to do this. I can input the vector values, but how do I insert them per the rules above into the matrix.

Quote:
OK, so now I have another question. I want to automatically create a matrix using a program that first inputs spring values k(i) into a vector V of size N, and then generates the square matrix K of size NxN that starts out populated with all zeros and then makes a tri-diagonal matrix based on the rules:

K(i,i-1) = -V(i) ; i-1>=1
K(i,i) = V(i)+V(i+1) ; i+1 <= N
K(i,i+1) = -V(i+1) ; i+1 <= N
K(N,N) = V(N)

Under RPL this was very easy to do with the GET and PUT commands, but I don't see anything in the Matrix programming commands to do this. I can input the vector values, but how do I insert them per the rules above into the matrix.

Let M0 be the matrix you want to create and M1 be the vector entered.

Assigning an element from M1 to M0:

M0(row, col):=M1(element)

You can also use the STO> arrow like this: M1(element)>M0(row, col)

Hope this helps you get started. Eddie

Thanks again, Eddie. Here's the program I've written. It's sloppy, but gets the job done. Note the commented lines at the end. I'm trying to run the solver within the program, but it crashes the calculator. It works fine in CAS, but I can't get it to work within a program. Any suggestions ?

```EXPORT Frequencies()
BEGIN
LOCAL I,J,N,X;
INPUT(N,"NDOF","N =","Enter N",0);
M7:=MAKEMAT(0,N);
M8:=MAKEMAT(0,N,N);
M9:=MAKEMAT(0,N,N);
FOR I FROM 1 TO N DO
INPUT(X,"Mass","Mass =","Enter Mass",0);
M8(I,I):=X
END;
FOR I FROM 1 TO N DO
INPUT(X,"Spring","Spring =","Enter Spring",0);
M7(I):=X
END;
FOR I FROM 1 TO N DO
FOR J FROM 1 TO N DO
IF I == J AND I < N THEN M9(I,J):=M7(I)+M7(I+1)
END;
IF I == J-1 THEN M9(I,J):=-M7(J)
END;
IF I == J+1 THEN M9(I,J):=-M7(I)
END;
END;
END;
M9(N,N):=M7(N);
//M7:=solve(DET(M9-X*M8)=0)
//PRINT("Frequencies = " + sqrt(M7)/(2*pi));
END;
```

It is not the easiest working with symbolic objects in Prime Programming. What I did was make a temporary object for DET(M9 - 'X'*M8) and then called up the zeros function to solve for 'X'. In programming (and in Home Mode), the CAS. must precede any CAS operation.

I chose zeros instead of solve, since CAS.zeros works in home mode, and CAS.solve doesn't do as well.

One more point: I suggest modifying input strings to make program clarify which mass and spring to ask for (Mass 1, Mass 2, etc.).

Program:

```EXPORT Frequencies()
BEGIN
LOCAL I,J,N,X,temp;
// Updated 10-25-2013
INPUT(N,"NDOF","N =","Enter N",0);
M7:=MAKEMAT(0,N);
M8:=MAKEMAT(0,N,N);
M9:=MAKEMAT(0,N,N);
FOR I FROM 1 TO N DO
// Modified Input Strings
INPUT(X,"Mass "+I,"Mass "+I+" =","Enter Mass",0);
M8(I,I):=X
END;
FOR I FROM 1 TO N DO
INPUT(X,"Spring "+I,"Spring "+I+" =","Enter Spring",0);
M7(I):=X
END;
FOR I FROM 1 TO N DO
FOR J FROM 1 TO N DO
IF I == J AND I < N THEN M9(I,J):=M7(I)+M7(I+1)
END;
IF I == J-1 THEN M9(I,J):=-M7(J)
END;
IF I == J+1 THEN M9(I,J):=-M7(I)
END;
END;
END;
M9(N,N):=M7(N);
// Some symbolic manipulation
temp:=DET(M9-'X'*M8);
M7:=CAS.zeros(temp);
PRINT("Frequencies = "+sqrt(M7)/(2*pi));
END;
```

Edited: 25 Oct 2013, 9:26 a.m.

I can't get it to work as you've written it. When I run an N=3 case, the result should be a 1x3 vector with 3 non-zero values, but it returns a 1x1 with the value zero. I tried CAS.zeros(temp=0) and got a "Bad argument type" error.

Hi! Hope you don't mind if I join here.

When I use "CAS.solve(temp=0)" instead of "CAS.zeros(temp)" it works for me and the Prime returns an 1xN-Vector with N non-zero values.

I still get "Error: Bad argument type"

What version do you have?

Software version 2013 8 13 (5106)

CAS version 1.1.0

Native OS Version V0.025.5106

Edited: 26 Oct 2013, 1:53 p.m.

Quote:
Software version 2013 8 13 (5106)

CAS version 1.1.0

Native OS Version V0.025.5106

How about replace the solve line with
M7:=CAS.zeros(temp, 'X');

I am going to try this on mine (Ver. 5106)

This is seems to work in a test program I did:

```EXPORT TEST1004(A,B,C)
LOCAL temp;
temp:=A*'X^2'+B*'X'+C;
M1:=CAS.zeros(temp,'X');
RETURN M1;
END;
```

TEST1004(1,4,4) returns [-2]
TEST1004(1,0,-4) returns [-2*i, 2*i]

Edited: 26 Oct 2013, 2:02 p.m.

It runs w/o an error, but again gives me Frequencies = [0.00] instead of [4.47,9.38,14.88], which is what I get if I run the solver manually inside the CAS screen.

Are you using CAS.solve or CAS.zeros?

If I am using CAS.solve(temp,'X') I am getting . Using CAS.zeros(temp,'X') I am getting results.

Trying CAS.solve(temp=0,'X') gives me an error.

Edited: 26 Oct 2013, 2:14 p.m.

I did it with CAS.zeros.

We are getting inconsistent results. Is X being single quoted? Does it have to be single quoted?

If anyone else can help, please jump in.

I have CAS.zeros(temp,'X');

I will have to download Frequencies to my calculator and look at it again. So far it is only at my computer (and I am at a local coffee shop right now).

Edited: 26 Oct 2013, 2:40 p.m.

Quote:
Hi! Hope you don't mind if I join here.

When I use "CAS.solve(temp=0)" instead of "CAS.zeros(temp)" it works for me and the Prime returns an 1xN-Vector with N non-zero values.

I am getting an error with CAS.solve(temp=0). And we welcome that you jumped in.

Edited: 26 Oct 2013, 2:20 p.m.

If you do run Frequencies, use the following input values:

N = 3, M1 = 2, M2 = 1, M3 = 1, K1 = 6000, K2 = 4000, K3 = 2000

Also, in the CAS screen I can get it to work with the following:

sqrt(solve(DET(M9-x*M8)=0))/(2*pi)

However, it is necessary to uncheck Exact: in the settings screen.

Strange - now it doesn't work anymore. Without having changed anything...

Downloaded Frequencies onto my calculator (Ver. 5106) as is. I get  like you get.

Trying

```M7:=(CAS.solve(DET(M9-'X'*M8)=0,'X');
```

crashed my calculator and caused a restart.

```M7:=CAS.solve(temp=0,'X');
```

also crashed and caused a restart.

Trying to use CAS.simplify first doesn't do it either. We may be forced to this numerically.

Edited: 28 Oct 2013, 3:53 p.m.

I think Tim W. said that it's not possible to embed non-CAS math operations such as DET inside a CAS operation such as zeros when used in a program, which is why this approach is doomed to failure no matter how much monkeying around we do with temporary variables and quotes. I did manage to expand the matrix determinant into a third order polynomial to get the proper result using zeros and solve, although for some reason it gives me only one root inside a program, whereas it gives me all three in the CAS screen. I then tried the Polynomial root tool "proot" and it returned all the roots inside the program as well as in CAS. Anyways, I'm not sure it's worth going through all these contortions to solve a problem that I've already programmed effectively on my 50g with RPL.

I sounds like proot is the way to go. So far I have only tested it outside the program and I see you have success with it.

Well, here's the program that works fine with proot. Unfortunately, this program is specific to a 3 DOF system only, whereas the original scheme worked with any N DOF system. So, if I want to solve a problem other than N=3, I will have to generate the matrices in the program and then manually run the solver within CAS.

```EXPORT Freq3(K1,K2,K3,M1,M2,M3)
BEGIN
LOCAL C1,C2,C3,C4;
A:=K1;
B:=K2;
C:=K3;
D:=M1;
E:=M2;
F:=M3;
C1:=-D*E*F;
C2:=A*E*F+B*D*F+B*E*F+C*D*E+C*D*F;
C3:=-A*B*F-A*C*E-A*C*F-B*C*D-B*C*E-B*C*F;
C4:=A*B*C;
M0:=CAS.proot({C1,C2,C3,C4});
RETURN sqrt(M0)/(2*pi);
END;
```

Edited: 29 Oct 2013, 12:51 p.m.

Quote:
It runs w/o an error, but again gives me Frequencies = [0.00] instead of [4.47,9.38,14.88], which is what I get if I run the solver manually inside the CAS screen.

This is late in the conversation, but I noticed that you are using DET() without specifying that the determinant be computed in the CAS environment. Try either:

CAS("temp:=DET(M9-'X'*M8)")

or

temp:=CAS.DET(M9-'X'*M8);

What makes it even more nuanced is the fact that you're using a global variable 'X' which has various implications in the home view. Since 'X' is just a dummy variable, why not use 'x' (lower case) which would not have any issues because it is not a built-in variable.

CAS("temp:=DET(M9-'x'*M8)");

Edited: 30 Oct 2013, 11:46 p.m.

Thanks for your detailed input. I found that temp:=CAS.det(M9-'X'*M8) does work, although CAS.solve(temp,'X') only returns one root, whereas running solve in CAS gives me all the roots. Using a lower case "x" causes the calculator to crash and reboot and using upper case letters for "det" results in a syntax error. Also, using CAS("temp:=DET(M9-'X'*M8)") results in a runtime error "Bad argument type".

Update. I was experimenting with using X as a local variable, and when I ran the program it crashed and reset the calculator erasing the CAS history, resetting all the home and CAS settings to the defaults, deleting all my user variables and zeroing all my stored data. Gosh, I just love this thing !

Edited: 31 Oct 2013, 2:23 a.m.

Quote:
Thanks for your detailed input. I found that temp:=CAS.det(M9-'X'*M8) does work, although CAS.solve(temp,'X') only returns one root, whereas running solve in CAS gives me all the roots. Using a lower case "x" causes the calculator to crash and reboot and using upper case letters for "det" results in a syntax error. Also, using CAS("temp:=DET(M9-'X'*M8)") results in a runtime error "Bad argument type".

Update. I was experimenting with using X as a local variable, and when I ran the program it crashed and reset the calculator erasing the CAS history, resetting all the home and CAS settings to the defaults, deleting all my user variables and zeroing all my stored data. Gosh, I just love this thing !

Here's how I got it to work. There is a real big issue with how CAS and Home interact (and in this case, how they seem to fail to interact). Please keep in mind that local variables are not recognized by the current CAS. There is likely a lot of issues with solve() as well. Anyway, because you declared temp and X local variables, they will not work in CAS. I got rid of temp and X from the LOCAL declaration.

`CAS("M7:=zeros(DET(M9-'x'*M8)=0,'x')"); // can change 'x' to 'X' if you want`

Edited: 31 Oct 2013, 9:01 a.m.

Thanks again for all your help. Your detailed explanation in the linked thread makes all of this a lot clearer, although I can see a very long learning curve ahead. It turns out that solve works just fine using this approach and I was having the same problems before with zeros. The reason I wanted to use solve instead of zeros is that solve gives the results in ascending order. The only change I made to your code was to replace the output matrix (vector) M7 with a list L0, which yields the same format as if I ran the problem directly in the CAS and there is no need to use the capabilities of a matrix here. The code is shown below, where the subroutine Spring_Mass() generates the matrices M8 and M9.

```EXPORT Frequencies()
BEGIN
Spring_Mass();
CAS("L0:=solve(DET(M9-'x'*M8)=0,'x')");
RETURN sqrt(L0)/(2*pi);
END;
```

```EXPORT Spring_Mass()
BEGIN
LOCAL I,J,X,N;
INPUT(N,"NDOF","N =","Enter N",0);
M7:=MAKEMAT(0,N);
M8:=MAKEMAT(0,N,N);
M9:=MAKEMAT(0,N,N);
FOR I FROM 1 TO N DO
INPUT(X,"Mass ","Mass "+I+" =","Enter Mass",0);
M8(I,I):=X
END;
FOR I FROM 1 TO N DO
INPUT(X,"Spring ","Spring "+I+" =","Enter Spring",0);
M7(I):=X
END;
FOR I FROM 1 TO N DO
FOR J FROM 1 TO N DO
IF I == J AND I < N THEN M9(I,J):=M7(I)+M7(I+1)
END;
IF I == J-1 THEN M9(I,J):=-M7(J)
END;
IF I == J+1 THEN M9(I,J):=-M7(I)
END;
END;
END;
M9(N,N):=M7(N);
END;
```

Quote:
I can see a very long learning curve ahead.

It has been a pretty steep learning curve for me as well -- and it still is. Maybe one day the interface will be more unified (Home vs CAS) and these issues will no longer be issues.