Algebraic equation in matrix form

Hi, I want to solve a Algebraic equation in matrix form:


Actually, I need to solve AX=0, that is to say I need to solve det(A)=0 with respect to k’ when k in (-1,1).
This is easy in matlab by running:

...
syms k;
kk = double(solve(det(A)));
...

But it takes a few hours to solve the whole problem(need to solve 10^3 same eqs, and in matlab for in for, very slow). So I tried to do it in Julia. However, I have no idea how to do it in Julia.
So anyone knows how to do it?
Any suggestion will be helpful!

1 Like

crosspost on SO: julia - Algebraic equation in matrix form - Stack Overflow

(I believe the policy is that crossposting if fine, but should be declared)

Edit: replaced the wrong link

Do you have to solve it symbolically or would a numerical solution be okay? A numerical solution could likely be obtained in milliseconds.

Hi, a numerical solution is ok. It is true that a numerical solution could be fast. However, I need to solve the eq for every k in (-1,1) and at every moment t in (0,1000) and this will cost a lot of time in matlab. But I am not familiar with Julia, so I dont know how to do it in Julia.

Sorry, I dont understand…

This is a nonlinear eigenproblem; there are much more efficient ways to solve it than black-box root solvers for \det A. See e.g. NonlinearEigenproblems.jl.

Furthermore, if you are trying to find the roots as a function of the continuous parameter k, that is a numerical continuation problem, and there are better algorithms than re-starting the root-finding from scratch for each k. See e.g. BifurcationKit.jl.

1 Like

In fact, it’s even better than that — looking at your matrix, it’s actually a linear eigenproblem in -k'^2 (assuming all of your \eta_\ell factors are constants): if you just move the k'^2 terms to the other side, the equation you wrote takes the form Bc = -k'^2 c, so you can just use sqrt.(.- eigvals(B)) to find the roots k'.

If your matrix is 7 \times 7 as shown in your post, you should be able to solve 1000 of these eigenproblems in a fraction of a second.

3 Likes

Not to mention it’s diagonal plus convolution, so you can use iterative eigensolvers+FFTs+diagonal preconditioning to solve it at a cost linear in the number of unknowns. I imagine this comes from a band structure computation in quantum mechanics? You might be able to adapt GitHub - JuliaMolSim/DFTK.jl: Density-functional toolkit to your purposes (see eg Custom potential · DFTK.jl for how to incorporate a custom potential), although that might be overkill for this relatively simple problem.

2 Likes

What do you mean? This is my question, and I did not post on SO.

Sorry, copypasted the wrong link, edited.

1 Like

Thank you for your kind reply, It really helps a lot! Have a nice day!

Hi, I do tried to solve a band structure. Actually it is a band structure of a photonic crystal of 1d in normal incident case. Thank you for your kind reply, It helps a lot!

Have a nice day!

I tried to solve the problem directly as an algebraic equation and it took a lot of time in matlab. However time cost can be reduced greatly as I solved it in matlab as a eigenproblem as stenengj said.Total time was reduced to 20s-40s in matlab and 8s or so in Julia. Still a lot to learn…