非线性代数方程组零点计算Codim=1的实例[附代码]

mxlzhenzhu

坛子里那么多高手，却不愿意分享自己的代码，遮遮掩掩，小弟初学，愿意分享给大家一个代码，计算非线性代数方程的Fold的实例文件，Predictor-Corrector算法，只能处理Codim=1的情况，Codim=2的情况太难了，还在学习中：

1. %% Predictor-Solver Method verification
   
2. %% mxl.2018-12-6
   
3. %% 总结：本程序非常成功，测试下面的三个例子都没有问题.
   
4. 
   
5. clear
   
6. clc
   
7. %   %---------------------------------------------------------------------%
   
8. %   (1)设计非线性函数
   
9. %   %---------------------------------------------------------------------%
   
10. syms u Lambda
    
11. f=@(u,Lambda)(u-1).^3.*(2*u-5)+u.*Lambda;
    
12. Gu=@(u,Lambda)Lambda + 2*(u - 1)^3 + 3*(2*u - 5)*(u - 1)^2;
    
13. Gy=@(u)u;
    
14. 
    
15. % f=@(u,Lambda)(u.^2-1).*(u.^2-4)+Lambda*u.^2.*exp(0.1*u);
    
16. % 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;
    
17. % Gy=@(u)u.^2.*exp(0.1*u);
    
18. 
    
19. % f=@(u,Lambda)u.^2-4*u+Lambda;
    
20. % Gu=@(u,Lambda)2*u-4;
    
21. % Gy=@(u)1;
    
22. 
    
23. 
    
24. 
    
25. 
    
26. u0=1;Lambda0=0;
    
27. h=1e-3;%计算这个零点
    
28. 
    
29. ds=3e-3;
    
30. Maximum_dS=10*ds;%允许的超曲线最大增量
    
31. Ill_Condition_TOL=1e-5;%小于这个值认为是病态
    
32. TOL=1e-8;%小于这个值认为是收敛的.
    
33. MaxTimes=500;
    
34. direction='Positive';
    
35. N=floor(5/h);
    
36. 
    
37. 
    
38. 
    
39. x=zeros(N,2);%保存
    
40. x(1,:)=[u0,Lambda0];
    
41. Jx=[Gu(u0,Lambda0),Gy(u0)];
    
42. 
    
43. tv0=Compute_Tangent_Vector(Jx,[],direction);
    
44. % Step_Ratio=0.6;%不能太冒进
    
45. tv=[];
    
46. if abs(f(u0,Lambda0))>TOL
    
47.     error('Incorrect I.C.')
    
48. end
    
49. for loop=2:N
    
50.    
    
51.     %----Predictor
    
52.     Jx=[Gu(u0,Lambda0),Gy(u0)];
    
53.     tv=Compute_Tangent_Vector(Jx,tv,direction);
    
54.     if dot(tv,tv0)>0
    
55.         tv0=tv;
    
56.     else
    
57.         tv=-tv;
    
58.     end
    
59.     A=[Jx;tv.'];
    
60.     B=inv(A);%计算第一次
    
61.     dx=B*[-Gy(u0);ds];
    
62.     Step_Ratio=ds/norm(dx);
    
63. %     if Step_Ratio<1
    
64.         dx=Step_Ratio*dx/10;
    
65. %     end
    
66. %     if ~strcmp(direction,'Positive')
    
67. %         up=u0+dx(1:end-1);
    
68. %         Lambda1=Lambda0+dx(end);
    
69. %     else
    
70.         up=u0-dx(1:end-1);
    
71.         Lambda1=Lambda0-dx(end);
    
72. %     end
    
73.    
    
74.     if ~tv(end)
    
75.         disp('It''s a Fold.');
    
76.     end
    
77. 
    
78.     %----Corrector
    
79.     for loopi=1:MaxTimes
    
80.         if loopi>1
    
81.             xk0=[xk1(1);Lambda0];
    
82.             xk1=xk2;            
    
83.         else
    
84.             xk0=[u0;Lambda0];
    
85.             xk1=[up;Lambda1];            
    
86.         end     
    
87.         %----迭代计算xk2，收敛以后保存
    
88.         r=xk1-xk0;
    
89.         delta=[f(xk1(1),xk1(2))-f(xk0(1),xk0(2));r.'*tv-ds];%% 采用 常数h迭代的时候，在达到LP的时候，这里是不正确的，要么修改h，要么想法计算得到LP。
    
90. 
    
91.         temp=dot(r,r);
    
92.         if temp
    
93.            c=(delta-A*r)/temp;
    
94.            A=A+c*r.';
    
95.         else
    
96.            break;
    
97.         end
    
98.         temp=r.'*B*delta;
    
99.         if temp
    
100.            B=B+(r-B*delta)*r.'*B/temp;
     
101.         end
     
102.         v=B*[f(xk1(1),xk1(2));0];
     
103.         xk2=xk1-v;
     
104.         if norm(v)>Maximum_dS
     
105.             xk2=xk1-v;            
     
106.         end
     
107. 
     
108.         if max(norm(v),norm(f(xk2(1),xk2(2))))<TOL           
     
109.             break;
     
110.         end
     
111.       
     
112.     end % for 继续迭代
     
113.    
     
114.     if loopi<MaxTimes
     
115.         x(loop,:)=xk2.';
     
116.         u0=xk2(1:end-1);
     
117.         Lambda0=xk2(end);
     
118.     else
     
119.         warning('Diverged');
     
120.         %----当前迭代步无解，需要新的算法，改变分叉参数增量方向
     
121. %         error('Sub-routine1')
     
122.         if strcmp(direction,'Positvie')   
     
123.             direction='Negative';
     
124.         else
     
125.             direction='Positvie';
     
126.         end
     
127. %         error('Diverged')
     
128.     end
     
129.    
     
130. end
     
131. plot(x(:,end),x(:,1),'b')
     
132. hold on
     
133. scatter(x(:,end),x(:,1),'r')
     
134. hold off

复制代码

上面这个是主程序，再贴一个计算切向向量的子程序：

1. function tv=Compute_Tangent_Vector(Jx,tv,varargin)
   
2. %% 根据Jx矩阵和初始化或者上一步切向向量Tv，计算新的切向向量.
   
3. %% 该程序只适用于亏秩1的矩阵，Example: Jx=[0 0]会失败
   
4. %% 迭代次数很少的。
   
5. [~,n]=size(Jx);% n-1 X n
   
6. 
   
7. if isempty(tv)
   
8.     tv=rand(n,1);
   
9. end
   
10. 
    
11. if nargin==3
    
12.     direction=varargin{1};
    
13. else
    
14.     direction='Positive';
    
15. end
    
16. 
    
17. b=zeros(n,1);
    
18. b(end)=1;
    
19. for loop=1:3
    
20.     A=[Jx;tv.'];
    
21.    
    
22.     tv=A\b;
    
23.     tv=tv./norm(tv);
    
24. end
    
25. 
    
26. if det(A)<0
    
27.         tv=-tv;
    
28. elseif det(A)==0
    
29.     warning('Matrix A is singular.')
    
30. else
    
31.     % nothing, keep det(A) positive
    
32. end
    
33. 
    
34. if ~strcmp(direction,'Positive')
    
35.     tv=-tv;
    
36. end
    
37. 
    
38. end

复制代码

哈哈哈，就这么点吧，希望版主加精。





补充内容 (2019-1-17 21:49):
1）定义了匿名函数f; 2)定义了函数的导数Gu；3）定义了弧长参数导数Gy; 未知数为 u, 增广向量为 x=[u,Lambda]， Lambda是弧长参数。程序功能比较简单，能够跟踪计算非线性方程的零点。

补充内容 (2019-1-17 21:51):
这个Script，要求Gy&和Gu是已知的; 如果不能显式求导，你必需写代码，用有限差分算法获取导数Gu&Gy。

补充内容 (2019-1-17 21:52):
这里并没有包含分叉探测，以及分支转换等代码内容，请自己往上面添加。这里只是简单的Predictor-Corrector算法。

补充内容 (2019-1-17 21:56):
Jx=[Gu,Gy]是一个Jacobian矩阵，n x (n+1)的矩阵，因此有左右零空间向量。计算切向向量的那个函数，可以输入direction='Positive'或者其他，也就是可以计算正反两个方向的切向向量。

jdx_2007

楼主无私奉献精神值得尊重！必须顶！

HLAN

赞！相互学习，共同进步。