Valores maximos funcion sinusoidal
Publicado por O (2 intervenciones) el 02/04/2019 21:04:59
Buenas tardes, soy nueva en el mundo de la programación y estoy teniendo una duda existencial en cuanto a mi código, que para ustedes seguro es una chorrada.
Mi intención es realizar un diagrama de bifurcación mediante el metodo runge-kutta 4; para ello debo coger los valores máximos de la función y plotearlos en función de la amplitud, el problema es que no se cómo coger los valores máximos para formar el gráfico adecuadamente...
Aquí os dejo mi código, muchisimas gracias por adelantado:
Mi intención es realizar un diagrama de bifurcación mediante el metodo runge-kutta 4; para ello debo coger los valores máximos de la función y plotearlos en función de la amplitud, el problema es que no se cómo coger los valores máximos para formar el gráfico adecuadamente...
Aquí os dejo mi código, muchisimas gracias por adelantado:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
h= 0.001;
n= 40000;%valores iniciales
x(1)= 0.1;
y(1)= 0.1;
t=[0:h:h*n];
%constantes
a= 0.7;
b= 0.8;
c= 12.5;
T= 1.125;
m=T/h;A=[0.1:0.0001:0.4];
%algoritmo runge-kutta
for z = 1:length(A)
for i= 1:5000 %iteración del valor inicial
kx1= c*(-y(1) + x(1) - (x(1)^3)/3 + A(z)*sin((2*pi*t(1))/T));
ky1= x(1) - b*y(1) + a;
kx2= c*(-(y(1) + ky1*h/2) + (x(1) + kx1*h/2) - ((x(1) + kx1*h/2)^3)/3 +A(z)*sin(2*pi*(t(1)+h/2)/T));
ky2= x(1) + kx1*h/2 - b*(y(1) + ky1*h/2) + a;
kx3= c*(-(y(1) + ky2*h/2) + (x(1) + kx2*h/2) - (x(1) + kx2*h/2)^3/3 +A(z)*sin((2*pi*(t(1)+h/2))/T));
ky3= x(1) + kx2*h/2 - b*(y(1) + ky2*h/2) + a;
kx4= c*(-(y(1) + ky3*h)+ (x(1) + kx3*h) - (x(1) + kx3*h)^3/3 +A(z)*sin((2*pi*(t(1)+h))/T));
ky4= x(1) + kx3*h - b*(y(1) + ky3*h) + a;
x(1)= x(1) + h/6*(kx1 + 2*kx2 + 2*kx3 + kx4);
y(1)= y(1) + h/6*(ky1 + 2*ky2 + 2*ky3 + ky4);
t(1) = t(1) + h;
end
filas=1;for i= 1:n-1 %runge kutta
kx1= c*(-y(i) + x(i) - (x(i)^3)/3 + A(z)*sin((2*pi*t(i))/T));
ky1= x(i) - b*y(i) + a;
kx2=c*(-(y(i) + ky1*h/2) + (x(i) + kx1*h/2) - (x(i) + kx1*h/2)^3/3 +A(z)*sin(2*pi*(t(i)+h/2)/T));
ky2=x(i) + kx1*h/2 - b*(y(i) + ky1*h/2) + a;
kx3=c*(-(y(i) + ky2*h/2) + (x(i) + kx2*h/2) - (x(i) + kx2*h/2)^3/3 +A(z)*sin((2*pi*(t(i)+h/2))/T));
ky3=x(i) + kx2*h/2 - b*(y(i) + ky2*h/2) + a;
kx4=c*(-(y(i) + ky3*h)+ (x(i) + kx3*h) - (x(i) + kx3*h)^3/3 +A(z)*sin((2*pi*(t(i)+h))/T));
ky4=x(i) + kx3*h - b*(y(i) + ky3*h) + a;
x(i+1)= x(i) + h/6*(kx1 + 2*kx2 + 2*kx3 + kx4);
y(i+1)= y(i) + h/6*(ky1 + 2*ky2 + 2*ky3 + ky4);
t(i+1) = t(i) + h;
if(mod(i,m)==0)
XX(filas,z) = max(x(i+1));
filas= filas+1;
end
end
end
plot(A, X, '.', 'markersize', 3);
title('Bifurcation diagram');
xlabel('A'); ylabel('x_n');
set(gca, 'xlim', [0.15 0.25]);
Valora esta pregunta
0