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!
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.
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.
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.
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.
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!
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…