newton shooting method ( 牛顿打靶法） 不收敛

zfhnlg

clear all    %the method of newton for two-points BVP
clc
e1=0.0001;       %迭代精度
e=1;
tm=2;
tspan=[0 tm];     %积分区间
% x0=[0.3236826 0.5779974 1.001114]';
x0=[-0.538 0.288 0.49883]';

while e>e1
y0=[1.076 0 0 x0'];  
[t1,y]=ode45('shooting_li2',tspan,y0);
zb=[y(end,1)-0 y(end,2)-0.576 y(end,3)-0.997661]';     %靶点偏差函数
e=max(abs(zb))

pfpz=[1 0 0 0 0 0;0 1 0 0 0 0;0 0 1 0 0 0];

y1_end=y(end,1);
y2_end=y(end,2);
y3_end=y(end,3);
% y1_end=y(1,1);
% y2_end=y(1,2);
% y3_end=y(1,3);

syms z1 z2 z3 z4 z5 z6
k=1;
A=jacobian([z4;z5;z6;-k*z1./(z1.^2+z2.^2+z3.^2).^(3/2);-k*z2./(z1.^2+z2.^2+z3.^2).^(3/2);-k*z3./(z1.^2+z2.^2+z3.^2).^(3/2)],[z1 z2 z3 z4 z5 z6]);

z1=y1_end;            %代入终点函数值，以获得雅克比数值矩阵
z2=y2_end;
z3=y3_end;
A_zi=subs(A);        %方程雅可比数值矩阵

a14=A_zi(1,4);
a25=A_zi(2,5);
a36=A_zi(3,6);
a41=A_zi(4,1);       %将雅克比数值矩阵与微分方程系数建立关系
a42=A_zi(4,2);
a43=A_zi(4,3);
a51=A_zi(5,1);
a52=A_zi(5,2);
a53=A_zi(5,3);
a61=A_zi(6,1);
a62=A_zi(6,2);
a63=A_zi(6,3);

X0=zeros(6,6);
for i=1:3
t0=[0 tm];
x00=[0 0 0 1 0 0;0 0 0 0 1 0;0 0 0 0 0 1];
[t2 x]=ode45(@(t,x) shooting_li2_newton_jacobi(t,x,a14,a25,a36,a41,a42,a43,a51,a52,a53,a61,a62,a63),t0,x00(i,:));
x11=x(end,1);
x12=x(end,2);
x13=x(end,3);
x14=x(end,4);
x15=x(end,5);
x16=x(end,6);

X0(:,i)=[x11 x12 x13 x14 x15 x16];
end
J=pfpz*X0;
de_x=-inv(J)*zb;
x0=x0+de_x;
while e<e1

    break
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function hs=shooting_li2(t,y)
k=1;
hs(1)=y(4);
hs(2)=y(5);
hs(3)=y(6);
hs(4)=-k*y(1)./(y(1).^2+y(2).^2+y(3).^2).^(3/2);
hs(5)=-k*y(2)./(y(1).^2+y(2).^2+y(3).^2).^(3/2);
hs(6)=-k*y(3)./(y(1).^2+y(2).^2+y(3).^2).^(3/2);
hs=hs(:);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function hs=shooting_li2_newton_jacobi(t,x,a14,a25,a36,a41,a42,a43,a51,a52,a53,a61,a62,a63)

hs(1)=a14*x(4);
hs(2)=a25*x(5);
hs(3)=a36*x(6);
hs(4)=a41*x(1)+a42*x(2)+a43*x(3);
hs(5)=a51*x(1)+a52*x(2)+a53*x(3);
hs(6)=a61*x(1)+a62*x(2)+a63*x(3);
hs=hs(:);
end

zfhnlg

迫切需求专家指导，可能就某一点没有通，但可惜......

zfhnlg

我自己怀疑问题可能出在误差函数（zb ）的雅克比矩阵计算上，请求专家指导。。。。

zfhnlg

这个程序是通的，但是无论改变积分区间、初始靶点，程序都不收敛，这个程序的步骤是按照 1991年 李庆扬 老师 编著的“刚性常微分方程及边值问题”（名字记不清）编写的。

zfhnlg

不好意思，这段程序是可以收敛。经过仔细查看，这段程序没有问题，只要将tm（积分区间上限）改为小一些，程序就收敛了。注：X0改为6X 3矩阵

zfhnlg

将积分区间放小可以收敛，而增大积分区间就不收敛，针对较大的积分区间上限，需要采用continuation technique (连续技术）。

请问，有学长或老师了解的吗？能否指导一下连续技术。谢谢