Solving for real and imaginary eigenvalues

[Eric Kangas]

Hi,

First of all I am using a coupled 2sd ODE.

First step was to declare variables:

a,l,x,y,u,v,xdot,ydot,udot,vdot = var('a,l,x,y,u,v,xdot,ydot,udot,vdot', domain = RR)

Second step turn the two 2sd ODE into four 1st ODE.

xdot = u;

ydot = v;

udot = -(1-a)*((x-a)/(sqrt((x-a)^2 + y^2))^3)-a*((x+1-a)/(sqrt((x+1-a)^2 + y^2))^3)+x+2*v;

vdot = -(1-a)*(y/(sqrt((x-a)^2 + y^2))^3)-a*(y/(sqrt((x+1-a)^2 + y^2))^3)+y-2*u;

Third I had to process the Jacobian by long hand:

N = matrix([(diff(xdot,x)-l,diff(xdot,y),diff(xdot,u),diff(xdot,v)),(diff(ydot,x),diff(ydot,y)-l,diff(ydot,u),diff(ydot,v)),(diff(udot,x),diff(udot,y),diff(udot,u)-l,diff(udot,v)),(diff(vdot,x),diff(vdot,y),diff(vdot,u),diff(vdot,v)-l)]);

Fourth, I took the determinate by long hand. Sage wouldn't allow me to process it by using the built in determinate function:

dN = N[0,0]*(N[1,1]*(N[2,2]*N[3,3]-N[2,3]*N[3,2]) + N[1,2]*(N[1,2]*N[3,3]-N[3,2]*N[1,3]) + N[1,3]*(N[1,2]*N[2,3]-N[2,2]*N[1,3])) + N[1,0]*(N[0,1]*(N[2,2]*N[3,3]-N[3,2]*N[2,3]) + N[2,1]*(N[0,2]*N[3,3]-N[3,2]*N[0,3]) + N[3,1]*(N[0,2]*N[2,3]-N[2,2]*N[0,3])) + N[2,0]*(N[1,0]*(N[1,2]*N[3,3]-N[3,2]*N[1,3]) + N[1,1]*(N[0,2]*N[3,3]-N[3,2]*N[0,3]) + N[1,3]*(N[0,2]*N[1,3]-N[1,2]*N[0,3])) + N[3,0]*(N[0,1]*(N[1,2]*N[2,3]-N[2,2]*N[1,3]) + N[1,1]*(N[0,2]*N[2,3]-N[2,2]*N[0,3]) + N[2,1]*(N[0,2]*N[1,3]-N[1,2]*N[0,3]));

Fifth I had to expand the determinate out to get the coefficients, and to create the determinate in terms of the eigenvalue to solve for the eigenvalue:

dN1 = expand(dN);

dN2 = dN1.coeffs(l); f1 = (dN2[1][0])*l^(dN2[1][1])+(dN2[2][0])*l^(dN2[2][1])+(dN2[3][0])*l^(dN2[3][1])

f2 = dN2[0][0];

soln = solve([f1 == 0], l);

Last I then created the solutions which turned out to be four solutions:

rsol0(a,x,y) = sol0(a,x,y).real()

isol0(a,x,y) = sol0(a,x,y).imag()

rsol1(a,x,y) = sol1(a,x,y).real()

isol1(a,x,y) = sol1(a,x,y).imag()

rsol2(a,x,y) = sol2(a,x,y).real()

isol2(a,x,y) = sol2(a,x,y).imag()

rsol3(a,x,y) = sol3(a,x,y).real()

isol3(a,x,y) = sol3(a,x,y).imag()

However, the computer hangs, and doesn't even finish the computation. Unless I assign values to either x,y, and/or a. Do I have to first assign values then take the real, and imaginary parts? If so that will end up being an issue with what I want to do with these eigenvalues.

[Eric Kangas]

Missed a step between the last and step before last step:

sol0(a,x,y) = soln[0].rhs()+f2;

sol1(a,x,y) = soln[1].rhs()+f2;

sol2(a,x,y) = soln[2].rhs()+f2;

sol3(a,x,y) = soln[3].rhs()+f2;

[kcrisman]

On Tuesday, June 26, 2012 7:11:50 AM UTC-4, Eric Kangas wrote:

Hi,

First of all I am using a coupled 2sd ODE.

First step was to declare variables:

a,l,x,y,u,v,xdot,ydot,udot,vdot = var('a,l,x,y,u,v,xdot,ydot,udot,vdot', domain = RR)

Second step turn the two 2sd ODE into four 1st ODE.

xdot = u;

ydot = v;

udot = -(1-a)*((x-a)/(sqrt((x-a)^2 + y^2))^3)-a*((x+1-a)/(sqrt((x+1-a)^2 + y^2))^3)+x+2*v;

vdot = -(1-a)*(y/(sqrt((x-a)^2 + y^2))^3)-a*(y/(sqrt((x+1-a)^2 + y^2))^3)+y-2*u;

At this point you should be able to do

V(x,y,u,v) = [xdot,ydot,udot,vdot]

V.diff()

to get the Jacobian matrix.  It's true that we don't implement eigenvalues for such generic symbolic matrices - presumably in general the computation could take arbitrarily long or be too difficult for any system?

Also, I don't know why the determinant didn't work for you.

expand((V.diff()-diagonal_matrix([l,l,l,l])).det()).coeffs(l)

gets you go dN2, though for whatever reason I get only three things there instead of your four.

Can you do this same process with a less complex system to show what *does* work?  The coefficients here are ridiculously huge and so I'm not sure what's happening.  At any rate trying this with some random expressions (didn't do the solving) doesn't give this hang.

[Eric Kangas]

Interesting that there are only 3 coefficients insterad of 4 by going through the process you had mentioned. I originally started this project on mathematica 6.0 student edition back in college, and I now need to recreate it in sage so I can publish a paper.

[kcrisman]

On Monday, July 30, 2012 6:37:00 PM UTC-4, Eric Kangas wrote:

Interesting that there are only 3 coefficients insterad of 4 by going through the process you had mentioned. I originally started this project on mathematica 6.0 student edition back in college, and I now need to recreate it in sage so I can publish a paper.

Do you still get the hang this way?  I might expect one with "solve", since Maxima could take very long to solve some relations, but no idea on the .real() bit.

[Eric Kangas]

There is a real and imaginary solution for each eigenvalue. So a total of 8 eigenvalues from when I did it on mathematica.

[kcrisman]

On Wednesday, August 1, 2012 10:08:56 AM UTC-4, Eric Kangas wrote:

There is a real and imaginary solution for each eigenvalue. So a total of 8 eigenvalues from when I did it on mathematica.

Okay, but do you get the hang?  That was the original question.

I'm not going to be able to do much Sage support the next couple weeks, unfortunately, so hopefully someone else will pick up this thread :)