Usage of roots and solve in MATLAB

Question

I have an equation which goes like this:

Here a, b, c, d, e, f, g, h are constants. I want to vary k from say .6 to 10 with 0.1 interval and to find w. There are 2 ways to do this in MATLAB.

One way is convert this equation into equation of the form (something)w^8-(something else)w^6....-(something else again)w^0=0, and then make use of command 'roots' in MATLAB (Method 1).

Another way is that defining symbolic functions and then executing the program. When you are using the this method, you may not need simplify the expression any further, you can just put it in first form itself (Method 2).

Both ways are shown in the script below:

%%% defining values
clear; clc;
a=0.1500;
b=0.20;
c=0.52;
d=0.5;
e=6;
f=30;
g=18;
h=2;

%% Method 1: varying k using roots
tic
i=0;
for k=.6:.1:10
    i=i+1;
    t8=a;
    t7=0;
    t6=-(1+e+a*(c+g))*(k^2) ;
    t5=0;
    t4=(k^2*(b+f+(c*e+g)*k^2)-a*(d+h-c*g*k^4));
    t3=0;
    t2=k^2*(d*(e+a*g)+h+a*c*h-(c*f+b*g)*k^2);
    t1=0;
    t0=(a*d*h)-(d*f+b*h)*k^2;
    q=[t8 t7 t6 t5 t4 t3 t2 t1 t0];
    r(i,:)=roots(q);
end
krho1(:,1)=.6:.1:10;
r_real=real(r);
r_img=imag(r);
dat1=[krho1 r_real(:,1) r_real(:,2) r_real(:,3) r_real(:,4) r_real(:,5) r_real(:,6) r_real(:,7) r_real(:,8)];
fnameout=('stack_using_roots.dat');
fid1=fopen(fnameout,'w+');
fprintf(fid1,'krho\t RR1\t RR2\t RR3\t RR4\t RR5\t RR6\t RR7\t RR8\t \r');
fprintf(fid1,'%6.4f %7.10f %7.10f %7.10f %7.10f %7.10f %7.10f %7.10f %7.10f \n',dat1');
fclose(fid1);
plot(krho1, r_real(:,1),krho1, r_real(:,2),krho1, r_real(:,3),krho1, r_real(:,4),krho1, r_real(:,5),krho1, r_real(:,6),krho1, r_real(:,7),krho1, r_real(:,8))
toc

%% Method 2: varying k using solve
tic
syms w k
i=0;
for k=.6:.1:10
    i=i+1;
    first=a/k^2;
    second=(w^2-b)/(w^4-k^2*c*w^2-d) ;
    third=(e*w^2-f)/(w^4-k^2*g*w^2-h);
    n(i,:)=double(solve(first-second-third, w));
end
krho1(:,1)=.6:.1:10;
r_real=real(n);
r_img=imag(n);

dat1=[krho1 r_real(:,1) r_real(:,2) r_real(:,3) r_real(:,4) r_real(:,5) r_real(:,6) r_real(:,7) r_real(:,8)];
fnameout=('stack_using_solve.dat');
fid1=fopen(fnameout,'w+');
fprintf(fid1,'krho\t RR1\t RR2\t RR3\t RR4\t RR5\t RR6\t RR7\t RR8\t \r');
fprintf(fid1,'%6.4f %7.10f %7.10f %7.10f %7.10f %7.10f %7.10f %7.10f %7.10f \n',dat1');
fclose(fid1);
figure;
plot(krho1, r_real(:,1),krho1, r_real(:,2),krho1, r_real(:,3),krho1, r_real(:,4),krho1, r_real(:,5),krho1, r_real(:,6),krho1, r_real(:,7),krho1, r_real(:,8))
toc

Method 1 uses roots command and Method 2 uses symbolic and solve. My question are as follows:

1. You can see that first section plots are coming in a short time, where the second one is taking a greater time. Is there any way to increase the speed?
2. The plots of both the section seems very different, and you may be forced to believe that I have made mistakes while carrying out the calculation from (a/k^2)-((w^2-b)/(w^4-k^2*cw^2-d))-((ew^2-f)/(w^4-k^2*g*w^2-h)) to (something)w^8-(something else)w^6....-(something else again)w^0. I can assure you that I have put it correctly. You can see what really happens, if you look for any particular value of krho in both the dat file (stack_using_roots and stack_using_solve). For, lets say, krho=3.6, the roots is the same in both the dat files, but the way in which it is 'written' is not in a proper way. That is why the plots looks awkward. In short, while using 'roots' command, the solutions are given in a orderd format, on the other hand while using 'solve', it is getting shifted randomly. What is really happening? Is there any way to get around this problem?

3. I have ran the program with

   i) syms w along with n(i,:)=double(solve(first-second-third==0, w));

   ii) syms w k along with n(i,:)=double(solve(first-second-third==0, w));

   iii) syms w k along with n(i,:)=double(solve(first-second-third, w));

In all these 3 cases, results seem to be same. Then what is the thing that we have to define as symbolic? And when do we use and do not use the expression '==0'?

Answer 1

1. Are there any ways to increase the speed?

Several. Some trivial speed improvements would come from defining variables before the loop. The big bottleneck is solve. Unfortunately, there isn't an obvious analytical solution to your problem without knowing k beforehand, so there's no obvious way to pull solve outside the for loop.

2. In short, while using 'roots' command, the solutions are given in a ordered format, on the other hand while using 'solve', it is getting shifted randomly. Why is that?

It is not really getting "shifted". Your function is symmetric about w = 0. So, for every root r there is another root at -r. Every time you call solve, it gives you the first, second, third then fourth roots, and then the same thing but this time the roots are multiplied by -1.

Sometimes solve chooses to take out -1 as a common factor. In these cases, it first gives you the roots multiplied by -1, then the positive roots. Why it sometimes takes out -1, sometimes doesn't, I don't know, but in your case (since you don't care about the imaginary part) you can fix this by replacing double(solve(first-second-third, w)) with sort(real(double(solve(first-second-third, w)))). The order of the roots won't be the same as in Method 1, but you won't get the weird switching behaviour.

3. In all these 3 cases, results seem to be same. Then what is the thing that we have to define as symbolic? And when do we use and do not use the expression '==0'?

syms w k vs. syms w doesn't make a difference because you redefine k as a numeric value (0.6, 0.7,... etc). Only w needs to be symbolic.

The reference page for solve dictates how the equation should be specified. If you scroll down to the section regarding the input variable eqns, it states

If any elements of eqns are symbolic expressions (without the right side), solve equates the element to 0.

This is why it makes no difference whether you write first-second-third==0 or first-second-third as the first input to solve.