16/09/2015
5
Nguyen Quang Hoang
Department of Applied Mechanics
17
>> [t, y] = ode45(@(t,y)vibr1dof(t,y, 1, 0.0,16,1, 4.3), [0:0.01:70], [0,0]);
>> plot(t,y(:,1)), ylabel('y_1 [mm]'), xlabel('t[s]')
-1
-0.5
0
0.5
1
-4
-3
-2
-1
0
1
2
3
4
0
10
20
30
40
50
60
70
-4
-2
0
2
4
y
2
[
m
m
/s
]
t[s]
0
10
20
30
40
50
60
70
-1
-0.5
0
0.5
1
y
1
[
m
m
]
t[s]
Nguyen Quang Hoang
Department of Applied Mechanics
18
-1
-0.5
0
0.5
1
-4
-3
-2
-1
0
1
2
3
4
0
10
20
30
40
50
60
70
-4
-2
0
2
4
y
2
[
m
m
/s
]
t[s]
0
10
20
30
40
50
60
70
-1
-0.
5
0
0.5
1
y
1 [
mm
]
t[s
]
Nguyen Quang Hoang
Department of Applied Mechanics
Giải hệ phương trình vi phân cấp 1
19
2
,
dx
dy
x
y
x
xy
dt
dt
(0)
0, (0) 1
x
y
2
1
2
1
2
1
2
1
1 2
,
,
dx
dx
x
x x
y
x
x
x
x x
dt
dt
Trước hết xét việc giải hệ hai phương trình vi phân cấp một sau
điều kiện đầu
Để giải hệ trên bằng ode45, ta viết một m-file thể hiện vế phải của hệ trên
function
xdot = eqx(t,x);
xdot = zeros(2,1);
xdot(1) = -x(1)^2 + x(2);
xdot(2) = -x(1) - x(1)* x(2);
end
% save with file name eqx.m
Viết lại hệ dạng véc tơ như sau:
Nguyen Quang Hoang
Department of Applied Mechanics
Giải hệ phương trình vi phân cấp 1
20
>> [t,x] = ode45('eqx',[0 10],[0,1]);
Nghiệm cần tìm x được ghi lại bằng hai véctơ cột, cột thứ nhất x1 = x(:,1)
và cột thứ hai x2 = x(:,2). Với lệnh vẽ
>> plot(t,x(:,1),t,x(:,2),'--'),xlabel('t'),
>> axis([0 10
–1.12 1.12])
0
2
4
6
8
10
- 1
- 0.5
0
0.5
1
t
- 0.6
- 0.4
- 0.2
0
0.2
0.4
0.6
- 0.4
- 0.2
0
0.2
0.4
0.6
0.8
1
>> plot(x(:,1),x(:,2), 'k-'),
