我还有问题,我的编程思想是这样的:对二阶非线性微分方程通过赋予不同的初值求出方程的解,一旦初始值赋新值就求解,然后将所得的新解与前面的结果作比较,如果结果相同就输出本次初始值。我编了一个程序大家看看能不能优化。我通过profile分析出每次解方程都要费很多时间。
X=zeros(4,10000);
X(1,:)=kron(1:10,ones(1,1000));
X(2,:)=kron(kron(ones(1,10),1:10),ones(1,100));
X(3,:)=kron(kron(ones(1,100),1:10),ones(1,10));
X(4,:)=kron(ones(1,1000),1:10);
refine =10;
G=[];
T=[];
n=[];
Z=[];
Y=[];
y0=[];
for j=1:10000
o=1;
w=1;
k=1;
aa=X(1,j);bb=X(2,j);cc=X(3,j);dd=X(4,j);
while o==k
T(:,o)=[aa;bb;cc;dd] ;
s=aa;
p=bb;
ss=cc;
yy=dd;
y0=[ (s-1/2)*h1+xmin (p-1/2)*h2+ymin (ss-1/2)*h3+zmin (yy-1/2)*h4+umin];
options = odeset('Events',@events,'AbsTol',1e-10,'OutputSel',1,...
'Refine',refine);
for i=1:1
[t,y,te,ye,ie] = ode23(@f,[tstart tfinal],y0,options);%此处费时最多
te
if isempty(ye)==1
break;
end
n=size(ye);
alfa=abs(ye(n(1),1)-ye(n(1),2))/2;%%%%
dy=collision( [ ye(n(1),3) ye(n(1),4) ] );
y0=[ ye(n(1),2) ye(n(1),1) dy(1) dy(2) ];
Z(w:w+3)=y0;
end
if isempty(ye)==1
aa=0;
break;
end
n=size(y0);
aa=round((y0(n(1),1)-xmin)/h1+1/2;
bb=round((y0(n(1),2)-ymin)/h2+1/2);
cc=round((y0(n(1),3)-zmin)/h3+1/2);
dd=round((y0(n(1),4)-umin)/h4+1/2);
o=o+1;
if aa<1|aa>10| bb<1|bb>10|cc<1|cc>10|dd<1|dd>10
aa=0;
else
T(:,o)=[aa;bb;cc;dd] ;
end
if aa==0
break;
else
k=1;
w=1;
G(:,k)=X(:,j);
while T(1,o)~=G(1,k)|T(2,o)~=G(2,k)|T(3,o)~=G(3,k)|T(4,o)~=G(4,k)
k=k+1;
y0=Z(w:w+3);
n=size(y0)
r=round((y0(n(1),1)-xmin)/h1+1/2);%bao ying she
q=round((y0(n(1),2)-ymin)/h2+1/2);%%%%%
rr=round((y0(n(1),3)-zmin)/h3+1/2);
qq=round((y0(n(1),4)-umin)/h4+1/2);
G(:,k)=[r;q;rr;qq];
w=w+4;
end
end
end
if aa~=0
pp=pp+1;
Y(1,pp)=(X(1,j)-1/2)*h1+xmin;
Y(2,pp)=(X(2,j)-1/2)*h2+ymin;
Y(3,pp)=(X(3,j)-1/2)*h3+zmin;
Y(4,pp)=(X(4,j)-1/2)*h4+umin;
Y(:,pp)
end
end |