Borrar filtros
Borrar filtros

Write a loop to plot a model graph

1 visualización (últimos 30 días)
Anitha Limann
Anitha Limann el 21 de Oct. de 2021
Comentada: Anitha Limann el 21 de Oct. de 2021
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

Respuesta aceptada

Alan Stevens
Alan Stevens el 21 de Oct. de 2021
Editada: Alan Stevens el 21 de Oct. de 2021
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 comentarios
Alan Stevens
Alan Stevens el 21 de Oct. de 2021
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 el 21 de Oct. de 2021
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

Iniciar sesión para comentar.

Más respuestas (0)

Categorías

Más información sobre Startup and Shutdown en Help Center y File Exchange.

Etiquetas

Community Treasure Hunt

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

Start Hunting!

Translated by