有限元 弧长法 matlab,有没有有关弧长法的程序！！！

% % % This is a code of Pseudo-Arc Length Continuation Method

%  the method can be used for solving the nonlinear equations

%  the principle can re found in some text books

clc

clear

% ------------------------开始预测--------------

k=1;

for ds=0.02:0.02:9

a0=1;

x0=3;

y0=2.6458;

a=a0;

x=x0;

y=y0;

step=5;

for i=1:1

DF(1,1)=2*a*x;

DF(1,2)=2*y;

DF(1,3)=x^2;

DF(2,1)=2*x;

DF(2,2)=-2*y;

DF(2,3)=0;

DF;

J1=DF(:,2:end);

J2=DF(:,1);

K=DF(:,3);

J2=[J2,K];

J3=DF(:,1:2);

J1=((-1)^(1+1))*det(J1);

J2=((-1)^(2+1))*det(J2);

J3=((-1)^(3+1))*det(J3);

vv=[J1,J2,J3]';

v=vv/norm(vv);

end

% --------------预测结束-------------------

% -------------开始校正----------------------

x=3;

y=3;

a=2;

p=[x,y,a]';

F(1,1)=a*x^2+y^2-16;

F(1,2)= x^2-y^2-2;

F(1,3)=(x-x0)*v(1,1)+(y-y0)*v(2,1)+(a-a0)*v(3,1)-ds;

delta=1e-6;

epsilon=1e-6;

for i=1:100

J=zeros(3,3);

J(1,1)=2*a*x; J(1,2)=2*y; J(1,3)=x^2;

J(2,1)=2*x; J(2,2)=-2*y;  J(2,3)=0;

J(3,1)=v(1,1); J(3,2)=v(2,1);  J(3,3)=v(3,1);

q=p-(J\F');

x=q(1);

y=q(2);

a=q(3);

F(1,1)=a*x^2+y^2-16;

F(1,2)= x^2-y^2-2;

F(1,3)=(x-x0)*v(1,1)+(y-y0)*v(2,1)+(a-a0)*v(3,1)-ds;

F;

Z=F;

err=norm(q-p);

relerr=err/(norm(q)+eps);

p=q;

Y=Z;

if (err

q

break

end

end

q

ds;

x0=q(1);

y0=q(2);

a0=q(3);

c(k)=a0;

b(k)=x0;

k=k+1;

a=q(3,1);

end

b

c

plot(c,b)

hold on

% a=c

% x=sqrt(18./(a+1))

% a=2:-0.1:-1

% y=sqrt(18./(a+1)-2)

% plot(a,x,'o')

% hold on

%

a=-1:0.1:1

x=sqrt(18./(a+1))

plot(a,x,'o')

hold on

我找到了这个程序，不知道对不对