Write a loop to plot a model graph

1 回表示 (過去 30 日間)
Anitha Limann
Anitha Limann 2021 年 10 月 21 日
コメント済み: Anitha Limann 2021 年 10 月 21 日
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 件)

カテゴリ

Help Center および File ExchangeStartup and Shutdown についてさらに検索

タグ

Community Treasure Hunt

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

Start Hunting!

Translated by