zfhnlg 发表于 2013-4-26 16:07
急求........
以下是三自由度两点边值问题的拟牛顿迭代打耙法,这个算法是传统的,因为只具有一个初始条件;对于两个初始条件的两点边值问题,大家遇到过吗?
%%%%%%%%%%%%%%%%%%%%主程序
clear all %采用Fx(差值)作为精度依据
clc
tspan=0:0.1:2;
x0=[-0.538 0.288 0.49883]'; %定义初始耙点
x=x0;
norm(x);
y0=[1.076 0 0 x']'; %定义初始条件
A0=eye(3);
e1=0.0000001;
[t,y]=ode45('shooting_li2',tspan,y0);
Fx0=[y(end,1)-0;y(end,2)-0.576;y(end,3)-0.997661];
%s=-inv(A)*Fx; % 由此可以看出,非点乘,点乘则要求矩阵的行列是一样的
c0=norm(Fx0,1);
c=c0;
%Y=zeros();
A=A0; %定义初始迭代矩阵(单位矩阵)
Fx=Fx0; %定义初始精度条件
while c>e1
A;
s=-inv(A)*Fx;
x=x+s; %耙点迭代
%Y=[Y' c]';
y0=[1.076 0 0 x']';
[t,y]=ode45('shooting_li2',tspan,y0); %解初始问题
Fx2=[y(end,1)-0;y(end,2)-0.576;y(end,3)-0.997661]; %精度条件
y(:,1);
c=norm(Fx2,1);
y=Fx2-Fx;
A=A+(y-A*s)*s'/(s'*s); %矩阵迭代
Fx=Fx2;
%C=[t y(end,1) y(end,2) y(end,3)];
while c<=e1
[t,y]=ode45('shooting_li2',tspan,y0); %解初始问题
y1=y(:,1);
%y2=y(2)
%y3=y(3)
C=[t y(:,1) y(:,2) y(:,3)] %输入结果
break %防止不断循环
end
end
%%%%%%%%%%%%%%%%%%%%%%%三自由度、6维函数
function hs=shooting_li2(t,y)
hs(1)=y(4); %y1'=y3
hs(2)=y(5); %y2'=y4
hs(3)=y(6);
hs(4)=-y(1)./(y(4).^2+y(5).^2+y(6).^2).^(3/2); %y3'=
hs(5)=-y(2)./(y(4).^2+y(5).^2+y(6).^2).^(3/2); %y4'=
hs(6)=-y(3)./(y(4).^2+y(5).^2+y(6).^2).^(3/2);
hs=hs(:);
end |