Checking if my program calculates what is supposed to.

Hello all, how are you?

I am trying to write a program that calculates aeroelasticiy stuff. I am not sure if the program is calculating what is supposed to.

The plan is:

1- Guess a value of k between 0 and 1

2- Calculate Ck

3- Calculate f11

4- Calculate f12

5- Calculate f21

6- Calculate f22

7- Calculate the polynomial roots of (p)

8- Pick one of the roots of the polynomial p (let´s say for example the first one)

9- The new value of k is equal to the real part of the root of the polynomial p (first root)

10- Return to step 2 and calculate everthing until k gets convergence.

>> %inputs

a = -1/5;

mu = 2; % 2 ≤ µ ≤ 20

rs = 1/4;

sigma = 1/sqrt(10);

xtheta = 0; % xtheta = 0.0 and 0.2

V = 0.05;

%equations and pre-defining vectors/matrices

k = 0:0.01:1;

r1 = zeros(1,101) ;

r2 = zeros(1,101) ;

r3 = zeros(1,101) ;

r4 = zeros(1,101) ;

for j = 2 : 101;

k(j);

Ck = (0.01365+0.2808.*k.*1i-k.^2./2) / (0.01365+0.3455.*k.*1i-k.^2);

f11 = sigma^2./V.^2 - k.^2./mu + (2*1i.*Ck)/mu;

f12 = (k.*(1i + a.*k) + (2 + 1i.*k.*(1-2*a)).*Ck)/mu;

f21 = (a.*k.^2 - 1i.*k.*(1 + 2*a).*Ck)/mu;

f22 = ((8*mu*rs)./V^2 + 4*1i*(1 + 2*a)*(2*1i - k.*(1 - 2*a)).*Ck - k.*(k - 4*1i + 8*a*(1i + a.*k)))/8*mu;

p = [(rs-xtheta^2) 0 (f22 + f11.*rs - f21.*xtheta - f12.*xtheta) 0 (f11.*f22 - f12.*f21)];

rts = roots(p);

r1(j) = rts(1);

r2(j) = rts(2);

r3(j) = rts(3);

r4(j) = rts(4);

k = imag(r1)

gamma = real (r1)./k;

end

plot(k,gamma);

hold on;

ylabel('γ');

xlabel('k');

hold off;

is this correct ?

Answer 1

[Mathieu NOE]

hello

yes your code needed some fix

at least I believe I have fixed most recurrence issues in the loop

now there are other issues , like

you said : 9- The new value of k is equal to the real part of the root of the polynomial p (first root)

in your code you're taking the iamginary part , not the real part

if I go with real, k tends to converge to zero whatever the initial value

if I go with imag, k bounces between two values

I guess none of these behaviour is ok

% >> %inputs
a = -1/5;
mu = 2; % 2 ≤ µ ≤ 20
rs = 1/4;
sigma = 1/sqrt(10);
xtheta = 0; % xtheta = 0.0 and 0.2
V = 0.05;
%equations and pre-defining vectors/matrices
k = 0.75;
k_visu(1) = k;

n = 10;% nb of loops
% r1 = zeros(1,n) ; 
% r2 = zeros(1,n) ;
% r3 = zeros(1,n) ;
% r4 = zeros(1,n) ;
for j = 2 : n;
    Ck = (0.01365+0.2808.*k.*1i-k.^2./2) / (0.01365+0.3455.*k.*1i-k.^2);
    f11 = sigma^2./V.^2 - k.^2./mu + (2*1i.*Ck)/mu;
    f12 = (k.*(1i + a.*k) + (2 + 1i.*k.*(1-2*a)).*Ck)/mu;
    f21 = (a.*k.^2 - 1i.*k.*(1 + 2*a).*Ck)/mu;
    f22 = ((8*mu*rs)./V^2 + 4*1i*(1 + 2*a)*(2*1i - k.*(1 - 2*a)).*Ck - k.*(k - 4*1i + 8*a*(1i + a.*k)))/8*mu;
    p = [(rs-xtheta^2) 0 (f22 + f11.*rs - f21.*xtheta - f12.*xtheta) 0 (f11.*f22 - f12.*f21)];
    rts = roots(p); 
%     r1(j) = rts(1);
%     r2(j) = rts(2);
%     r3(j) = rts(3);
%     r4(j) = rts(4);
    k = real(rts(1));
    k_visu(j) = k;
%     gamma(j) = real (rts(1))./k;
end

plot(1:n,k_visu);

[Rodrigo Pena]

Thank you so much !

Maybe the value converges to zero for the first root.

for the second root, the value may converge to another number.

[Rodrigo Pena]

Hello, I have got an idea:

I would like to get the number in red and throw him as my new value of k. how to do this ?

[Mathieu NOE]

hello seems tyhe idea is working

now we have k converging towards : -0.0688

code :

% >> %inputs
a = -1/5;
mu = 2; % 2 ≤ µ ≤ 20
rs = 1/4;
sigma = 1/sqrt(10);
xtheta = 0; % xtheta = 0.0 and 0.2
V = 0.05;
%equations and pre-defining vectors/matrices
k = 0.5;
k_visu(1) = k;
n = 10;% nb of loops

for j = 2 : n;
    Ck = (0.01365+0.2808.*k.*1i-k.^2./2) / (0.01365+0.3455.*k.*1i-k.^2);
    f11 = sigma^2./V.^2 - k.^2./mu + (2*1i.*Ck)/mu;
    f12 = (k.*(1i + a.*k) + (2 + 1i.*k.*(1-2*a)).*Ck)/mu;
    f21 = (a.*k.^2 - 1i.*k.*(1 + 2*a).*Ck)/mu;
    f22 = ((8*mu*rs)./V^2 + 4*1i*(1 + 2*a)*(2*1i - k.*(1 - 2*a)).*Ck - k.*(k - 4*1i + 8*a*(1i + a.*k)))/8*mu;
    p = [(rs-xtheta^2) 0 (f22 + f11.*rs - f21.*xtheta - f12.*xtheta) 0 (f11.*f22 - f12.*f21)];
    rts = roots(p); 
    rts_real_sorted = sort(real(rts));
    k = rts_real_sorted(1);
    k_visu(j) = k;
%     gamma(j) = real (rts(1))./k;
end

plot(1:n,k_visu);

[Rodrigo Pena]

Yaay ! thank you !

I have got the same results as well.

I expanded the same theory to the other roots.

The challenge now is store these converged values.

%inputs

a = -1/5;

mu = 2; % 2 ≤ µ ≤ 20

rs = 1/4;

sigma = 1/sqrt(10);

xtheta = 0; % xtheta = 0.0 and 0.2

V = 0.05;

%equations and pre-defining vectors/matrices

k = 0.1;

k1_visu(1) = k

n = 10; % nb of loops

for j = 2 : n;

Ck = (0.01365+0.2808.*k.*1i-k.^2./2) / (0.01365+0.3455.*k.*1i-k.^2);

f11 = sigma^2./V.^2 - k.^2./mu + (2*1i.*Ck)/mu;

f12 = (k.*(1i + a.*k) + (2 + 1i.*k.*(1-2*a)).*Ck)/mu;

f21 = (a.*k.^2 - 1i.*k.*(1 + 2*a).*Ck)/mu;

f22 = ((8*mu*rs)./V^2 + 4*1i*(1 + 2*a)*(2*1i - k.*(1 - 2*a)).*Ck - k.*(k - 4*1i + 8*a*(1i + a.*k)))/8*mu;

p = [(rs-xtheta^2) 0 (f22 + f11.*rs - f21.*xtheta - f12.*xtheta) 0 (f11.*f22 - f12.*f21)];

rts = roots(p)

K = sort(rts, 'ComparisonMethod','real')

k1 = real(K(1))

k1_visu(j) = k1

end

k = 0.1;

k2_visu(1) = k

for j = 2 : n;

Ck = (0.01365+0.2808.*k.*1i-k.^2./2) / (0.01365+0.3455.*k.*1i-k.^2);

f11 = sigma^2./V.^2 - k.^2./mu + (2*1i.*Ck)/mu;

f12 = (k.*(1i + a.*k) + (2 + 1i.*k.*(1-2*a)).*Ck)/mu;

f21 = (a.*k.^2 - 1i.*k.*(1 + 2*a).*Ck)/mu;

f22 = ((8*mu*rs)./V^2 + 4*1i*(1 + 2*a)*(2*1i - k.*(1 - 2*a)).*Ck - k.*(k - 4*1i + 8*a*(1i + a.*k)))/8*mu;

p = [(rs-xtheta^2) 0 (f22 + f11.*rs - f21.*xtheta - f12.*xtheta) 0 (f11.*f22 - f12.*f21)];

rts = roots(p)

K = sort(rts, 'ComparisonMethod','real')

k2 = real(K(2))

k2_visu(j) = k2

end

k = 0.1;

k3_visu(1) = k

for j = 2 : n;

Ck = (0.01365+0.2808.*k.*1i-k.^2./2) / (0.01365+0.3455.*k.*1i-k.^2);

f11 = sigma^2./V.^2 - k.^2./mu + (2*1i.*Ck)/mu;

f12 = (k.*(1i + a.*k) + (2 + 1i.*k.*(1-2*a)).*Ck)/mu;

f21 = (a.*k.^2 - 1i.*k.*(1 + 2*a).*Ck)/mu;

f22 = ((8*mu*rs)./V^2 + 4*1i*(1 + 2*a)*(2*1i - k.*(1 - 2*a)).*Ck - k.*(k - 4*1i + 8*a*(1i + a.*k)))/8*mu;

p = [(rs-xtheta^2) 0 (f22 + f11.*rs - f21.*xtheta - f12.*xtheta) 0 (f11.*f22 - f12.*f21)];

rts = roots(p)

K = sort(rts, 'ComparisonMethod','real')

k3 = real(K(3))

k3_visu(j) = k3

end

k = 0.1;

k4_visu(1) = k

for j = 2 : n;

Ck = (0.01365+0.2808.*k.*1i-k.^2./2) / (0.01365+0.3455.*k.*1i-k.^2);

f11 = sigma^2./V.^2 - k.^2./mu + (2*1i.*Ck)/mu;

f12 = (k.*(1i + a.*k) + (2 + 1i.*k.*(1-2*a)).*Ck)/mu;

f21 = (a.*k.^2 - 1i.*k.*(1 + 2*a).*Ck)/mu;

f22 = ((8*mu*rs)./V^2 + 4*1i*(1 + 2*a)*(2*1i - k.*(1 - 2*a)).*Ck - k.*(k - 4*1i + 8*a*(1i + a.*k)))/8*mu;

p = [(rs-xtheta^2) 0 (f22 + f11.*rs - f21.*xtheta - f12.*xtheta) 0 (f11.*f22 - f12.*f21)];

rts = roots(p)

K = sort(rts, 'ComparisonMethod','real')

k4 = real(K(4))

k4_visu(j) = k4

end

plot(1:n,k1_visu) %first root

hold on;

plot(1:n,k2_visu) %second root

plot(1:n,k3_visu) %third root

plot(1:n,k4_visu) %fourth root

hold off;

[Mathieu NOE]

here my code suggestion

as you reuse the same loop , I put it in a function with parameters as arguments

% % >> %inputs
params.a = -1/5;
params.mu = 2; % 2 ≤ µ ≤ 20
params.rs = 1/4;
params.sigma = 1/sqrt(10);
params.xtheta= 0; % xtheta = 0.0 and 0.2;
params.V = 0.05;
%equations and pre-defining vectors/matrices
params.k = 0.5;
params.loops = 10;% nb of loops;% nb of loops
roots = 4;

k_visu = zeros(roots,params.loops);

% main loop / code
for ck = 1:roots
    k_visu(ck,:) = loop(params, ck);
    legend_str{ck} = (['root # : ' num2str(ck)]);
end

plot(1:params.loops,k_visu);grid
legend(legend_str)
xlabel('iteration');

%%

function k_visu = loop(params, root_index)

a = params.a;
mu = params.mu; % 2 ≤ µ ≤ 20
rs = params.rs;
sigma = params.sigma;
xtheta = params.xtheta; % xtheta = 0.0 and 0.2
V = params.V;
%equations and pre-defining vectors/matrices
k = params.k;
loops = params.loops;% nb of loops

k_visu = zeros(1,loops);

    for j = 2 : loops

        Ck = (0.01365+0.2808.*k.*1i-k.^2./2) / (0.01365+0.3455.*k.*1i-k.^2);
        f11 = sigma^2./V.^2 - k.^2./mu + (2*1i.*Ck)/mu;
        f12 = (k.*(1i + a.*k) + (2 + 1i.*k.*(1-2*a)).*Ck)/mu;
        f21 = (a.*k.^2 - 1i.*k.*(1 + 2*a).*Ck)/mu;
        f22 = ((8*mu*rs)./V^2 + 4*1i*(1 + 2*a)*(2*1i - k.*(1 - 2*a)).*Ck - k.*(k - 4*1i + 8*a*(1i + a.*k)))/8*mu;
        p = [(rs-xtheta^2) 0 (f22 + f11.*rs - f21.*xtheta - f12.*xtheta) 0 (f11.*f22 - f12.*f21)];
        rts = roots(p); 
        rts_real_sorted = sort(real(rts));
        k = rts_real_sorted(root_index);
        k_visu(j) = k;
    %     gamma(j) = real (rts(1))./k;
    end

end

[Rodrigo Pena]

Hello !

When I paste your code on matlab, it returns:

unrecognized variable "loop"

[Mathieu NOE]

NB : in my code there is no variable "loop" (without s at the end)

so have you retrieved by mistake the trailing "s" to variable "loops"

or this is because my subfunction was called loop (I changed it to mainloop)

if you have an older matlab the subfunction must be stored in a separate file and not at the end of your main script

[Rodrigo Pena]

my matlab is 2020b.

here is the print:

[Mathieu NOE]

verry funny because I'm running 2020b also and without any error

what if you do a clear all first ?

do you paste the code in the matlab command window or did you save it in a script ?

[Mathieu NOE]

ok

got the same error if you copy / paste in the matlab command window

that method does not work if you have a nested function in your code

save the code in a m file and run it from the editor.

[Rodrigo Pena]

I did clear all, I simply paste the code in the matlab command window.

I am using student version from university, I don´t know if this means something.

That´s weird.

[Rodrigo Pena]

aah... let me see how I can do that... just a minute

[Rodrigo Pena]

I believe I don´t have editor:

This is my screen

[Rodrigo Pena]

Nevermind... I found the thing. LOL

[Rodrigo Pena]

Why we have different results ?

[Mathieu NOE]

ok

good that you find the editor !! (would have surprised me that your release has no editor...)

I got the same plot as you

the first version of the code with 1st root correspond to the lower blue trace

[Rodrigo Pena]

I mean, the first and the fourth root are the same value, but the second and the tird are different.

[Mathieu NOE]

no

your polynomial p is like p = a X^4 + b X^2 + constant

so you have 2 pairs of roots , within each pair the same real / imag values but different sign

here the result of rts while looping :

rts =

-0.0324 +39.9570i

0.0324 -39.9570i

0.0466 - 6.3277i

-0.0466 + 6.3277i

rts =

0.0095 +39.9452i

-0.0095 -39.9452i

0.0719 - 6.3143i

-0.0719 + 6.3143i

rts =

-0.0115 -39.9476i

0.0115 +39.9476i

-0.0684 + 6.3128i

0.0684 - 6.3128i

and so on

so at the end the 4 lines in ypur graph will be symmetrical (except the starting point)

[Rodrigo Pena]

aaah ok !

I got it !

Thank you !

[Rodrigo Pena]

I just implemented gamma, I have made some mistakes, k = imag(rts(1)) not k = real(rts(1)).

I just need to do one more thing... the code above is giving converged values of k and gamma. Now I would like to plot the variation of k and gamma in function of V ( which was equal to 0.05 ). V has to vary from 0.05 up to 5 in increments of 0.05.

Can you help me please?

%inputs

a = -1/5;

mu = 2; % 2 ≤ µ ≤ 20

rs = 1/4;

sigma = 1/sqrt(10);

xtheta = 0; % xtheta = 0.0 and 0.2

V = 0.05;

%equations and pre-defining vectors/matrices

k = 0.1;

k1_visu(1) = k;

gamma1_visu(1) = k;

n = 10; % nb of loops

for j = 2 : n

Ck = (0.01365+0.2808.*k.*1i-k.^2./2) / (0.01365+0.3455.*k.*1i-k.^2);

f11 = sigma^2./V.^2 - k.^2./mu + (2*1i.*Ck)/mu;

f12 = (k.*(1i + a.*k) + (2 + 1i.*k.*(1-2*a)).*Ck)/mu;

f21 = (a.*k.^2 - 1i.*k.*(1 + 2*a).*Ck)/mu;

f22 = ((8*mu*rs)./V^2 + 4*1i*(1 + 2*a)*(2*1i - k.*(1 - 2*a)).*Ck - k.*(k - 4*1i + 8*a*(1i + a.*k)))/8*mu;

p = [(rs-xtheta^2) 0 (f22 + f11.*rs - f21.*xtheta - f12.*xtheta) 0 (f11.*f22 - f12.*f21)];

rts = roots(p);

K = sort(rts, 'ComparisonMethod','real');

k1 = imag(K(1));

k1_visu(j) = k1

gamma1 = real(K(1))/k1;

gamma1_visu(j) = gamma1

end

k = 0.1;

k2_visu(1) = k;

gamma2_visu(1) = k;

for j = 2 : n

Ck = (0.01365+0.2808.*k.*1i-k.^2./2) / (0.01365+0.3455.*k.*1i-k.^2);

f11 = sigma^2./V.^2 - k.^2./mu + (2*1i.*Ck)/mu;

f12 = (k.*(1i + a.*k) + (2 + 1i.*k.*(1-2*a)).*Ck)/mu;

f21 = (a.*k.^2 - 1i.*k.*(1 + 2*a).*Ck)/mu;

f22 = ((8*mu*rs)./V^2 + 4*1i*(1 + 2*a)*(2*1i - k.*(1 - 2*a)).*Ck - k.*(k - 4*1i + 8*a*(1i + a.*k)))/8*mu;

p = [(rs-xtheta^2) 0 (f22 + f11.*rs - f21.*xtheta - f12.*xtheta) 0 (f11.*f22 - f12.*f21)];

rts = roots(p);

K = sort(rts, 'ComparisonMethod','real');

k2 = imag(K(2));

k2_visu(j) = k2

gamma2 = real(K(2))/k2;

gamma2_visu(j) = gamma2

end

k = 0.1;

k3_visu(1) = k;

gamma3_visu(1) = k;

for j = 2 : n

Ck = (0.01365+0.2808.*k.*1i-k.^2./2) / (0.01365+0.3455.*k.*1i-k.^2);

f11 = sigma^2./V.^2 - k.^2./mu + (2*1i.*Ck)/mu;

f12 = (k.*(1i + a.*k) + (2 + 1i.*k.*(1-2*a)).*Ck)/mu;

f21 = (a.*k.^2 - 1i.*k.*(1 + 2*a).*Ck)/mu;

f22 = ((8*mu*rs)./V^2 + 4*1i*(1 + 2*a)*(2*1i - k.*(1 - 2*a)).*Ck - k.*(k - 4*1i + 8*a*(1i + a.*k)))/8*mu;

p = [(rs-xtheta^2) 0 (f22 + f11.*rs - f21.*xtheta - f12.*xtheta) 0 (f11.*f22 - f12.*f21)];

rts = roots(p);

K = sort(rts, 'ComparisonMethod','real');

k3 = imag(K(3));

k3_visu(j) = k3

gamma3 = real(K(3))/k3;

gamma3_visu(j) = gamma3

end

k = 0.1;

k4_visu(1) = k;

gamma4_visu(1) = k;

for j = 2 : n

Ck = (0.01365+0.2808.*k.*1i-k.^2./2) / (0.01365+0.3455.*k.*1i-k.^2);

f11 = sigma^2./V.^2 - k.^2./mu + (2*1i.*Ck)/mu;

f12 = (k.*(1i + a.*k) + (2 + 1i.*k.*(1-2*a)).*Ck)/mu;

f21 = (a.*k.^2 - 1i.*k.*(1 + 2*a).*Ck)/mu;

f22 = ((8*mu*rs)./V^2 + 4*1i*(1 + 2*a)*(2*1i - k.*(1 - 2*a)).*Ck - k.*(k - 4*1i + 8*a*(1i + a.*k)))/8*mu;

p = [(rs-xtheta^2) 0 (f22 + f11.*rs - f21.*xtheta - f12.*xtheta) 0 (f11.*f22 - f12.*f21)];

rts = roots(p);

K = sort(rts, 'ComparisonMethod','real');

k4 = imag(K(4));

k4_visu(j) = k4

gamma4 = real(K(4))/k4;

gamma4_visu(j) = gamma4

end

k_converged = [k1_visu(10), k2_visu(10), k3_visu(10), k4_visu(10)]

gamma_converged = [gamma1_visu(10), gamma2_visu(10), gamma3_visu(10), gamma4_visu(10)]

figure(1)

plot(1:n,k1_visu) %first imaginary root

hold on;

plot(1:n,k2_visu) %second imaginary root

plot(1:n,k3_visu) %third imaginary root

plot(1:n,k4_visu) %fourth imaginary root

grid on;

hold off;

figure(2)

plot(1:n,gamma1_visu)

hold on;

plot(1:n,gamma2_visu)

plot(1:n,gamma3_visu)

plot(1:n,gamma4_visu)

grid on;

hold off;

[Mathieu NOE]

hi

before doing one more iteration, I was surprised by the output of the code.

i believe there is a mistake when you have created the k1, k2, k3 ,k4 varaibles , because these varaibles are not used inside the for loop to update the equation : it uses only "k" and not k1, k2, k3 ,k4

see the lines i have commented and what I have corrected is just below

so this would be my correction - after we will see how to continue.

clc
clear all

%inputs
a = -1/5;
mu = 2; % 2 ≤ µ ≤ 20
rs = 1/4;
sigma = 1/sqrt(10);
xtheta = 0; % xtheta = 0.0 and 0.2
V = 0.05;
%equations and pre-defining vectors/matrices
k = 0.1;
k1_visu(1) = k;
gamma1_visu(1) = k;
n = 10; % nb of loops

for j = 2 : n
    Ck = (0.01365+0.2808.*k.*1i-k.^2./2) / (0.01365+0.3455.*k.*1i-k.^2);
    f11 = sigma^2./V.^2 - k.^2./mu + (2*1i.*Ck)/mu;
    f12 = (k.*(1i + a.*k) + (2 + 1i.*k.*(1-2*a)).*Ck)/mu;
    f21 = (a.*k.^2 - 1i.*k.*(1 + 2*a).*Ck)/mu;
    f22 = ((8*mu*rs)./V^2 + 4*1i*(1 + 2*a)*(2*1i - k.*(1 - 2*a)).*Ck - k.*(k - 4*1i + 8*a*(1i + a.*k)))/8*mu;
    p = [(rs-xtheta^2) 0 (f22 + f11.*rs - f21.*xtheta - f12.*xtheta) 0 (f11.*f22 - f12.*f21)];
    rts = roots(p);
    K = sort(rts, 'ComparisonMethod','real');
%     k1 = imag(K(1)); % it cannot be k1 because this loop is working on
%     "k" and not k1 (otherwise there is no recursion on k value
%     k1_visu(j) = k1
%     gamma1 = real(K(1))/k1;
    k = imag(K(1)); % my correction
    k1_visu(j) = k; % my correction
    gamma1 = real(K(1))/k; % my correction
    gamma1_visu(j) = gamma1;
    
   end
k = 0.1;
k2_visu(1) = k;
gamma2_visu(1) = k;
for j = 2 : n
    Ck = (0.01365+0.2808.*k.*1i-k.^2./2) / (0.01365+0.3455.*k.*1i-k.^2);
    f11 = sigma^2./V.^2 - k.^2./mu + (2*1i.*Ck)/mu;
    f12 = (k.*(1i + a.*k) + (2 + 1i.*k.*(1-2*a)).*Ck)/mu;
    f21 = (a.*k.^2 - 1i.*k.*(1 + 2*a).*Ck)/mu;
    f22 = ((8*mu*rs)./V^2 + 4*1i*(1 + 2*a)*(2*1i - k.*(1 - 2*a)).*Ck - k.*(k - 4*1i + 8*a*(1i + a.*k)))/8*mu;
    p = [(rs-xtheta^2) 0 (f22 + f11.*rs - f21.*xtheta - f12.*xtheta) 0 (f11.*f22 - f12.*f21)];
    rts = roots(p);
    K = sort(rts, 'ComparisonMethod','real'); 
%     k2 = imag(K(2)); % same comment as above
%     k2_visu(j) = k2
%     gamma2 = real(K(2))/k2;
    k = imag(K(2)); % my correction
    k2_visu(j) = k; % my correction
    gamma2 = real(K(2))/k; % my correction
    gamma2_visu(j) = gamma2;
    
   end
k = 0.1;
k3_visu(1) = k;
gamma3_visu(1) = k;
for j = 2 : n
    Ck = (0.01365+0.2808.*k.*1i-k.^2./2) / (0.01365+0.3455.*k.*1i-k.^2);
    f11 = sigma^2./V.^2 - k.^2./mu + (2*1i.*Ck)/mu;
    f12 = (k.*(1i + a.*k) + (2 + 1i.*k.*(1-2*a)).*Ck)/mu;
    f21 = (a.*k.^2 - 1i.*k.*(1 + 2*a).*Ck)/mu;
    f22 = ((8*mu*rs)./V^2 + 4*1i*(1 + 2*a)*(2*1i - k.*(1 - 2*a)).*Ck - k.*(k - 4*1i + 8*a*(1i + a.*k)))/8*mu;
    p = [(rs-xtheta^2) 0 (f22 + f11.*rs - f21.*xtheta - f12.*xtheta) 0 (f11.*f22 - f12.*f21)];
    rts = roots(p);
    K = sort(rts, 'ComparisonMethod','real'); 
%     k3 = imag(K(3)); % same comment as above
%     k3_visu(j) = k3
%     gamma3 = real(K(3))/k3;
    k = imag(K(3)); % my correction
    k3_visu(j) = k % my correction
    gamma3 = real(K(3))/k; % my correction
    gamma3_visu(j) = gamma3;
    
   end
k = 0.1;
k4_visu(1) = k;
gamma4_visu(1) = k;
for j = 2 : n
    Ck = (0.01365+0.2808.*k.*1i-k.^2./2) / (0.01365+0.3455.*k.*1i-k.^2);
    f11 = sigma^2./V.^2 - k.^2./mu + (2*1i.*Ck)/mu;
    f12 = (k.*(1i + a.*k) + (2 + 1i.*k.*(1-2*a)).*Ck)/mu;
    f21 = (a.*k.^2 - 1i.*k.*(1 + 2*a).*Ck)/mu;
    f22 = ((8*mu*rs)./V^2 + 4*1i*(1 + 2*a)*(2*1i - k.*(1 - 2*a)).*Ck - k.*(k - 4*1i + 8*a*(1i + a.*k)))/8*mu;
    p = [(rs-xtheta^2) 0 (f22 + f11.*rs - f21.*xtheta - f12.*xtheta) 0 (f11.*f22 - f12.*f21)];
    rts = roots(p);
    K = sort(rts, 'ComparisonMethod','real'); 
%     k4 = imag(K(4)); % same comment as above
%     k4_visu(j) = k4
%     gamma4 = real(K(4))/k4;
    k = imag(K(4)); % my correction
    k4_visu(j) = k; % my correction
    gamma4 = real(K(4))/k; % my correction
    gamma4_visu(j) = gamma4
    
   end
k_converged = [k1_visu(10), k2_visu(10), k3_visu(10), k4_visu(10)]
gamma_converged = [gamma1_visu(10), gamma2_visu(10), gamma3_visu(10), gamma4_visu(10)]
figure(1)
plot(1:n,k1_visu) %first imaginary root
hold on;
plot(1:n,k2_visu) %second imaginary root
plot(1:n,k3_visu) %third imaginary root
plot(1:n,k4_visu) %fourth imaginary root
grid on;
hold off;
figure(2)
plot(1:n,gamma1_visu)
hold on;
plot(1:n,gamma2_visu)
plot(1:n,gamma3_visu)
plot(1:n,gamma4_visu)
grid on;
hold off;

[Mathieu NOE]

so after this correction, your code and my version are equivalent

the first plot are identical, but no global convergence for all 4 roots

so what do we do ?

my code below :

% I am trying to write a program that calculates aeroelasticiy stuff. I am not sure if the program is calculating what is supposed to. 
% The plan is:
% 1- Guess a value of k between 0 and 1
% 2- Calculate Ck
% 3- Calculate f11
% 4- Calculate f12
% 5- Calculate f21
% 6- Calculate f22
% 7- Calculate the polynomial roots of (p)
% 8- Pick one of the roots of the polynomial p (let´s say for example the first one)
% 9- The new value of k is equal to the real part of the root of the polynomial p (first root)
% 10- Return to step 2 and calculate everthing until k gets convergence.
clc
clear all
% %inputs
% a = -1/5;
% mu = 2; % 2 ≤ µ ≤ 20
% rs = 1/4;
% sigma = 1/sqrt(10);
% xtheta = 0; % xtheta = 0.0 and 0.2
% V = 0.05;
% %equations and pre-defining vectors/matrices
% k = 0.1;
% k1_visu(1) = k;
% gamma1_visu(1) = k;

n = 10; % nb of loops
% % >> %inputs
params.a = -1/5;
params.mu = 2; % 2 ≤ µ ≤ 20
params.rs = 1/4;
params.sigma = 1/sqrt(10);
params.xtheta= 0; % xtheta = 0.0 and 0.2;
params.V = 0.05;
%equations and pre-defining vectors/matrices
params.k = 0.1;
params.loops = 10;% nb of loops;% nb of loops
roots = 4;

k_visu = zeros(roots,params.loops);

% main loop / code
for ck = 1:roots
    k_visu(ck,:) = mainloop(params, ck);
    legend_str{ck} = (['root # : ' num2str(ck)]);
end

plot(1:params.loops,k_visu);grid
legend(legend_str)
xlabel('iteration');

%%

function k_visu = mainloop(params, root_index)

a = params.a;
mu = params.mu; % 2 ≤ µ ≤ 20
rs = params.rs;
sigma = params.sigma;
xtheta = params.xtheta; % xtheta = 0.0 and 0.2
V = params.V;
%equations and pre-defining vectors/matrices
k = params.k;
loops = params.loops;% nb of loops

k_visu = zeros(1,loops);

    for j = 2 : loops
        
    Ck = (0.01365+0.2808.*k.*1i-k.^2./2) / (0.01365+0.3455.*k.*1i-k.^2);
    f11 = sigma^2./V.^2 - k.^2./mu + (2*1i.*Ck)/mu;
    f12 = (k.*(1i + a.*k) + (2 + 1i.*k.*(1-2*a)).*Ck)/mu;
    f21 = (a.*k.^2 - 1i.*k.*(1 + 2*a).*Ck)/mu;
    f22 = ((8*mu*rs)./V^2 + 4*1i*(1 + 2*a)*(2*1i - k.*(1 - 2*a)).*Ck - k.*(k - 4*1i + 8*a*(1i + a.*k)))/8*mu;
    p = [(rs-xtheta^2) 0 (f22 + f11.*rs - f21.*xtheta - f12.*xtheta) 0 (f11.*f22 - f12.*f21)];
    rts = roots(p);
       
    K = sort(rts, 'ComparisonMethod','real');
    k = imag(K(root_index));

    k_visu(j) = k;

    %     gamma(j) = real (rts(1))./k;
    end

end

[Rodrigo Pena]

oh gosh...

[Mathieu NOE]

exactly....

[Rodrigo Pena]

how about raising the number of interactions ? let´s say 20

[Rodrigo Pena]

let´s get rid of the other three roots and just work with 1. When we manage to solve the problem, we expand for the rest of the roots

[Mathieu NOE]

so this is the code for looping over V

% I am trying to write a program that calculates aeroelasticiy stuff. I am not sure if the program is calculating what is supposed to. 
% The plan is:
% 1- Guess a value of k between 0 and 1
% 2- Calculate Ck
% 3- Calculate f11
% 4- Calculate f12
% 5- Calculate f21
% 6- Calculate f22
% 7- Calculate the polynomial roots of (p)
% 8- Pick one of the roots of the polynomial p (let´s say for example the first one)
% 9- The new value of k is equal to the real part of the root of the polynomial p (first root)
% 10- Return to step 2 and calculate everthing until k gets convergence.
clc
clear all

n = 10; % nb of loops
% % >> %inputs
params.a = -1/5;
params.mu = 2; % 2 ≤ µ ≤ 20
params.rs = 1/4;
params.sigma = 1/sqrt(10);
params.xtheta= 0; % xtheta = 0.0 and 0.2;
params.V = 0.05;
%equations and pre-defining vectors/matrices
params.k = 0.1;
params.loops = 10;% nb of loops;% nb of loops
roots = 4;

%%% loop with root number %%%%
% k_visu = zeros(roots,params.loops);
% 
% % main loop / code
% for ck = 1:roots
%     k_visu(ck,:) = mainloop(params, ck);
%     legend_str{ck} = (['root # : ' num2str(ck)]);
% end
% 
% plot(1:params.loops,k_visu);grid
% legend(legend_str)
% xlabel('iteration');

%% loop vs V parameter %%
% V has to vary from 0.05 up to 5 in increments of 0.05.
V = (0.05:0.05:5);
k1_converged_visu = zeros(length(V),1);
k2_converged_visu = zeros(length(V),1);
k3_converged_visu = zeros(length(V),1);
k4_converged_visu = zeros(length(V),1);

% main loop / code
for ck = 1:length(V)
    params.V = V(ck);

    k1_visu = mainloop(params, 1);
    k1_converged_visu(ck) = k1_visu(end);
    
    k2_visu = mainloop(params, 2);
    k2_converged_visu(ck) = k2_visu(end);
    
    k3_visu = mainloop(params, 3);
    k3_converged_visu(ck) = k3_visu(end);
    
    k4_visu = mainloop(params, 4);
    k4_converged_visu(ck) = k4_visu(end);
end
xaxis_data = 1:length(V);
plot(xaxis_data,[k1_converged_visu k2_converged_visu k3_converged_visu k4_converged_visu]);grid
xlabel('V');

%%

function k_visu = mainloop(params, root_index)

a = params.a;
mu = params.mu; % 2 ≤ µ ≤ 20
rs = params.rs;
sigma = params.sigma;
xtheta = params.xtheta; % xtheta = 0.0 and 0.2
V = params.V;
%equations and pre-defining vectors/matrices
k = params.k;
loops = params.loops;% nb of loops

k_visu = zeros(1,loops);

    for j = 2 : loops
        
    Ck = (0.01365+0.2808.*k.*1i-k.^2./2) / (0.01365+0.3455.*k.*1i-k.^2);
    f11 = sigma^2./V.^2 - k.^2./mu + (2*1i.*Ck)/mu;
    f12 = (k.*(1i + a.*k) + (2 + 1i.*k.*(1-2*a)).*Ck)/mu;
    f21 = (a.*k.^2 - 1i.*k.*(1 + 2*a).*Ck)/mu;
    f22 = ((8*mu*rs)./V^2 + 4*1i*(1 + 2*a)*(2*1i - k.*(1 - 2*a)).*Ck - k.*(k - 4*1i + 8*a*(1i + a.*k)))/8*mu;
    p = [(rs-xtheta^2) 0 (f22 + f11.*rs - f21.*xtheta - f12.*xtheta) 0 (f11.*f22 - f12.*f21)];
    rts = roots(p);
       
    K = sort(rts, 'ComparisonMethod','real');
    k = imag(K(root_index));

    k_visu(j) = k;

    %     gamma(j) = real (rts(1))./k;
    end

end

[Rodrigo Pena]

Ok, thank you so much !

I will talk to my professor, I don´t know why we are not getting converged values. What the reck that code was doing when we started to get converged values? I was so happy I felt I was going in the right direction but now, everything is gone ):

Thank you again !

[Mathieu NOE]

No , it's not lost

the problem lies around real or imag part of the roots

it was converging with real part of roots;

now that you mentionned : I have made some mistakes, k = imag(rts(1)) not k = real(rts(1)).

we change to imag and it's not converging anymore

that's the point to double check

[Rodrigo Pena]

Ok, let me see here

[Rodrigo Pena]

That is my assingment. I belive I wrote everything correct. Thank you so much for helping me !!!

[Rodrigo Pena]

Why k1_visu is starting already full in the first interaction?

[Rodrigo Pena]

Hello !

I have great news !

if you change flow velocity (V) to 0.52, you will see converged values of k !!

So now, our value of V will start from 0.52 up to 5 in increments of 0.01.

[Mathieu NOE]

OK

so we keep the code as last modification (k = imag and not real part of the roots)

    K = sort(rts, 'ComparisonMethod','real');
    k = imag(K(root_index));

but only the first root seems to converge, the 3 others keep oscillating , you confirm ?

I increased the number of loops to 15

[Mathieu NOE]

and if this is validated then, the final code for plotting k vs V will be :

% I am trying to write a program that calculates aeroelasticiy stuff. I am not sure if the program is calculating what is supposed to. 
% The plan is:
% 1- Guess a value of k between 0 and 1
% 2- Calculate Ck
% 3- Calculate f11
% 4- Calculate f12
% 5- Calculate f21
% 6- Calculate f22
% 7- Calculate the polynomial roots of (p)
% 8- Pick one of the roots of the polynomial p (let´s say for example the first one)
% 9- The new value of k is equal to the real part of the root of the polynomial p (first root)
% 10- Return to step 2 and calculate everthing until k gets convergence.
clc
clear all

% % >> %inputs
params.a = -1/5;
params.mu = 2; % 2 ≤ µ ≤ 20
params.rs = 1/4;
params.sigma = 1/sqrt(10);
params.xtheta= 0; % xtheta = 0.0 and 0.2;
%equations and pre-defining vectors/matrices
params.k = 0.1;
params.loops = 15;% nb of loops;% nb of loops


%% loop vs V parameter %%
% V has to vary from 0.52 up to 5 in increments of 0.05.
V = (0.52:0.05:5);
k1_converged_visu = zeros(length(V),1);
% k2_converged_visu = zeros(length(V),1);
% k3_converged_visu = zeros(length(V),1);
% k4_converged_visu = zeros(length(V),1);

% main loop / code
for ck = 1:length(V)
    params.V = V(ck);

    k1_visu = mainloop(params, 1); % 1 here means 1st root is computed
    k1_converged_visu(ck) = k1_visu(end);
    
%     k2_visu = mainloop(params, 2);
%     k2_converged_visu(ck) = k2_visu(end);
%     
%     k3_visu = mainloop(params, 3);
%     k3_converged_visu(ck) = k3_visu(end);
%     
%     k4_visu = mainloop(params, 4);
%     k4_converged_visu(ck) = k4_visu(end);
end
% plot(V,[k1_converged_visu k2_converged_visu k3_converged_visu k4_converged_visu]);grid
plot(V,k1_converged_visu);grid
title(' k vs V ');
ylabel('k (1st root)');
xlabel('V');

%%

function k_visu = mainloop(params, root_index)

a = params.a;
mu = params.mu; % 2 ≤ µ ≤ 20
rs = params.rs;
sigma = params.sigma;
xtheta = params.xtheta; % xtheta = 0.0 and 0.2
V = params.V;
%equations and pre-defining vectors/matrices
k = params.k;
loops = params.loops;% nb of loops

k_visu = zeros(1,loops);

    for j = 2 : loops
        
    Ck = (0.01365+0.2808.*k.*1i-k.^2./2) / (0.01365+0.3455.*k.*1i-k.^2);
    f11 = sigma^2./V.^2 - k.^2./mu + (2*1i.*Ck)/mu;
    f12 = (k.*(1i + a.*k) + (2 + 1i.*k.*(1-2*a)).*Ck)/mu;
    f21 = (a.*k.^2 - 1i.*k.*(1 + 2*a).*Ck)/mu;
    f22 = ((8*mu*rs)./V^2 + 4*1i*(1 + 2*a)*(2*1i - k.*(1 - 2*a)).*Ck - k.*(k - 4*1i + 8*a*(1i + a.*k)))/8*mu;
    p = [(rs-xtheta^2) 0 (f22 + f11.*rs - f21.*xtheta - f12.*xtheta) 0 (f11.*f22 - f12.*f21)];
    rts = roots(p);
       
    K = sort(rts, 'ComparisonMethod','real');
    k = imag(K(root_index));

    k_visu(j) = k;

    %     gamma(j) = real (rts(1))./k;
    end

end

and this is the plot I got :

[Rodrigo Pena]

Yes, if you increase the flow velocity (V) to 0.783 you will see converged values for the 3 roots. The only problem is, as we increase the value of V, the value of k for the fourth root does not converge however, gamma for the fourth root converges.

That´s weird because gamma should converge automatically when k converges.

Thank you so much !

[Mathieu NOE]

you're welcome

good luck for the future

[Rodrigo Pena]

Hello Mathieu,

Sorry to bother you again.

Can you please inplement the %gamma(j) = real (rts(1))./k; ?

Thank you !

[Mathieu NOE]

hello

yes of course

but I also printed your assignment

i believe for f22 there is an error when you divide by (8*mu)

so I send you back the two corrected codes

you write :

f22 = ((8*mu*rs)./V^2 + 4*1i*(1 + 2*a)*(2*1i - k.*(1 - 2*a)).*Ck - k.*(k - 4*1i + 8*a*(1i + a.*k)))/8*mu

whereas it should be, if I read correctly your assignment :

f22 = ((8*mu*rs)./V^2 + 4*1i*(1 + 2*a)*(2*1i - k.*(1 - 2*a)).*Ck - k.*(k - 4*1i + 8*a*(1i + a.*k)))/(8*mu);

[Rodrigo Pena]

Thank you very much ! I will talk to my professor and explain what is going on, k4 still does not converge, while the other roots converges.

[Rodrigo Pena]

Hello Mathieu,

Sorry to bother you again.

Can you do this for the rest of the roots ? I belive that you did just for the first one.

I tried to implement to the other roots but I keep obtaining the same curve with the exacly values of the first root.

Thank you so much !

[Mathieu NOE]

hello Rodrigo

here you are :

% I am trying to write a program that calculates aeroelasticiy stuff. I am not sure if the program is calculating what is supposed to. 
% The plan is:
% 1- Guess a value of k between 0 and 1
% 2- Calculate Ck
% 3- Calculate f11
% 4- Calculate f12
% 5- Calculate f21
% 6- Calculate f22
% 7- Calculate the polynomial roots of (p)
% 8- Pick one of the roots of the polynomial p (let´s say for example the first one)
% 9- The new value of k is equal to the real part of the root of the polynomial p (first root)
% 10- Return to step 2 and calculate everthing until k gets convergence.
clc
clear all

% % >> %inputs
params.a = -1/5;
params.mu = 2; % 2 ≤ µ ≤ 20
params.rs = 1/4; % NB : rs means r squared (no need to put ^2 in root equation
params.sigma = 1/sqrt(10);
params.xtheta= 0; % xtheta = 0.0 and 0.2;
%equations and pre-defining vectors/matrices
params.k = 0.1;
params.loops = 15;% nb of loops;% nb of loops
root_index = 1:4;  % 1:4; %  root


%% loop vs V parameter %%
% V has to vary from 0.52 up to 5 in increments of 0.05.
V = (0.52:0.05:5);

for ci = 1:length(root_index)
    
    % main loop / code
    for ck = 1:length(V)
        params.V = V(ck);

        [k_visu,gamma] = mainloop(params, root_index(ci)); 
        k_converged_visu(ci,ck) = k_visu(end);
        gamma_converged_visu(ci,ck) = gamma(end);

    end
    
end



figure(1),plot(V,k_converged_visu);grid
title(' k vs V ');
ylabel('k (1st root)');
xlabel('V');

figure(2),plot(V,gamma_converged_visu);grid
title(' gamma vs V ');
ylabel('gamma (1st root)');
xlabel('V');


%%

function [k_visu,gamma] = mainloop(params, root_index)

a = params.a;
mu = params.mu; % 2 ≤ µ ≤ 20
rs = params.rs;
sigma = params.sigma;
xtheta = params.xtheta; % xtheta = 0.0 and 0.2
V = params.V;
%equations and pre-defining vectors/matrices
k = params.k;
loops = params.loops;% nb of loops

k_visu = zeros(1,loops);

    for j = 2 : loops
        
    Ck = (0.01365+0.2808.*k.*1i-k.^2./2) / (0.01365+0.3455.*k.*1i-k.^2);
    f11 = sigma^2./V.^2 - k.^2./mu + (2*1i.*Ck)/mu; %ok
    f12 = (k.*(1i + a.*k) + (2 + 1i.*k.*(1-2*a)).*Ck)/mu; %ok
    f21 = (a.*k.^2 - 1i.*k.*(1 + 2*a).*Ck)/mu; %ok
%     f22 = ((8*mu*rs)./V^2 + 4*1i*(1 + 2*a)*(2*1i - k.*(1 - 2*a)).*Ck - k.*(k - 4*1i + 8*a*(1i + a.*k)))/8*mu;  % not ok - denominator is missing ( ) 
    f22 = ((8*mu*rs)./V^2 + 4*1i*(1 + 2*a)*(2*1i - k.*(1 - 2*a)).*Ck - k.*(k - 4*1i + 8*a*(1i + a.*k)))/(8*mu); %ok
    p = [(rs-xtheta^2) 0 (f22 + f11.*rs - f21.*xtheta - f12.*xtheta) 0 (f11.*f22 - f12.*f21)]; %ok
   rts = roots(p);
       
    K = sort(rts, 'ComparisonMethod','real');
    k = imag(K(root_index));

    k_visu(j) = k;
%     gamma(j) = real (rts(root_index))./k;

    gamma(j) = real (K(root_index))./k;
    end

end