Write a loop to plot a model graph

Hello,

I wrote this code build a model to find the velocities.

The observed data are plotted in blue and model is in red asteric marks. The model graph also should be a curve that lies closer to observed data with the same shape, however for some reason i get a staright line for this.

Could someone please check this one out and help me figure out my mistake.

I have attached the graph of VL vs T. and instructions i used to write this code.

%% Observed Wave dispersion data;
Tob = (10:5:90); % s (time)
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 (observed velocities)
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

Answer 1

[Alan Stevens]

A little more like this. You need a smaller value of dv and a while loop. The fit isn't very good!

%% Observed Wave dispersion data;
Tob = (10:5:90); % s (time)
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 (observed velocities)
figure (2);
plot(Tob,VLob,'b+') %observed graph
hold on;
grid
%%
T0=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.0005; % used instead of 0.005 of condition VL = VL+0.005 
itmax = 1000;
for T=T0:dt:Tmax
        its = 0;
        delta = 1;
        while delta>0.001 && its<itmax
            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);
            VL = VL+dv;
            its = its+1;
            delta = abs(A-B);
        end
        plot(T,VL,'r*')

end

[Anitha Limann]

hello,

Thank you very much for your kind response.

I used your code with slightly defferent observed data, and non-rounded off vs1, vs2 and z values. However I get the graph as in following image. I have copy and pasted the code i used. Could you please tell me if im missing something?

Also I'd be really grateful if you could explain me why you used "its = 0" and "itmax=1000".

%% Observed Love wave dispersion data;
Tob = (10:5:90); % s
VLob = [3.69 3.78 3.95 4.07 4.21 4.33 4.42 4.49 4.56 4.61 4.65 4.67 4.71 4.75 4.78 4.82 4.86]; %km/s
figure (2);
plot(Tob,VLob,'b+')
hold on;
grid

%%
T0=10; % s 
Tmax = 90; %s
VL = 3.61; %km
D1=2.98;
D2=3.28;
z = 39.9968241272153
vs1 = 3.60843918243516
vs2 = 4.70614826532137
dt=5;
dv=0.0005; % used instead of 0.005 of condition VL = VL+0.005 
itmax = 1000;
for T=T0:dt:Tmax
        its = 0;
        delta = 1;
        while delta>0.001 && its<itmax
            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);
            VL = VL+dv;
            its = its+1;
            delta = abs(A-B);
        end
        plot(T,VL,'r*')
end

[Alan Stevens]

You need an even smaller value of dv!

its represents the number of iterations taken by the while loop. itmax is the maximum allowed number of such iterations. If you don't limit the number of iterations it's possible that the while loop might run forever, never reaching the specified convergence.

I've included a section in the following that displays any times for which its = itmax. you can't trust the result if this occurs, because the convergence criterion probably hasn't been satisfied.

%% Observed Love wave dispersion data;
Tob = (10:5:90); % s
VLob = [3.69 3.78 3.95 4.07 4.21 4.33 4.42 4.49 4.56 4.61 4.65 4.67 4.71 4.75 4.78 4.82 4.86]; %km/s
figure (2);
plot(Tob,VLob,'b+')
hold on;
grid

%%
T0=10; % s 
Tmax = 90; %s
VL = 3.61; %km
D1=2.98;
D2=3.28;
z = 39.9968241272153;
vs1 = 3.60843918243516;
vs2 = 4.70614826532137;
dt=5;
dv=0.0001; % used instead of 0.005 of condition VL = VL+0.005 
itmax = round(1/dv);
for T=T0:dt:Tmax
        its = 0;
        delta = 1;
        while delta>0.001 && its<itmax
            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);
            VL = VL+dv;
            its = its+1;
            if its==itmax  % Can't trust that abs(A-B)<0.001 if its = itmax
                disp(T)
            end
            delta = abs(A-B);
        end
        plot(T,VL,'r*')
end

[Anitha Limann]

Hello,

Thank you so much. That does make sence for why i would only get the expected graph for rounded off values.

This helps a lot.

Dilini