Vectorizing a function looking for eigenvalues of a matrix

4 次查看（过去 30 天）
Ferdy Ellen
Ferdy Ellen 2022-10-25
评论： Walter Roberson 2022-10-25
In short, I have a function that defines a matrix based on an input value, and then returns the eigenvalue of that matrix that is closest to 1. I want to plot this for different inputs, but I get the following warning:
'Warning: Function behaves unexpectedly on array inputs. To improve performance, properly vectorize your function to return an output with the same size and shape as the input arguments.'
Is it possible to vectorize this function in order to plot it properly (and without waiting a very long time)?
My code looks somewhat like this (if necessary I can post the entire code but it's pretty long).
eigenvalue = @eigenvalue_func
fplot(eigenvalue)
function eta = eigenvalue_func(z)
p = 1:1:100;
q = p;
weights = rand(100,1); %in reality I'm using a distribution for p, q and the weights based on Gaussian Quadrature
[Q, P] = meshgrid(q,p);
[W, ~] = meshgrid(weights, p);
v_0 = -30./(P.*Q.*(P.^2-Q.^2).*2.*pi.^2).*(Q.*cos(Q).*sin(P)-P.*cos(P).*sin(Q));
A = -4.*pi.*v_0.*(Q.^2 - z ).^(-1).*Q.^2.*W;
for j = 1:100 %the diagonal has to be defined differently since the above returns NaN/Inf (I take the limit here)
v_0(j,j) = -30./(2.*pi^2.*P(j,j).^2).*(0.5-(sin(2.*P(j,j)))./(4.*P(j,j)));
A(j,j) = -4.*pi.*v_0(j,j).*(Q(j,j).^2 - z).^(-1).*Q(j,j).^2.*W(j,j);
end
eta_full = eigs(A, 20);
[~, index_eigenvalue] = min(abs(eigs(A, 20)-1));
eta = eta_full(index_eigenvalue);
end
  6 个评论
显示 4更早的评论
Bruno Luong
Bruno Luong 2022-10-25
编辑：Bruno Luong 2022-10-25
The bottleneg is probably calling EIGS (have you made any profiler of your code?) which cannot be vectorized at the moment, and possibly fplot waste a lot of time to find the interval and resolution to plot the function.
If you want faster plot, change fplot and eigs.
For instant it is totally iefficient to call twice eigs with 20 eigenvalues
eta_full = eigs(A, 20);
[~, index_eigenvalue] = min(abs(eigs(A, 20)-1));
eta = eta_full(index_eigenvalue)
where you can get the same result with
eta = eigs(A,1,1);
Walter Roberson
Walter Roberson 2022-10-25
Since eigs is not internally vectorized, it might make sense to define a function that accepted a vector of z values, and used background pool to evaluate multiple z in parallel.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

Walter Roberson
Walter Roberson 2022-10-25
wrapper = @(z) arrayfun(eigenvalue, z)
fplot(wrapper)

类别

在 Help Center 和 File Exchange 中查找有关 Linear Algebra 的更多信息

标签

产品


版本

R2019b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by