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

Question

Anitha Limann 2021-10-14

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1563571-how-to-write-a-loop-that-finds-the-best-fit-value-for-each-given-input-pair-and-plot-all-the-resulte

评论： Mathieu NOE 2021-10-28

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.

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)

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

6 个评论
显示 4更早的评论隐藏 4更早的评论

Mathieu NOE 2021-10-14

在 MATLAB Online 中打开

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 2021-10-14

geo.jpg

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 2021-10-15

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1563571-how-to-write-a-loop-that-finds-the-best-fit-value-for-each-given-input-pair-and-plot-all-the-resulte#answer_809526

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

6 个评论
显示 4更早的评论隐藏 4更早的评论

Mathieu NOE 2021-10-18

在 MATLAB Online 中打开

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
    

Mathieu NOE 2021-10-25

在 MATLAB Online 中打开

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');