Utilisation des racines et résolution dans MATLAB

Question

J'ai une équation qui va comme ceci :

Ici, a, b, c, d, e, f, g, h sont des constantes. Je veux faire varier k de 0,6 à 10 avec un intervalle de 0,1 et trouver w. Il y a deux façons de le faire dans MATLAB.

Une façon de procéder est de convertir cette équation en une équation de la forme (quelque chose)w^8-(quelque chose d'autre)w^6....-(quelque chose d'autre encore)w^0=0, puis d'utiliser la commande 'roots' dans MATLAB (Méthode 1).

Une autre méthode consiste à définir des fonctions symboliques, puis à exécuter le programme. Lorsque vous utilisez cette méthode, vous n'avez pas besoin de simplifier davantage l'expression, vous pouvez simplement la mettre sous sa forme initiale (méthode 2).

Les deux façons sont montrées dans le script ci-dessous :

%%% 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

La méthode 1 utilise la commande roots et la méthode 2 utilise la commande symbolique et solve. Mes questions sont les suivantes :

1. Vous pouvez voir que les intrigues de la première section arrivent en peu de temps, alors que la seconde prend plus de temps. Y a-t-il un moyen d'augmenter la vitesse ?
2. Les tracés des deux sections semblent très différents, et vous pouvez être forcé de croire que j'ai fait des erreurs en effectuant le calcul à partir de (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)) à (quelque chose)w^8-(autre chose)w^6....-(autre chose encore)w^0 . Je peux vous assurer que je l'ai formulé correctement. Vous pouvez voir ce qui se passe réellement, si vous recherchez une valeur particulière de krho dans les deux fichiers de données (stack_using_roots et stack_using_solve). Pour, disons, krho=3.6, les racines sont les mêmes dans les deux fichiers de données, mais la façon dont elles sont "écrites" n'est pas correcte. C'est pourquoi les tracés sont maladroits. En bref, en utilisant la commande 'roots', les solutions sont données dans un format ordonné, d'autre part en utilisant 'solve', elles sont décalées de façon aléatoire. Que se passe-t-il réellement ? Existe-t-il un moyen de contourner ce problème ?

3. J'ai lancé le programme avec

   i) symétriquement à w avec n(i, :)=double(solve(first-second-third==0, w)) ;

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

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

Dans ces 3 cas, les résultats semblent être les mêmes. Alors quelle est la chose que nous devons définir comme symbolique ? Et quand utilise-t-on et n'utilise-t-on pas l'expression '==0' ?

Answer 1

1. Y a-t-il des moyens d'augmenter la vitesse ?

Plusieurs. Certaines améliorations triviales de la vitesse viendraient de la définition des variables avant la boucle. Le gros goulot d'étranglement est solve . Malheureusement, il n'y a pas de solution analytique évidente à votre problème sans savoir k à l'avance, il n'y a donc pas de moyen évident de tirer solve en dehors de la boucle for.

2. En bref, en utilisant la commande 'roots', les solutions sont données dans un format ordonné, alors qu'en utilisant 'solve', elles sont déplacées de façon aléatoire. Comment cela se fait-il ?

Ce n'est pas vraiment être "décalé". Votre fonction est symétrique autour de w = 0. Donc, pour chaque racine r, il y a une autre racine à -r. Chaque fois que vous appelez solve, cela vous donne la première, la deuxième, la troisième puis la quatrième racine, et ensuite la même chose mais cette fois les racines sont multiplié par -1.

Parfois, solve choisit de retirer -1 comme facteur commun. Dans ces cas, il vous donne d'abord les racines multipliées par -1, puis les racines positives. Je ne sais pas pourquoi il enlève parfois -1 et parfois non, mais dans votre cas (puisque vous ne vous souciez pas de la partie imaginaire), vous pouvez résoudre ce problème en remplaçant double(solve(first-second-third, w)) avec sort(real(double(solve(first-second-third, w)))) . L'ordre des racines ne sera pas le même que dans la méthode 1, mais vous n'aurez pas le comportement bizarre de commutation.

3. Dans ces 3 cas, les résultats semblent être les mêmes. Alors quelle est la chose que nous devons définir comme symbolique ? Et quand utilise-t-on et n'utilise-t-on pas l'expression '==0' ?

syms w k vs. syms w ne fait pas de différence car vous redéfinissez k comme une valeur numérique (0.6, 0.7,... etc). Seul w doit être symbolique.

Le site page de référence para solve dicte comment l'équation doit être spécifiée. Si vous faites défiler vers le bas jusqu'à la section concernant la variable d'entrée eqns Elle stipule

Si des éléments de eqns sont des expressions symboliques (sans le côté droit), solve égalise l'élément à 0.

C'est pourquoi cela ne fait aucune différence que vous écriviez first-second-third==0 o first-second-third comme première entrée à solve .