Speeding up computation using fsolve

Hi,

I am using fsolve to solve a transcendental eigenvalue problem and since I need to solve it for multiple frequencies, I am trying to find a way to speed up the calculation. I have no experience with optimization and therefore I'm not sure where to look for answers.

The code I am trying to optimize is :

omegas=linspace(100*2*pi,1e4*2*pi);
kguess=[1e-4 5e-4 1e-3 5e-3 1e-2 5e-2 0.1 0.5 1:20 100 200 500];

for iW=1:length(omegas)
   w=omegas(iW);
    for ik=1:length(kguess)
        k0=kguess(ik);
        f = @(k)fun_transcendental(k,w,K,M,T_per,Lx,Ly,theta);
       options=optimset('MaxFunEvals',1000,'MaxIter',1000,'Display','off');
       x0(ik) = fsolve(f,k0,options);
     end
  %post-processing
  X(iW,:)=x0_postprocessed;
end

where the function fun_transcendental is solving det(Kr-w^Mr)=0 where Kr and Mr are functions of exp(-i*cst*k). My goal is to solve for the complex eigenvalues k(w).

A couple comments:

1) I keep on getting the following error although it does find the correct solution. I am not sure why.

No solution found.
fsolve stopped because the relative size of the current step is less than the default value of the step size tolerance squared, but the vector of function values is not near zero as measured by the default value of the function tolerance.

2) Obviously my kguess is quite large and I get repeated solutions. I post-process the vector x0 to only keep the physical eigenvalues (i.e., eigenvalues that return an eigenvector when using the null function), and then only keep of of the repeated eigenvalues. In the end, for each w, I should have 3 eigenvalues. I have tried using the 3 eigenvalues of the solution (the real part and their absolute value) at iW-1 as initial guess for the iW loop but it failed (it only found 2/3 solutions for the second loop and 1/3 for the third).

3) I know that I could use a parallel for-loop for the w loop. I actually have tried and it speeds up things quite a lot. However, I might not be able to use the parallel toolbox, so I'd really like to understand if and how I can speed up this code using the for-loop only.

Thanks.

[Walter Roberson]

Could you post fun_transcendental ?

[Clemence]

Sure.

function F=fun_transcendental(k,w,K,M,T_per,Lx,Ly,theta)

lx=exp(-1i*k*Lx*cos(theta));
ly=exp(-1i*k*Ly*sin(theta));

T=T_per.t0+lx*T_per.tx+ly*T_per.ty+lx*ly*T_per.txy;
Kr=T'*K*T;
Mr=T'*M*T;

Fa=(Kr-w.^2*Mr);
F=det(Fa);

The T_per.ti matrices enforce periodicity and only contain 1s and 0s..

[Walter Roberson]

It appears we need values for M, T_per, Lx, Ly, theta to execute this.

[Clemence]

Could I email them to you directly? K and M are 45x45 matrices.

[Walter Roberson]

You can attach a .mat file

[Clemence]

Sorry, I didn't know you could do that. Thanks for your help.

[Walter Roberson]

Your posted code writes the fsolve() answers into x0, but then stores the undefined variable x0_postprocessed ?

[Walter Roberson]

You are giving real-valued guesses, but all of your solutions are complex-valued ?

[Walter Roberson]

I had hoped that I might be able to get somewhere by using a symbolic k and doing a symbolic solve for the determinant. Unfortunately the matrices involved are too large (40x40) for that to be practical, even taking into account that they have a fair number of 0s in them.

I am not clear as to whether you are expecting to be able to find a real-valued k for the solution?

The last few versions of MATLAB are happy to explore imaginary solutions in fsolve if you give them only a single guess; earlier versions would get lost and could not possibly find a solution once even a single complex result occurred.

If you are looking for real-valued k then you should use a two-element guess instead of a single element; that confines the search range. Unfortunately though when you do that, there is no way to say, "but within that range, try looking near _ first"

[Clemence]

Thanks Walter. With the matrices I gave, I am expecting mostly real-valued solutions, although some of them are complex with small imaginary parts. For now, I am only interested in the real part of the solution but in the future, I'll be looking at complex solutions. So I'd like my code to be as generic as possible.

By displaying the iteration details, I noticed that the matrix Kr-w^2Mr is ill-conditioned, leading to infinite determinant instead of zero. I am not sure if I could try to do something on that side so the solution converges faster, or if it doesn't really matter.