MATLAB Answers

plotting real solution of matrix system - time issue

1 view (last 30 days)
Francesca Lucia Margherita Fabiano
Answered: Dana on 29 Sep 2020
Hi,
I am trying to plot a diagram from the solutions of a system, that is made of real and imaginary solutions.
I am just looking for the plot of the real ones.
Now that I am plotting all the solutions together, it takes several minutes to matlab to plot everything.
Here is the code:
*parameters*
S = [k(1,1)-M*w^2 k(1,2) 0 0;
k(1,2) k(2,2)-Idgear*w^2 0 i*Ipgear*w^2*Omega;
0 0 k(1,1)-M*w^2 k(1,2);
0 -i*Ipgear*w*Omega k(1,2) k(2,2)-Idgear*w^2]
speeds = solve(det(S) ,w)
figure
hold on
ylim([0 2000])
xlim([-1500 1500])
fplot(speeds(1,1))
fplot(speeds(2,1))
fplot(speeds(3,1))
fplot(speeds(4,1))
fplot(speeds(7,1))
fplot(speeds(5,1))
fplot(speeds(6,1))
fplot(speeds(8,1))
xline(360)
ylim([0 2000])
xlim([-1500 1500])

  0 Comments

Sign in to comment.

Answers (1)

Dana
Dana on 29 Sep 2020
I may be misunderstanding things, but it looks to me like the equation you're trying (det(S)=0) to solve is an 8th-order polynomial in w. Algebraic solutions to 8th-order polynomials don't exist except in very special cases. When you apply solve to a case like this, the "solution" isn't typically actually an explicit algebraic solution. Instead, it's an instruction on how to get an explicit numeric solution if necessary. Specifically, the solutions found by solve will look something like:
speeds(1) = root(<polynomial in z>,z,1)
speeds(2) = root(<polynomial in z>,z,2)
.
.
speeds(8) = root(<polynomial in z>,z,8)
where root(<polynomial in z>,z,j) means the j-th root of the specified polynomial in z. If and when MATLAB needs to find an actual explicit solution numerically, it finds these roots numerically.
In your case, I'm guessing the cofficients of the above polynomial in z also contain one (and only one) generic symbolic variable, with everything else numeric. For the sake of the explanation, let's suppose M is the generic symbolic. When you call fplot(speeds(j)), MATLAB first creates a numeric grid for M, and then for each point in that grid, calls a numeric root-finding algorithm to get the desired root.
This is quite inefficient, for (at least) two reasons. First, MATLAB uses variable-precision arithmetic to get numeric solutions to symbolic equations, and that can be quite slow. Unless you need a very high degree of accuracy (e.g., solutions determined to 32 significant digits), you can likely achieve speed improvements by solving using regular double-precision arithmetic.
Second, there's some duplication of effort here. For example, consider the first gridpoint of M, call it m1. When you do fplot(speeds(1)), MATLAB first sets M=m1 in the polynomial and then finds the 1st root. When you do fplot(speeds(2)), MATLAB does the same thing except finds the 2nd root, fplot(speeds(3)) the 3rd root, etc. This means you are finding each of the roots of the polynomial with M=m1 separately. Since nearly all of the work involved in finding all the roots of the polynomial with M=m1 will already have been done when you call fplot(speeds(1)), there are wasted calculations when you then do fplot(speeds(2)), and then fplot(speeds(3)), etc.
Instead of how you've done it, I would do something like the following (again, for the sake of argument assuming M is the only non-numeric variable in S other than w):
*parameters*
S = [k(1,1)-M*w^2 k(1,2) 0 0;
k(1,2) k(2,2)-Idgear*w^2 0 i*Ipgear*w^2*Omega;
0 0 k(1,1)-M*w^2 k(1,2);
0 -i*Ipgear*w*Omega k(1,2) k(2,2)-Idgear*w^2]
% We'll think of detS as a polynomial in w
detS = simplify(det(S));
% Get vector of coefficients of the polynomial in w, and flip them
% so that the first element corresponds to the highest degree of w.
c = fliplr(coeffs(detS,w));
Mbds = [-1500 1500]; % bounds of M range
nM = 100; % number of values of M to get roots at
Mvec = linspace(Mbds(1),Mbds(2),nM)'; % vector of M values
rts = nan(nM,8); % create matrix to hold roots
for j = 1:nM % for each M gridpoint
% substitute in for M, then convert to numeric (double) vector
cj = double(subs(c,M,Mvec(j)));
rtj = roots(cj); % get all 8 roots of polynomial
rtj = rtj(imag(rtj)==0); % drop complex roots
rtj = sort(rtj); % sort remaining roots
nrtj = numel(rtj); % how many real roots there are
% Assign result to roots matrix. Note that if there were complex
% roots, some elements of rts(j,:) won't be assigned and will
% stay as NaN. So if you see an NaN in the end, you'll know that
% entry is not actually a root.
rts(j,1:nrtj) = rtj';
end
figure
hold on
plot(Mvec,rts)
xline(360)
ylim([0 2000])
xlim([-1500 1500])

  0 Comments

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!