具有雅可比稀疏模式的大型非线性方程组
此示例说明如何提供雅可比稀疏模式来有效地求解大型稀疏方程组。该示例使用为 n 方程组定义的目标函数,
要解的方程是 、。该示例使用 。本示例末尾的 nlsf1a 辅助函数实现了目标函数 。
在求解同一系统的示例 具有雅可比矩阵的大型稀疏非线性方程组 中,目标函数具有明确的雅可比矩阵。然而,有时您无法明确计算雅可比矩阵,但您可以确定雅可比矩阵中哪些元素是非零的。在这个示例中,雅可比矩阵是三对角的,因为在 的定义中出现的唯一变量是 、 和 。因此,雅可比矩阵的唯一非零部分位于主对角线及其两个相邻的对角线。雅可比稀疏模式是一个矩阵,其非零元素对应于雅可比矩阵中(潜在的)非零元素。
创建一个稀疏的 × 三对角矩阵,表示雅可比稀疏模式。(有关此代码的解释,请参阅 spdiags。)
n = 1000; e = ones(n,1); Jstr = spdiags([e e e],-1:1,n,n);
查看 Jstr 的结构。
spy(Jstr)

设置 fsolve 选项以使用 'trust-region' 算法,这是唯一可以使用雅可比稀疏模式的 fsolve 算法。
options = optimoptions('fsolve','Algorithm','trust-region','JacobPattern',Jstr);
将初始点 x0 设置为-1 个条目的向量。
x0 = -1*ones(n,1);
求解。
[x,fval,exitflag,output] = fsolve(@nlsf1a,x0,options);
Equation solved. fsolve completed because the vector of function values is near zero as measured by the value of the function tolerance, and the problem appears regular as measured by the gradient. <stopping criteria details>
检查得到的残差范数和 fsolve 所采用的函数计算的次数。
disp(norm(fval))
1.0523e-09
disp(output.funcCount)
25
残差范数接近于零,表明 fsolve 正确地求解了该问题。函数计算的数量相当少,只有大约二十几个。将此函数计算的数量与未提供雅可比稀疏模式的数量进行比较。
options = optimoptions('fsolve','Algorithm','trust-region'); [x,fval,exitflag,output] = fsolve(@nlsf1a,x0,options);
Equation solved. fsolve completed because the vector of function values is near zero as measured by the value of the function tolerance, and the problem appears regular as measured by the gradient. <stopping criteria details>
disp(norm(fval))
1.0657e-09
disp(output.funcCount)
5005
求解器基本上可以得出与之前相同的解,但需要进行数千次函数计算才能实现。
辅助函数
以下代码会创建 nlsf1a 辅助函数。
function F = nlsf1a(x) % Evaluate the vector function n = length(x); F = zeros(n,1); i = 2:(n-1); F(i) = (3-2*x(i)).*x(i)-x(i-1)-2*x(i+1) + 1; F(n) = (3-2*x(n)).*x(n)-x(n-1) + 1; F(1) = (3-2*x(1)).*x(1)-2*x(2) + 1; end