function Test_ode45
%calculate the motion
in 20 secs
tspan
= [0 20];
%x
=4, v
=0 when t
=0
x0
= [4 0];
% [t
, x
] = ode45
(@hamofun
, tspan
, x0
);
% plot
(t
,x
);
% hold on
[t
, x
] = ode45
(@dampfun
, tspan
, x0
);
plot
(t
,x
);
end
function xd
=hamofun
(t
, x
)
%Example
: hamonic motion
% ddot
(x
)+w
^2*x
=0
%---------------------------------------------
%xd
(n
): n th order derivative of x
xd
=zeros
(2,1);
%x
(1)-x
; x
(2)-v
: depend on
input
xd
(1)=x
(2);
xd
(2)=-2^2*x
(1);
end
function xd
=dampfun
(t
, x
)
%Example
: damped vibration
% ddot
(x
)=-w
^2*x
-a
*v
...
%---------------------------------------------
a
=0.4;
xd
=zeros
(2,1);
xd
(1)=x
(2);
xd
(2)=-2^2*x
(1)-a
*x
(2);
end
备注
解了两个方程，第一个是谐振子，第二个是受迫振动，均为二阶常微分方程。解之前设置初值，这是必须的。 我自己在写微分方程时，变量替换的时候有些晕，知道其物理意义可能理解起来快一些：代码中的x是两列的矩阵，第一列x(1)是位移，第二列x(2)是速度；位移的一阶导xd(1)是速度，二阶导xd(2)是加速度。 plot(t,x)时自然画出两条曲线，一条是位移随时间变化曲线，另一条是速度。
