matlab求解非线性方程组

非线性方程组求解
1.mulStablePoint用不动点迭代法求非线性方程组的一个根
function [r,n]=mulStablePoint(F,x0,eps)
%非线性方程组：f
%初始解：a
%解的精度：eps
%求得的一组解：r
%迭代步数：n
if nargin==2
eps=1.0e-6;
end
x0 = transpose(x0);
n=1;
tol=1;
while tol>eps
r= subs(F,findsym(F),x0); %迭代公式
tol=norm(r-x0); %注意矩阵的误差求法，norm为矩阵的欧几里德范数
n=n+1;
x0=r;
if(n>100000) %迭代步数控制
disp('迭代步数太多，可能不收敛！');
return;
end
end
2.mulNewton用牛顿法法求非线性方程组的一个根
function [r,n]=mulNewton(F,x0,eps)
if nargin==2
eps=1.0e-4;
end
x0 = transpose(x0);
Fx = subs(F,findsym(F),x0);
var = findsym(F);
dF = Jacobian(F,var);
dFx = subs(dF,findsym(dF),x0);
r=x0-inv(dFx)*Fx;
n=1;
tol=1;
while tol>eps
x0=r;
Fx = subs(F,findsym(F),x0);
dFx = subs(dF,findsym(dF),x0);
r=x0-inv(dFx)*Fx; %核心迭代公式
tol=norm(r-x0);
n=n+1;
if(n>100000) %迭代步数控制 disp('迭代步数太多，可能不收敛！');
return;
end
end
3.mulDiscNewton用离散牛顿法法求非线性方程组的一个根
function [r,m]=mulDiscNewton(F,x0,h,eps)
format long;
if nargin==3
eps=1.0e-8;
end
n = length(x0);
fx = subs(F,findsym(F),x0);
J = zeros(n,n);
for i=1:n
x1 = x0;
x1(i) = x1(i)+h(i);
J(:,i) = (subs(F,findsym(F),x1)-fx)/h(i);
end
r=transpose(x0)-inv(J)*fx;
m=1;
tol=1;
while tol>eps
xs=r;
fx = subs(F,findsym(F),xs);
J = zeros(n,n);
for i=1:n
x1 = xs;
x1(i) = x1(i)+h(i);
J(:,i) = (subs(F,findsym(F),x1)-fx)/h(i);
end
r=xs-inv(J)*fx; %核心迭代公式
tol=norm(r-xs);
m=m+1;
if(m>100000) %迭代步数控制 disp('迭代步数太多，可能不收敛！');
return;
end
end
format short;
4.mulMix用牛顿-雅可比迭代法求非线性方程组的一个根function [r,m]=mulMix(F,x0,h,l,eps)
if nargin==4
eps=1.0e-4;
end
n = length(x0);
J = zeros(n,n);
Fx = subs(F,findsym(F),x0);
for i=1:n
x1 = x0;
x1(i) = x1(i)+h(i);
J(:,i) = (subs(F,findsym(F),x1)-Fx)/h(i);
end
D = diag(diag(J));
C =
D - J;
inD = inv(D);
H = inD*C;
Hm = eye(n,n);
for i=1:l-1
Hm = Hm + power(H,i);
end
dr = Hm*inD*Fx;
r = transpose(x0)-dr; m=1;
tol=1;
while tol>eps
x0=r;
Fx = subs(F,findsym(F),x0);
J = zeros(n,n);
for i=1:n
x1 = x0;
x1(i) = x1(i)+h(i);
J(:,i) = (subs(F,findsym(F),x1)-Fx)/h(i);
end
D = diag(diag(J));
C =
D - J;
inD = inv(D);
H = inD*C;
Hm = eye(n,n);
for i=1:l-1
Hm = Hm + power(H,i);
end
dr = Hm*inD*Fx;
r = x0-dr; %核心迭代公式
tol=norm(r-x0);
m=m+1;
if(m>100000) %迭代步数控制 disp('迭代步数太多，可能不收敛！');
return;
end
end
5.mulNewtonSOR用牛顿-SOR迭代法求非线性方程组的一个根
function [r,m]=mulNewtonSOR(F,x0,w,h,l,eps)
if nargin==5
eps=1.0e-4;
end
n = length(x0);
J = zeros(n,n);
Fx = subs(F,findsym(F),x0);
for i=1:n
x1 = x0;
x1(i) = x1(i)+h(i);
J(:,i) = (subs(F,findsym(F),x1)-Fx)/h(i);
end
D = diag(diag(J));
L = -tril(J-D);
U = -triu(J-D);
inD = inv(D-w*L);
H = inD*(D - w*D+w*L);;
Hm = eye(n,n);
for i=1:l-1
Hm = Hm + power(H,i);
end
dr = w*Hm*inD*Fx;
r = transpose(x0)-dr;
m=1;
tol=1;
while tol>eps
x0=r;
Fx = subs(F,findsym(F),x0);
J = zeros(n,n);
for i=1:n
x1 = x0;
x1(i) = x1(i)+h(i);
J(:,i) = (subs(F,findsym(F),x1)-Fx)/h(i);
end
D = diag(diag(J));
L = -tril(J-D);
U = -triu(J-D);
inD = inv(D-w*L);
H = inD*(D - w*D+w*L);;
Hm = eye(n,n);
for i=1:l-1
Hm = Hm + power(H,i);
end
dr = w*Hm*inD*Fx;
r = x0-dr; %核心迭代公式
tol=norm(r-x0);
m=m+1;
if(m>100000) %迭代步数控制 disp('迭代步数太多，可能不收敛！');
return;
end
end
6.mulDNewton用牛顿下山法求非线性方程组的一个根
function [r,m]=mulDNewton(F,x0,eps)
%非线性方程组：F
%初始解：x0
%解的精度：eps
%求得的一组解：r
%迭代步数：n
if nargin==2
eps=1.0e-4;
end
x0 = transpose(x0);
dF = Jacobian(F);
m=1;
tol=1;
while tol>eps
ttol=1;
w=1;
Fx = subs(F,findsym(F),x0);
dFx = subs(dF,findsym(dF),x0);
F1=norm(Fx);
while ttol>=0 %下面的循环是选取下山因子w的过程
r=x0-w*inv(dFx)*Fx; %核心的迭代公式
Fr = subs(F,findsym(F),r);
ttol=norm(Fr)-F1;
w=w/2;
end
tol=norm(r-x0);
m=m+1;
x0=r;
if(m>100000) %迭代步数控制
disp('迭代步数太多，可能不收敛！');
return;
end
end
7.mulGXF1用两点割线法的第一种形式求非线性方程组的一个根
function [r,m]=mulGXF1(F,x0,x1,eps)
format long;
if nargin==3
eps=1.0e-4;
end
x0 = transpose(x0);
x1 = transpose(x1);
n = length(x0);
fx = subs(F,findsym(F),x0);
fx1 = subs(F,findsym(F),x1);
h = x0 - x1;
J = zeros(n,n);
for i=1:n
xt = x1;
xt(i) = x0(i);
J(:,i) = (subs(F,findsym(F),xt)-fx1)/h(i);
end
r=x1-inv(J)*fx1;
m=1;
tol=1;
while tol>eps
x0 = x1;
x1 = r;
fx = subs(F,findsym(F),x0);
fx1 = subs(F,findsym(F),x1);
h = x0 - x1;
J = zeros(n,n);
for i=1:n
xt = x1;
xt(i) = x0(i);
J(:,i) = (subs(F,findsym(F),xt)-fx1)/h(i);
end
r=x1-inv(J)*fx1;
tol=norm(r-x1);
m=m+1;
if(m>100000) %迭代步数控制 disp('迭代步数太多，可能不收敛！');
return;
end
end
format short;
8.mulGXF2用两点割线法的第二种形式求非线性方程组的一个根
function [r,m]=mulGXF2(F,x0,x1,eps)
format long;
if nargin==3
eps=1.0e-4;
end
x0 = transpose(x0);
x1 = transpose(x1);
n = length(x0);
fx = subs(F,findsym(F),x0);
fx1 = subs(F,findsym(F),x1);
h = x0 - x1;
J = zeros(n,n);
xt = x1;
xt(1) = x0(1);
J(:,1) = (subs(F,findsym(F),xt)-subs(F,findsym(F),x1))/h(1);
for i=2:n
xt = x1;
xt(1:i) = x0(1:i);
xt_m = x1;
xt_m(1:i-1) = x0(1:i-1);
J(:,i) = (subs(F,findsym(F),xt)-subs(F,findsym(F),xt_m))/h(i);
end
r=x1-inv(J)*fx1;
m=1;
tol=1;
while tol>eps
x0 = x1;
x1 = r;
fx = subs(F,findsym(F),x0);
fx1 = subs(F,findsym(F),x1);
h = x0 - x1;
J = zeros(n,n);
xt = x1;
xt(1) = x0(1);
J(:,1) = (subs(F,findsym(F),xt)-subs(F,findsym(F),x1))/h(1);
for i=2:n
xt = x1;
xt(1:i) = x0(1:i);
xt_m = x1;
xt_m(1:i-1) = x0(1:i-1);
J(:,i) = (subs(F,findsym(F),xt)-subs(F,findsym(F),xt_m))/h(i);
end
r=x1-inv(J)*fx1;
tol=norm(r-x1);
m=m+1;
if(m>100000) %迭代步数控制 disp('迭代步数太多，可能不收敛！');
return;
end
end
format short;
9.mulVNewton用拟牛顿法求非线性方程组的一组解
function [r,m]=mulVNewton(F,x0,A,eps)
%方程组：F
%方程组的初始解：x0
% 初始A矩阵：A
%解的精度：eps
%求得的一组解：r
%迭代步数：m
if nargin==2
A=eye(length(x0)); %A取为单位阵
eps=1.0e-4;
else
if nargin==3
eps=1.0e-4;
end
end
x0 = transpose(x0);
Fx = subs(F, findsym(F),x0);
r=x0-A\Fx;
m=1;
tol=1;
while tol>eps
x0=r;
Fx = subs(F, findsym(F),x0);
r=x0-A\Fx;
y=r-x0;
Fr = subs(F, findsym(F),r);
z= Fr-Fx;
A1=A+(z-A*y)*transpose(y)/norm(y); %调整A A=A1;
m=m+1;
if(m>100000) %迭代步数控制 disp('迭代步数太多，可能不收敛！');
return;
end
tol=norm(r-x0);
end
10.mulRank1用对称秩1算法求非线性方程组的一个根
function [r,n]=mulRank1(F,x0,A,eps)
if nargin==2
l = length(x0);
A=eye(l); %A取为单位阵
eps=1.0e-4;
else
if nargin==3
eps=1.0e-4;
end
end
fx = subs(F,findsym(F),x0);
r=transpose(x0)-inv(A)*fx;
n=1;
tol=1;
while tol>eps
x0=r;
fx = subs(F,findsym(F),x0);
r=x0-inv(A)*fx;
y=r-x0;
fr = subs(F,findsym(F),r);
z = fr-fx;
A1=A+ fr *transpose(fr)/(transpose(fr)*y); %调整A A=A1;
n=n+1;
if(n>100000) %迭代步数控制
disp('迭代步数太多，可能不收敛！');
return;
end
tol=norm(r-x0);
end
11.mulDFP用D-F-P算法求非线性方程组的一组解
function [r,n]=mulDFP(F,x0,A,eps)
if nargin==2
l = length(x0);
B=eye(l); %A取为单位阵
eps=1.0e-4;
else
if nargin==3
eps=1.0e-4;
end
end
fx = subs(F,findsym(F),x0);
r=transpose(x0)-B*fx;
n=1;
tol=1;
while tol>eps
x0=r;
fx = subs(F,findsym(F),x0);
r=x0-B*fx;
y=r-x0;
fr = subs(F,findsym(F),r);
z = fr-fx;
B1=B+ y*y'/(y'*z)-B*z*z'*B/(z'*B*z); %调整A
B=B1;
n=n+1;
if(n>100000) %迭代步数控制
disp('迭代步数太多，可能不收敛！');
return;
end
tol=norm(r-x0);
end
12.mulBFS用B-F-S算法求非线性方程组的一个根
function [r,n]=mulBFS(F,x0,B,eps)
if nargin==2
l = length(x0);
B=eye(l); %B取为单位阵
eps=1.0e-4;
else
if nargin==3
eps=1.0e-4;
end
end
fx = subs(F,findsym(F),x0);
r=transpose(x0)-B*fx;
n=1;
tol=1;
while tol>eps
x0=r;
fx = subs(F,findsym(F),x0);
r=x0-B*fx;
y=r-x0;
fr = subs(F,findsym(F),r);
z = fr-fx;
u = 1 + z'*B*z/(y'*z);
B1= B+ (u*y*y'-B*z*y'-y*z'*B)/(y'*z); %调整B B=B1;
n=n+1;
if(n>100000) %迭代步数控制
disp('迭代步数太多，可能不收敛！');
return;
end
tol=norm(r-x0);
end
13.mulNumYT用数值延拓法求非线性方程组的一组解
function [r,m]=mulNumYT(F,x0,h,N,eps)
format long;
if nargin==4
eps=1.0e-8;
end
n = length(x0);
fx0 = subs(F,findsym(F),x0);
x0 = transpose(x0);
J = zeros(n,n);
for k=0:N-1
fx = subs(F,findsym(F),x0);
for i=1:n
x1 = x0;
x1(i) = x1(i)+h(i);
J(:,i) = (subs(F,findsym(F),x1)-fx)/h(i);
end
inJ = inv(J);
r=x0-inJ*(fx-(1-k/N)*fx0);
x0 = r;
end
m=1;
tol=1;
while tol>eps
xs=r;
fx = subs(F,findsym(F),xs);
J = zeros(n,n);
for i=1:n
x1 = xs;
x1(i) = x1(i)+h(i);
J(:,i) = (subs(F,findsym(F),x1)-fx)/h(i);
end
r=xs-inv(J)*fx; %核心迭代公式
tol=norm(r-xs);
m=m+1;
if(m>100000) %迭代步数控制 disp('迭代步数太多，可能不收敛！');
return;
end
end
format short;
14.DiffParam1用参数微分法中的欧拉法求非线性方程组的一组解
function r=DiffParam1(F,x0,h,N)
%非线性方程组：f
%初始解：x0
%数值微分增量步大小：h
%雅可比迭代参量:l
%解的精度：eps
%求得的一组解：r
%迭代步数：n
x0 = transpose(x0);
n = length(x0);
ht = 1/N;
Fx0 = subs(F,findsym(F),x0);
for k=1:N
Fx = subs(F,findsym(F),x0);
J = zeros(n,n);
for i=1:n
x1 = x0;
x1(i) = x1(i)+h(i);
J(:,i) = (subs(F,findsym(F),x1)-Fx)/h(i);
end
inJ = inv(J);
r = x0 - ht*inJ*Fx0;
x0 = r;
end
15.DiffParam2用参数微分法中的中点积分法求非线性方程组的一组解function r=DiffParam2(F,x0,h,N)
%非线性方程组：f
%初始解：x0
%数值微分增量步大小：h
%雅可比迭代参量:l
%解的精度：eps
%求得的一组解：r
%迭代步数：n
x0 = transpose(x0);
n = length(x0);
ht = 1/N;
Fx0 = subs(F,findsym(F),x0);
J = zeros(n,n);
for i=1:n
xt = x0;
xt(i) = xt(i)+h(i);
J(:,i) = (subs(F,findsym(F),xt)-Fx0)/h(i);
end
inJ = inv(J);
x1 = x0 - ht*inJ*Fx0;
for k=1:N
x2 = x1 + (x1-x0)/2;
Fx2 = subs(F,findsym(F),x2);
J = zeros(n,n);
for i=1:n
xt = x2;
xt(i) = xt(i)+h(i);
J(:,i) = (subs(F,findsym(F),xt)-Fx2)/h(i);
end
inJ = inv(J);
r = x1 - ht*inJ*Fx0;
x0 = x1;
x1 = r;
end
16.mulFastDown用最速下降法求非线性方程组的一组解
function [r,m]=mulFastDown(F,x0,h,eps)
format long;
if nargin==3
eps=1.0e-8;
end
n = length(x0);
x0 = transpose(x0);
m=1;
tol=1;
while tol>eps
fx = subs(F,findsym(F),x0);
J = zeros(n,n);
for i=1:n
x1 = x0;
x1(i) = x1(i)+h;
J(:,i) = (subs(F,findsym(F),x1)-fx)/h;
end
lamda = fx/sum(diag(transpose(J)*J));
r=x0-J*lamda; %核心迭代公式
fr = subs(F,findsym(F),r);
tol=dot(fr,fr);
x0 = r;
m=m+1;
if(m>100000) %迭代步数控制 disp('迭代步数太多，可能不收敛！');
return;
end
end
format short;
17.mulGSND用高斯牛顿法求非线性方程组的一组解
function [r,m]=mulGSND(F,x0,h,eps)
format long;
if nargin==3
eps=1.0e-8;
end
n = length(x0);
x0 = transpose(x0);
m=1;
tol=1;
while tol>eps
fx = subs(F,findsym(F),x0);
J = zeros(n,n);
for i=1:n
x1 = x0;
x1(i) = x1(i)+h;
J(:,i) = (subs(F,findsym(F),x1)-fx)/h;
end
DF = inv(transpose(J)*J)*transpose(J);
r=x0-DF*fx; %核心迭代公式
tol=norm(r-x0);
x0 = r;
m=m+1;
if(m>100000) %迭代步数控制 disp('迭代步数太多，可能不收敛！');
return;
end
end
format short;
18.mulConj用共轭梯度法求非线性方程组的一组解
function [r,m]=mulConj(F,x0,h,eps)
format long;
if nargin==3
eps=1.0e-6;
end
n = length(x0);
x0 = transpose(x0);
fx0 = subs(F,findsym(F),x0);
p0 = zeros(n,n);
for i=1:n
x1 = x0;
x1(i) = x1(i)*(1+h);
p0(:,i) = -(subs(F,findsym(F),x1)-fx0)/h;
end
m=1;
tol=1;
while tol>eps
fx = subs(F,findsym(F),x0);
J = zeros(n,n);
for i=1:n
x1 = x0;
x1(i) = x1(i)+h;
J(:,i) = (subs(F,findsym(F),x1)-fx)/h;
end
lamda = fx/sum(diag(transpose(J)*J));
r=x0+p0*lamda; %核心迭代公式
fr = subs(F,findsym(F),r);
Jnext = zeros(n,n);
for i=1:n
x1 = r;
x1(i) = x1(i)+h;
Jnext(:,i) = (subs(F,findsym(F),x1)-fr)/h;
end
abs1 = transpose(Jnext)*Jnext;
abs2 = transpose(J)*J;
v = abs1/abs2;
if (abs(det(v)) < 1)
p1 = -Jnext+p0*v;
else
p1 = -Jnext;
end
tol=norm(r-x0);
p0 = p1;
x0 = r;
m=m+1;
if(m>100000) %迭代步数控制 disp('迭代步数太多，可能不收敛！');
return;
end
end
format short;
19.mulDamp用阻尼最小二乘法求非线性方程组的一组解
function [r,m]=mulDamp(F,x0,h,u,v,eps)
format long;
if nargin==5
eps=1.0e-6;
end
FI = transpose(F)*F/2;
n = length(x0);
x0 = transpose(x0);
m=1;
tol=1;
while tol>eps
j = 0;
fx = subs(F,findsym(F),x0);
J = zeros(n,n);
for i=1:n
x1 = x0;
x1(i) = x1(i)+h;
afx = subs(F,findsym(F),x1);
J(:,i) = (afx-fx)/h;
end
FIx = subs(FI,findsym(FI),x0);
for i=1:n
x2 = x0;
x2(i) = x2(i)+h;
gradFI(i,1) = (subs(FI,findsym(FI),x2)-FIx)/h;
end
s=0;
while s==0
A = transpose(J)*J+u*eye(n,n);
p = -A\gradFI;
r = x0 + p;
FIr = subs(FI,findsym(FI),r);
if FIr<FIx
if j == 0
u = u/v;
j = 1;
else
s=1;
end
else
u = u*v;
j = 1;
if norm(r-x0)<eps
s=1;
end
end
end
x0 = r;
tol = norm(p);
m=m+1;
if(m>100000) %迭代步数控制disp('迭代步数太多，可能不收敛！');
return;
end
end
format short;。