马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
本帖最后由 煜宸0922 于 2011-5-24 11:35 编辑
各位大侠,这个程序运行时间太长,各位还有什么办法优化一下?
function [T,x,V]=changerk344(N)
x(:,1)=[0;0;0;0]; % 初始状态
h=0.01; % 初始步长,最大步长
AbsTol=1e-5; % 容差
p=4;
T=zeros(1,N); % 时间向量
T(1,1)=0; % 时间起点
p21=2^p-1;
V =[];
tic
for i=1:N-1
if abs(x(1,i))<0.09 % 定步长
h=0.01; % 固定步长
x(:,i+1)=rk144(x(:,i),h,T(i));
T(1,i+1) = T(1,i) + h; % 更新时间
else % 变步长
x1=rk144(x(:,i),h,T(1,i));
x2=rk144(x(:,i),h/2,T(1,i));
a=abs(x1-x2)/p21;
while a(1,1)> AbsTol % 选择合适步长
h=h/10;
x1=rk144(x(:,i),h,T(1,i));
x2=rk144(x(:,i),h/2,T(1,i));
a=abs(x1-x2)/p21;
end
x(:,i+1)=rk144(x(:,i),h,T(i)); % 根据选择的步长更新状态
if abs(x(1,i+1))>=0.1
V=x(3,i+1);
x(3,i+1)=-0.8*x(3,i+1);
x(1,i+1)=sign(x(1,i+1))*min(0.1,abs(x(1,i+1)));
end
T(1,i+1) = T(1,i) + h;
end
h=0.01;
end
toc
T = T(:);
x = x(:,1:N);
|