-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathOcode-showspring2d.m
101 lines (88 loc) · 2.36 KB
/
Ocode-showspring2d.m
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
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
k1 = 1;
k2 = 100;
l = 1;
m = 1;
Y(1) = 0.1;
Vx(1) = 0;
X(1) = 0.2;
Vy(1) = 0;
dt = 0.01;
i = 1;
for t = 1:100
a= 4*( X(i)^2 + Y(i)^2);
z1 = 0.5*sqrt((l^2)+4*(l*X(i))+ a);
l1 = z1-l/2;
z2 = 0.5*sqrt((l^2)-4*(l*X(i))+a);
l2 = z2-l/2;
b1= l/2 + X(i);
b2 = l/2 - X(i);
dvy1 = -(k1*l1*Y(i)/z1 + k2*l2*Y(i)/z2)*dt/m;
dvx1 = -((k1*l1*b1/z1) - (k2*l2*b2/z2))*dt/m;
dx1 = Vx(i)*dt;
dy1 = Vy(i)*dt;
z1 = 0.5*sqrt((l^2)+4*(l*(X(i)+dx1/2))+ a);
l1 = z1-l/2;
z2 = 0.5*sqrt((l^2)-4*(l*(X(i)+dx1/2))+a);
l2 = z2-l/2;
b1= l/2 + X(i)+dx1/2;
b2 = l/2 - X(i) - dx1/2;
dvy2 = -(k1*l1*(Y(i)+dy1/2)/z1 + k2*l2*(Y(i)+dy1/2)/z2)*dt/m;
dvx2 = -((k1*l1*b1/z1) - (k2*l2*b2/z2))*dt/m;
dx2 = (Vx(i)+dvx1/2)*dt;
dy2 = (Vy(i)+dvy1/2)*dt;
z1 = 0.5*sqrt((l^2)+4*(l*(X(i)+dx2/2))+ a);
l1 = z1-l/2;
z2 = 0.5*sqrt((l^2)-4*(l*(X(i)+dx2/2))+a);
l2 = z2-l/2;
b1= l/2 + X(i)+dx2/2;
b2 = l/2 - X(i) - dx2/2;
dvy3 = -(k1*l1*(Y(i)+dy2/2)/z1 + k2*l2*(Y(i)+dy2/2)/z2)*dt/m;
dvx3 = -((k1*l1*b1/z1) - (k2*l2*b2/z2))*dt/m;
dx3 = (Vx(i)+dvx2/2)*dt;
dy3 = (Vy(i)+dvy2/2)*dt;
z1 = 0.5*sqrt((l^2)+4*(l*(X(i)+dx3))+ a);
l1 = z1-l/2;
z2 = 0.5*sqrt((l^2)-4*(l*(X(i)+dx3))+a);
l2 = z2-l/2;
b1= l/2 + X(i)+dx3;
b2 = l/2 - X(i) - dx3;
dvy4 = -(k1*l1*(Y(i)+dy3)/z1 + k2*l2*(Y(i)+dy3)/z2)*dt/m;
dvx4 = -((k1*l1*b1/z1) - (k2*l2*b2/z2))*dt/m;
dx4 = (Vx(i)+dvx3)*dt;
dy4 = (Vy(i)+dvy3)*dt;
% z = 0.5*sqrt((l^2)+4*((Y(i)+dy1/2)^2));
% x = z-l/2;
% dy2 = (V(i)+dv1/2)*dt;
% dv2 = -((k1+k2)/m)*(x*(Y(i)+dy1/2)/z)*dt;
%
% z = 0.5*sqrt((l^2)+4*((Y(i)+dy2)^2));
% x = z-l/2;
% dy3 = (V(i)+dv2)*dt;
% dv3 = -((k1+k2)/m)*(x*(Y(i)+dy2)/z)*dt;
dy = (dy1+2*dy2+2*dy3+dy4)/6;
dvx = (dvx1+2*dvx2+2*dvx3+dvx4)/6;
dvy = (dvy1+2*dvy2+2*dvy3+dvy4)/6;
dx= (dx1+2*dx2+2*dx3+dx4)/6;
Y(i+1) = Y(i)+dy;
X(i+1) = X(i)+dx;
Vx(i+1) = Vx(i)+dvx;
Vy(i+1) = Vy(i)+dvy;
#V(i+1) = V(i) +dv
plot(X(i+1), Y(i+1));
hold on
for j = -l/2:0.01:X(i+1)
plot(j, (Y(i+1)/(X(i+1)+l/2))*(j+l/2), 'r')
hold on
end
for j = X(i+1):0.01:l/2
plot(j, (Y(i+1)/(X(i+1)-l/2))*(j-l/2), 'g')
hold on
end
# plot([(-l/2):X(i+1)],[0:Y(i+1)]);
# hold on
#clear
pause(0.01);
#clear;
hold off
i = i+1;
end