newton shooting method ( 牛顿打靶法） 不收敛
clear all %the method of newton for two-points BVPclc
e1=0.0001; %迭代精度
e=1;
tm=2;
tspan=; %积分区间
% x0=';
x0=[-0.538 0.288 0.49883]';
while e>e1
y0=;
=ode45('shooting_li2',tspan,y0);
zb='; %靶点偏差函数
e=max(abs(zb))
pfpz=;
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(,);
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=;
x00=;
=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)=;
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
迫切需求专家指导，可能就某一点没有通，但可惜...... 我自己怀疑问题可能出在误差函数（zb ）的雅克比矩阵计算上，请求专家指导。。。。 这个程序是通的，但是无论改变积分区间、初始靶点，程序都不收敛，这个程序的步骤是按照 1991年 李庆扬 老师 编著的“刚性常微分方程及边值问题”（名字记不清）编写的。 不好意思，这段程序是可以收敛。经过仔细查看，这段程序没有问题，只要将tm（积分区间上限）改为小一些，程序就收敛了。注：X0改为6X 3矩阵 将积分区间放小可以收敛，而增大积分区间就不收敛，针对较大的积分区间上限，需要采用continuation technique (连续技术）。
请问，有学长或老师了解的吗？能否指导一下连续技术。谢谢
页:
[1]