how to write a loop that finds the best fit value for each given input pair and plot all the resulted points?

Hello,

I have t and vl observed data that i plotted first and then find the theoritical points and plot them in the same graph.

To construct the model to find the theoritical values i have to write a code according to following steps.

1. select a starting pair for t1 and vl1. (in this case t=10 and vl = 3.61).
2. Calculate A
3. Calculate B
4. Calculate A-B
5. Test abs(A-B); if abs(A-B) < 0.001 stop, then go to 8. if not go to 6
6. increment vl=vl+0.005
7. goto 2.
8. increment t=t+5
9. goto 1
10. finally plot each pair for A=B (approximately, within given tolerance)

Below is the draft loop I wrote but i am not sure how to actually make it work. Please help me.

for T=10; % tmax = 90;
    for VL=3.61 
    A = atan((vs2.^2).*D2.*(sqrt(1-(VL.^2/vs2.^2)))./(vs1.^2.)*D1.*(sqrt(1-(VL.^2/vs1.^2))));
    B = (((2.*pi.*z)./(VL.*T)).*(sqrt((VL.^2/vs1.^2)-1)));
        while abs(A-B)<0.001
             T=T+5;
             A=A
             B=B
             break
         end
         while abs(A-B)>=0.001
             VL=VL+0.005;
             A=A
             B=B
         end
    end
end

[Mathieu NOE]

hello Anita

can you share the data ? ( t and vl observed data)

[Mathieu NOE]

also,

what are vs1, vs2 and z ?

tx

[Anitha Limann]

vs1, vs2, d1,d2 and z are just constant values.

vs1=3.6

vs2=4.7

d1=2.9

d2=3.2

z=40

Thank you!

[Anitha Limann]

sure,

T = [10:5:90];

VL = 3.6,3.7,3.9,4.1,4.2,4.3,4.4,4.4,4.5,4.6,4.65,4.7,4.7,4.75,4.8,4.82,4.85];

Thank you

[Mathieu NOE]

hello

I have a problem with your equations

see the two examples below . depending if vl>vs1 or vl<vs1 , one output is real and the other is complex , so there is no way to compare them

are you sure that there is no mistake (like missing abs somewhere) in the equations ?

code

clc
clearvars
% T = [10:5:90];
% VL = [3.6;3.7,3.9,4.1,4.2,4.3,4.4,4.4,4.5,4.6,4.65,4.7,4.7,4.75,4.8,4.82,4.85];

% vs1, vs2, d1,d2 and z are just constant values. 
vs1=3.6;
vs2=4.7;
d1=2.9;
d2=3.2;
z=40;


% select a starting pair for t1 and vl1. (in this case t=10 and vl = 3.61). 
% Calculate A
% Calculate B
% Calculate A-B
% Test abs(A-B);  if abs(A-B) < 0.001 stop, then go to 8. if not go to 6
% increment vl=vl+0.005
% goto 2.
% increment t=t+5
% goto 1
% finally plot each pair for A=B (approximately, within given tolerance)

% case 1 
t=10 
vl = 3.61 % is slightly above vs1=3.6;

    A = atan((vs2.^2).*d2.*(sqrt(1-(vl.^2/vs2.^2)))./(vs1.^2.)*d1.*(sqrt(1-(vl.^2/vs1.^2))))
    B = (((2.*pi.*z)./(vl.*t)).*(sqrt((vl.^2/vs1.^2)-1)))

% case 2 
t=10 
vl = 3.51 % is slightly below vs1=3.6;
    
    A = atan((vs2.^2).*d2.*(sqrt(1-(vl.^2/vs2.^2)))./(vs1.^2.)*d1.*(sqrt(1-(vl.^2/vs1.^2))))
    B = (((2.*pi.*z)./(vl.*t)).*(sqrt((vl.^2/vs1.^2)-1)))

gives :

case 1 :

t = 10

vl = 3.6100

A = 0.0000 + 0.9856i

B = 0.5193

case 2 :

t = 10

vl = 3.5100

A = 1.1665

B = 0.0000 + 1.5911i

[Anitha Limann]

Attached is the given question paper. I am supposed to use a value slightly higher than vs1 and lower than observed vL1.

Please see the attachment.

Answer 1

[Mathieu NOE]

hello Anitha

so yes there was a small bug in your equations

when you do the computation of A, you fliiped the terms

(vs1.^2.)*d1.*(sqrt(1-(VL.^2/vs1.^2)))

must be replaced by

(vs1^2)*d1*(sqrt((vl^2/vs1^2)-1))

that's the reason I got A complex numbers and not reals numbers because the ter under the square root was negative

I am still working on the code ....

[Anitha Limann]

hi thank you so much. please let me know how i can make this code work to plot model data.

[Mathieu NOE]

hello

this is the code

I added a plot to see how the error (err = abs(A-B)) evolves at we iterate

I found that the increment step for vl (vl = vl + 0.005) is too coarse im my opinion, so the minimal error will not be low enough to reach the condition abs(A-B) < 0.001

so I tried with a smaller increment and it works fine then

this is a plot of with vl = vl + 0.005

you see the min error is above 10^-3 so we do not meet abs(A-B) < 0.001 and the code will not iterate

With the reduced increment vl = vl + 0.0001 , the code works fine and we get the results for all 17 points

my code is :

clc
clearvars
close all

T = [10:5:90];
VL = [3.6,3.7,3.9,4.1,4.2,4.3,4.4,4.4,4.5,4.6,4.65,4.7,4.7,4.75,4.8,4.82,4.85];

% vs1, vs2, d1,d2 and z are just constant values. 
vs1=3.6;
vs2=4.7;
d1=2.9;
d2=3.2;
z=40;


% select a starting pair for t1 and vl1. (in this case t=10 and vl = 3.61). 
% Calculate A
% Calculate B
% Calculate A-B
% Test abs(A-B);  if abs(A-B) < 0.001 stop, then go to 8. if not go to 6
% increment vl=vl+0.005
% goto 2.
% increment t=t+5
% goto 1
% finally plot each pair for A=B (approximately, within given tolerance)

    
    %start with pair 1
    k = 1;
    t = T(k);
    vl = VL(k);
    
    % initialization 
    ci = 1;
    vl_svg = [];
    svg = 1;
    
    figure(1)
    loglog(svg)
    hold on

    % for debug only
    max_samples = 10000;
    
    %% main loop 
    
    while 1 % infinite loop until one "break" condition is reached 
        
                %safety break to avoid infinite loop forever in case of main loop malfunction
                % used simply to stop even when vl = vl + 0.005 is used and
                % will not meet the (abs(AminusB) < 0.001) condition
                if ci > max_samples % 
                    
                    %%%%%%%%%%% debug only
                   loglog(svg)
                     xlabel('loop iterations')
                   ylabel('err = abs(A-B)')                    
                    ci = 1; % reset counter
                    svg = 1; % reset svg vector
                    %%%%%%%%%%% debug only

                    break

                end

            % computa A, B , A-B
            numA = (vs2^2)*d2*(sqrt(1-(vl^2/vs2^2)));
            denA = (vs1^2)*d1*(sqrt((vl^2/vs1^2)-1)); % correction here
            A = atan(numA/denA);
            B = (((2.*pi.*z)./(vl.*t)).*(sqrt((vl.^2/vs1.^2)-1)));
            AminusB = A-B;
            
            % increment vl
%              vl = vl + 0.005; % increment 
            vl = vl + 0.0001; % I had to reduce the increment otherwise we could never meet abs(AminusB) < 0.001
         
            %%%%%%%%%%%%%% debug only %%%%%%%%%%%%%%
            svg(ci) = abs(AminusB); % just for fun ; plot the convergence curve for each run
            ci = ci+1; % increment loop index 
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            
            if abs(AminusB) < 0.001

                % save results (optimal value of vl is stored at k'th iteration in vl_svg(k)
                vl_svg(k) = vl;
                ci = 1; % reset counter
                
                  % debug only
                   loglog(svg)
                   xlabel('loop iterations')
                   ylabel('err = abs(A-B)')
                   svg = 1;% reset svg
                
                % update input data (next T and VL value)
    %             t = t+5; can be done by  : k = k+1;  and t = T(k); % pick from the list
                    k = k+1;
                    if k > length(T) % stop condition : exit (all iterations completed)
                        break
                    end
                    t = T(k);

            end
        end
        
    % display data in command window (if not empty)
    if ~isempty(vl_svg)
        vl_svg
    end

%     % results obtained are : 
%    
% vl_svg =
% 
%   Columns 1 through 9
% 
%     3.6745    3.7531    3.8484    3.9518    4.0549    4.1510    4.2358    4.3077    4.3674
% 
%   Columns 10 through 17
% 
%     4.4163    4.4562    4.4889    4.5159    4.5382    4.5568    4.5725    4.5858
    

[Anitha Limann]

Hello,

I wrote this code that i get somewhat similar to the expected graph. As you can see the model graph is plotted as red stars. The model graph also should be a curve that lies closer to observed with the same shape but for some reason i get a staright line for this. If you do not mind could you please check this one out and help me figure out my mistake.

The graph of VL vs T.

%% Observed Love wave dispersion data;
Tob = (10:5:90); % s
VLob = [3.6 3.7 3.9 4.0 4.2 4.3 4.4 4.45 4.5 4.6 4.65 4.7 4.75 4.8 4.82 4.85 4.9]; %km/s
figure (2);
plot(Tob,VLob,'b+') %observed graph
hold on;
grid

%%
T=10; % s 
Tmax = 90; %s
VL = 3.61; %km

vs1=3.6;
vs2=4.7;
d1=2.9;
d2=3.2;
z=40;

dt=5;
dv=0.05: % used instead of 0.005 of condition "VL = VL+0.005" 
for T=T:dt:Tmax
    A = atan((vs2^2)*D2*(sqrt(1-(VL^2/vs2^2)))/(vs1^2)*D1*(sqrt((VL^2/vs1^2)-1)));
    B = (((2*pi*z)/(VL*T))*(sqrt((VL^2/vs1^2)-1)));
    if abs(A-B)<0.001
        figure (2)
        plot(T,VL,'r*')
        hold on;
       break
    else 
        VL=VL+dv;
        figure (2)
        plot(T,VL,'r*')
    end
end

[Mathieu NOE]

hello

for the time being I simply added the plot of data to my code

I get this plot (not a straight line)

code :

clc
clearvars
close all

T = [10:5:90];
% VL = [3.6,3.7,3.9,4.1,4.2,4.3,4.4,4.4,4.5,4.6,4.65,4.7,4.7,4.75,4.8,4.82,4.85];
VL = [3.6 3.7 3.9 4.0 4.2 4.3 4.4 4.45 4.5 4.6 4.65 4.7 4.75 4.8 4.82 4.85 4.9]; %km/s

% vs1, vs2, d1,d2 and z are just constant values. 
vs1=3.6;
vs2=4.7;
d1=2.9;
d2=3.2;
z=40;


    
    %start with pair 1
    k = 1;
    t = T(k);
    vl = VL(k);
    
    % initialization 
    ci = 1;
    vl_svg = [];
    svg = 1;
    
    figure(1)
    loglog(svg)
    hold on

    % for debug only
    max_samples = 10000;
    
    %% main loop 
    
    while 1 % infinite loop until one "break" condition is reached 
        
                %safety break to avoid infinite loop forever in case of main loop malfunction
                % used simply to stop even when vl = vl + 0.005 is used and
                % will not meet the (abs(AminusB) < 0.001) condition
                if ci > max_samples % 
                    
                    %%%%%%%%%%% debug only
                   loglog(svg)
                     xlabel('loop iterations')
                   ylabel('err = abs(A-B)')                    
                    ci = 1; % reset counter
                    svg = 1; % reset svg vector
                    %%%%%%%%%%% debug only

                    break

                end

            % computa A, B , A-B
            numA = (vs2^2)*d2*(sqrt(1-(vl^2/vs2^2)));
            denA = (vs1^2)*d1*(sqrt((vl^2/vs1^2)-1)); % correction here
            A = atan(numA/denA);
            B = (((2.*pi.*z)./(vl.*t)).*(sqrt((vl.^2/vs1.^2)-1)));
            AminusB = A-B;
            
            % increment vl
%              vl = vl + 0.005; % increment 
             vl = vl + 0.0001; % I had to reduce the increment otherwise we could never meet abs(AminusB) < 0.001
         
            %%%%%%%%%%%%%% debug only %%%%%%%%%%%%%%
            svg(ci) = abs(AminusB); % just for fun ; plot the convergence curve for each run
            ci = ci+1; % increment loop index 
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            
            if abs(AminusB) < 0.001

                % save results (optimal value of vl is stored at k'th iteration in vl_svg(k)
                vl_svg(k) = vl;
                ci = 1; % reset counter
                
                  % debug only
                   loglog(svg)
                   xlabel('loop iterations')
                   ylabel('err = abs(A-B)')
                   svg = 1;% reset svg
                
                % update input data (next T and VL value)
    %             t = t+5; can be done by  : k = k+1;  and t = T(k); % pick from the list
                    k = k+1;
                    if k > length(T) % stop condition : exit (all iterations completed)
                        break
                    end
                    t = T(k);

            end
        end
        
    % display data in command window (if not empty)
    if ~isempty(vl_svg)
        vl_svg
    end

   
figure (2);
plot(T,VL,'b-+',T,vl_svg,'r-*') %observed and computed data
legend('observed data','computed data');

[Anitha Limann]

Thank you so much for your help.

This make more sense now.

[Mathieu NOE]

My pleasure !