非线性代数方程组零点计算Codim=1的实例[附代码]
本帖最后由 mxlzhenzhu 于 2019-1-7 11:06 编辑坛子里那么多高手，却不愿意分享自己的代码，遮遮掩掩，小弟初学，愿意分享给大家一个代码，计算非线性代数方程的Fold的实例文件，Predictor-Corrector算法，只能处理Codim=1的情况，Codim=2的情况太难了，还在学习中：
%% Predictor-Solver Method verification
%% mxl.2018-12-6
%% 总结：本程序非常成功，测试下面的三个例子都没有问题.
clear
clc
% %---------------------------------------------------------------------%
% (1)设计非线性函数
% %---------------------------------------------------------------------%
syms u Lambda
f=@(u,Lambda)(u-1).^3.*(2*u-5)+u.*Lambda;
Gu=@(u,Lambda)Lambda + 2*(u - 1)^3 + 3*(2*u - 5)*(u - 1)^2;
Gy=@(u)u;
% f=@(u,Lambda)(u.^2-1).*(u.^2-4)+Lambda*u.^2.*exp(0.1*u);
% Gu=@(u,Lambda)2*u*(u^2 - 1) + 2*u*(u^2 - 4) + 2*Lambda*u*exp(u/10) + (Lambda*u^2*exp(u/10))/10;
% Gy=@(u)u.^2.*exp(0.1*u);
% f=@(u,Lambda)u.^2-4*u+Lambda;
% Gu=@(u,Lambda)2*u-4;
% Gy=@(u)1;
u0=1;Lambda0=0;
h=1e-3;%计算这个零点
ds=3e-3;
Maximum_dS=10*ds;%允许的超曲线最大增量
Ill_Condition_TOL=1e-5;%小于这个值认为是病态
TOL=1e-8;%小于这个值认为是收敛的.
MaxTimes=500;
direction='Positive';
N=floor(5/h);
x=zeros(N,2);%保存
x(1,:)=;
Jx=;
tv0=Compute_Tangent_Vector(Jx,[],direction);
% Step_Ratio=0.6;%不能太冒进
tv=[];
if abs(f(u0,Lambda0))>TOL
error('Incorrect I.C.')
end
for loop=2:N
%----Predictor
Jx=;
tv=Compute_Tangent_Vector(Jx,tv,direction);
if dot(tv,tv0)>0
tv0=tv;
else
tv=-tv;
end
A=;
B=inv(A);%计算第一次
dx=B*[-Gy(u0);ds];
Step_Ratio=ds/norm(dx);
% if Step_Ratio<1
dx=Step_Ratio*dx/10;
% end
% if ~strcmp(direction,'Positive')
% up=u0+dx(1:end-1);
% Lambda1=Lambda0+dx(end);
% else
up=u0-dx(1:end-1);
Lambda1=Lambda0-dx(end);
% end
if ~tv(end)
disp('It''s a Fold.');
end
%----Corrector
for loopi=1:MaxTimes
if loopi>1
xk0=;
xk1=xk2;
else
xk0=;
xk1=;
end
%----迭代计算xk2，收敛以后保存
r=xk1-xk0;
delta=;%% 采用 常数h迭代的时候，在达到LP的时候，这里是不正确的，要么修改h，要么想法计算得到LP。
temp=dot(r,r);
if temp
c=(delta-A*r)/temp;
A=A+c*r.';
else
break;
end
temp=r.'*B*delta;
if temp
B=B+(r-B*delta)*r.'*B/temp;
end
v=B*;
xk2=xk1-v;
if norm(v)>Maximum_dS
xk2=xk1-v;
end
if max(norm(v),norm(f(xk2(1),xk2(2))))<TOL
break;
end
end % for 继续迭代
if loopi<MaxTimes
x(loop,:)=xk2.';
u0=xk2(1:end-1);
Lambda0=xk2(end);
else
warning('Diverged');
%----当前迭代步无解，需要新的算法，改变分叉参数增量方向
% error('Sub-routine1')
if strcmp(direction,'Positvie')
direction='Negative';
else
direction='Positvie';
end
% error('Diverged')
end
end
plot(x(:,end),x(:,1),'b')
hold on
scatter(x(:,end),x(:,1),'r')
hold off上面这个是主程序，再贴一个计算切向向量的子程序：
function tv=Compute_Tangent_Vector(Jx,tv,varargin)
%% 根据Jx矩阵和初始化或者上一步切向向量Tv，计算新的切向向量.
%% 该程序只适用于亏秩1的矩阵，Example: Jx=会失败
%% 迭代次数很少的。
[~,n]=size(Jx);% n-1 X n
if isempty(tv)
tv=rand(n,1);
end
if nargin==3
direction=varargin{1};
else
direction='Positive';
end
b=zeros(n,1);
b(end)=1;
for loop=1:3
A=;
tv=A\b;
tv=tv./norm(tv);
end
if det(A)<0
tv=-tv;
elseif det(A)==0
warning('Matrix A is singular.')
else
% nothing, keep det(A) positive
end
if ~strcmp(direction,'Positive')
tv=-tv;
end
end
哈哈哈，就这么点吧，希望版主加精{:4_74:}。
补充内容 (2019-1-17 21:49):
1）定义了匿名函数f; 2)定义了函数的导数Gu；3）定义了弧长参数导数Gy; 未知数为 u, 增广向量为 x=， Lambda是弧长参数。程序功能比较简单，能够跟踪计算非线性方程的零点。
补充内容 (2019-1-17 21:51):
这个Script，要求Gy&和Gu是已知的; 如果不能显式求导，你必需写代码，用有限差分算法获取导数Gu&Gy。
补充内容 (2019-1-17 21:52):
这里并没有包含分叉探测，以及分支转换等代码内容，请自己往上面添加。这里只是简单的Predictor-Corrector算法。
补充内容 (2019-1-17 21:56):
Jx=是一个Jacobian矩阵，n x (n+1)的矩阵，因此有左右零空间向量。计算切向向量的那个函数，可以输入direction='Positive'或者其他，也就是可以计算正反两个方向的切向向量。 楼主无私奉献精神值得尊重！必须顶！ 赞！相互学习，共同进步。
页:
[1]