How to remove the straight line coming at the x-axis at x = -0.11532 in my code
1 visualización (últimos 30 días)
Mostrar comentarios más antiguos
vetri veeran
el 21 de Abr. de 2017
Comentada: vetri veeran
el 22 de Abr. de 2017
Hi all,
In my code, I want to remove the straight line coming at the x-axis at x = -0.11532.
I attached my code below,
B = 1e-4;
sigma_on = 0.45;
x_on = 0.06;
sigma_p = 4e-5;
A = 1e-10;
sigma_off = 0.013;
x_off = 0.4;
G_m = 0.025;
a = 7.2e-6;
b = 4.7;
beta = 500;
rho = 1e-3;
v_m = 1;
k=50;
t = -1:0.001:1;
for x = 1:length(t)
G(x) = G_m*t(x)+a*exp(b*sqrt(v_m/(1+exp(-v_m/rho))-v_m/(1+exp(v_m/rho))))*(1-t(x));
f1(x) = A*sinh(v_m/sigma_off)*exp(-(x_off^2/t(x)^2))*exp(1/(1+beta*G(x)*v_m^2))*(1/(1+exp(k*v_m)));
f2(x) = B*sinh(v_m/sigma_on)*exp(-(t(x)^2/x_on^2))*exp(G(x)*v_m^2/sigma_p)*(1/(1+exp(-k*v_m)));
f(x) = f1(x) + f2(x);
end
semilogy(t,f);
xlabel('x','fontsize', 20);
ylabel('y','fontsize', 20);
Could someone help me.
0 comentarios
Respuesta aceptada
Walter Roberson
el 21 de Abr. de 2017
You cannot do much about it. Your expression for f1 includes
exp(1/(1+beta*G(x)*v_m^2))
That value can be arbitrarily high if (1+beta*G(x)*v_m^2) approaches 0; with your beta = 500 and v_m = 1, that is the condition that G(x) approximately equal -1/500 .
You can go to the line above,
G(x) = G_m*t(x)+a*exp(b*sqrt(v_m/(1+exp(-v_m/rho))-v_m/(1+exp(v_m/rho))))*(1-t(x));
and construct
-1/500 == G_m*T+a*exp(b*sqrt(v_m/(1+exp(-v_m/rho))-v_m/(1+exp(v_m/rho))))*(1-T);
and solve for T. You get a result of T about -0.115316250006499 . You are incrementing by 0.001 so when your T becomes -0.115 you have a near singularity .
If your equations are correct then the only way to avoid the singularity is not to calculate near it.
3 comentarios
Walter Roberson
el 21 de Abr. de 2017
Considering that you have 1/(1+beta*G(x)*v_m^2) and the denominator switches between positive and negative, what kind of output were you hoping for?
Más respuestas (0)
Ver también
Categorías
Más información sobre Loops and Conditional Statements en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!