Usage of roots and solve in MATLAB

I have an equation which goes like this:

(a/k^2)-((w^2-b)/(w^4-k^2*c*w^2-d))-((e*w^2-f)/(w^4-k^2*g*w^2-h)) (I don’t know how to put this in equation style).

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.

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.

I am giving the programs to do that in both ways:

    %%defining values
    clear all; clc;
    a=0.1500;
    b=0.20;
    c=0.52;
    d=0.5;
    e=6;
    f=30;
    g=18;
    h=2;
    %%for 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
    %%for 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);
    % 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

You can comment second section, run the program, and then afterwards comment first section and run the program.

First section (%% for varying k using roots) uses roots command and second section (%% for varying k using solve) 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*c*w^2-d))-((e*w^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

• syms w along with n(i,:)=double(solve(first-second-third==0, w));
• syms w k along with n(i,:)=double(solve(first-second-third==0, w));
• 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

[Karan Gill]

1. To understand the difference between numeric and symbolic arithmetic, please see: Choose Symbolic or Numeric Arithmetic

2. The ordering of roots depends on the internal algorithm. There's no way to order them.

3. From the "solve" doc page: " If eqn is a symbolic expression (without the right side), the solver assumes that the right side is 0, and solves the equation eqn == 0. "