homotopy.m

求解非线性方程组 的同伦曲线法

%文件名：Homotopy.m
function root = Homotopy(funcF,funcH,x,N,tol,Nmax)
% 功能: 求解非线性方程组的同伦算法
% 输入:
% funcF----- 原始方程函数句柄
% funcH----- 同伦方程函数句柄
% x       ---- 初始迭代点
% N(可选) ---- 同伦参数t的划分数（默认是10）
% tol(可选)---精度要求（默认是1e-4） 
% Nmax(可选)---最大迭代次数（默认100次）
% 输出:
% root--- -解向量 
if size(x,1) == 1; x = x'; end  
if nargin<6;Nmax=100;end
if nargin<5 | isempty(tol); tol=1e-4;end
if nargin<4 | isempty(N);N=10;end
k=0;
x0=x;
n=length(x);
dx=tol+1;
f0=tol+1;
for m=1:N
    t=m/N;
while (norm(dx)>tol & norm(f0)>tol)&  k<Nmax
    [jac,f0] = JacobianH(x,funcF,funcH,x0,t);
    dx = jac\(-f0);
    x = x + dx;  
    k=k+1;
end
dx=tol+1;
f0=tol+1;
root=x;
k=0;
end
 
function [jac,f0] = JacobianH(x,func1,func2,x0,t)
%计算同伦函数的Jacobian矩阵和函数值
h = 1.0e-6;
n = length(x);
jac = zeros(n);
f0 = feval(func2,x,func1,x0,t);
for i =1:n
    temp = x(i);
    x(i) = temp + h;
    f1 = feval(func2,x,func1,x0,t);
    x(i) = temp;
    jac(:,i) = (f1 - f0)/h;
end