Write a loop to plot a model graph

2 次查看(过去 30 天)
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

采纳的回答

Alan Stevens
Alan Stevens 2021-10-21
编辑:Alan Stevens 2021-10-21
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
  3 个评论
Alan Stevens
Alan Stevens 2021-10-21
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
Anitha Limann 2021-10-21
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

请先登录,再进行评论。

更多回答(0 个)

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by