同伦（Homotopy）算法是求解非线性方程组 F(x)=0 的一种强大且全局收敛的数值方法。它通过构造一个从简单问题 G(x)=0 到目标问题 F(x)=0 的连续形变路径，并沿着这条路径追踪解，从而有效地避开牛顿法等传统局部方法对初始值敏感的缺点。
核心思想与算法流程
其核心思想可以直观地理解为：假设你想从已知的'起点' G(x)=0 的解 x₀，走到未知的'终点' F(x)=0 的解 x*。同伦算法就是构造一条连接两者的平滑路径（同伦路径），然后小心翼翼地沿着这条路径走到终点。
最常见的同伦是凸同伦： H(x, t) = t * F(x) + (1 - t) * G(x) = 0，其中 t ∈ [0, 1] 是参数。
- 当 t=0 时，H(x, 0) = G(x) = 0，对应简单问题。
- 当 t=1 时，H(x, 1) = F(x) = 0，对应目标问题。
算法的主要步骤如下：
- 构造同伦方程：选择简单的初始函数 G(x)（例如 G(x)=x - x₀）和同伦参数 t。
- 设置起点：取 t=0，求解（或已知）简单方程 G(x)=0 的解 x₀。点 (x₀, 0) 是路径起点。
- 预测步：在当前点 (x_k, t_k)，沿路径切线方向移动一小步，得到预测点 (x̃, t̃)。
- 校正步：固定 t=t̃，以 x̃ 为初值，用牛顿法求解方程 H(x, t̃)=0，得到校正后的精确路径点 (x_{k+1}, t_{k+1})。
- 参数步进：重复预测 - 校正过程，使 t 从 0 逐渐增加至 1。当 t=1 时，对应的 x 即为目标方程 F(x)=0 的数值解。
具体流程逻辑如下：
- 开始：构造同伦方程 H(x,t)
- 起点：t=0，求解 x₀ 使得 H(x₀,0)=0
- 循环判断：是否满足 t ≥ 1？
- 否：执行预测步（沿切线方向计算 x̃, t̃）→ 校正步（固定 t=t̃，用牛顿法求解 H(x, t̃)=0）→ 更新路径点 → 继续循环
- 是：输出目标解 x* ≈ x_{k+1}，结束
MATLAB 代码实现示例
使用预测 - 校正（Parameter Continuation）同伦方法求解非线性方程组的 MATLAB 示例。
function[x_solution, iter, solution_path]=homotopy_solver(F, x0, varargin)
% 使用预测 - 校正同伦法求解非线性方程组 F(x) = 0
% 输入：
% F - 函数句柄，返回方程值（列向量）和雅可比矩阵 [F_val, J] = F(x)
% x0 - 同伦路径起点，应满足 G(x)=x-x0=0，即简单问题的解
% varargin - 可选参数：'MaxIter', 'Tol', 'StepSize'
% 输出：
% x_solution - F(x)=0 的近似解
% iter - 实际迭代次数
% solution_path - 追踪的路径 (t, x) 历史记录
% 设置默认参数
p = inputParser;
addParameter(p,'MaxIter',100,@isnumeric);
addParameter(p,'Tol',1e-10,@isnumeric);
addParameter(p,'StepSize',0.05,@isnumeric);
parse(p, varargin{:});
max_iter = p.Results.MaxIter;
tol = p.Results.Tol;
dt = p.Results.StepSize;
n = length(x0);
x = x0;
% 初始解
t = 0;
% 初始同伦参数
path = [t, x'];
% 记录路径
fprintf('开始同伦路径追踪...\n');
fprintf('t 从 0 -> 1, 步长 dt = %.3f\n', dt);
% --- 主循环：追踪同伦路径直到 t >= 1 ---
for iter = 1:max_iter
% 1. 预测步（欧拉法沿切线方向预测）
[H_val, J_H] = compute_homotopy(F, x, t);
% 计算关于 (x, t) 的全导数：J_H * dx/dt + dH/dt = 0
dH_dt = F(x);
% 凸同伦 H = t*F + (1-t)*(x-x0)，所以 dH/dt = F(x) - (x-x0)
% 求解切线方向 [dx/dt; 1]
tangent = -J_H \ dH_dt;
dx_dt = tangent(1:n);
% 只取关于 x 的部分
% 预测下一步
t_pred = min(t + dt, 1.0);
% 确保 t 不超过 1
x_pred = x + dx_dt * (t_pred - t);
% 2. 校正步（固定 t，用牛顿法求解 H(x, t_pred)=0）
[x_corr, newton_success] = newton_corrector(F, x_pred, t_pred, x0, tol);
if ~newton_success
warning('在校正步中牛顿法未收敛，尝试减小步长。');
dt = dt * 0.5;
% 步长减半
if dt < 1e-4
error('步长过小，路径追踪失败。');
end
continue;
% 不更新迭代计数，重试当前步
end
% 更新当前点和参数
x = x_corr;
t = t_pred;
path = [path; t, x'];
% 记录新点
% 打印进度信息
residual = norm(F(x));
fprintf('迭代 %3d: t = %.4f, 残差 = %.4e\n', iter, t, residual);
% 检查是否到达终点且满足精度
if t >= 1.0 - eps && residual < tol
fprintf('\n同伦路径追踪成功完成！\n');
fprintf('在 t = %.6f 处达到终点，最终残差 = %.4e\n', t, residual);
break;
end
% 自适应调整步长（可选：可根据校正步的迭代次数调整 dt）
end
if iter == max_iter && (t < 1 || residual > tol)
warning('达到最大迭代次数，但可能未完全收敛。');
end
x_solution = x;
solution_path = path;
end
function[H_val, J_H] = compute_homotopy(F, x, t, x0)
% 计算凸同伦 H(x,t) = t*F(x) + (1-t)*(x-x0) 的值和雅可比矩阵
if nargin < 4
x0 = zeros(size(x));
% 默认 G(x) = x
end
[F_val, J_F] = F(x);
% 用户函数需返回雅可比矩阵
H_val = t * F_val + (1-t)*(x - x0);
J_H = t * J_F + (1-t)*eye(length(x));
end
function[x_corr, success] = newton_corrector(F, x_init, t, x0, tol)
% 校正步：用牛顿法求解 H(x, t) = 0
max_newton_iter = 20;
x_corr = x_init;
success = false;
for n_iter = 1:max_newton_iter
[H_val, J_H] = compute_homotopy(F, x_corr, t, x0);
dx = -J_H \ H_val;
x_corr = x_corr + dx;
if norm(dx) < tol || norm(H_val) < tol
success = true;
break;
end
end
end
