Click the blue text above to subscribe!
In the vast field of scientific and engineering computation, the problem of solving nonlinear equation systems shines like a brilliant pearl, illuminating the combination of mathematical wisdom and computational technology. From simulating multi-body systems in physics to analyzing market equilibrium in economics, from parameter optimization in machine learning to stability studies in control systems, nonlinear equation systems are ubiquitous, and the quality of their solutions directly affects the reliability and accuracy of the entire system.
The general mathematical form of a nonlinear equation system can be expressed as finding a vector such that , where is a nonlinear vector-valued function. Specifically, it can be expanded as:
The essential challenge of such problems lies in their nonlinear characteristics, which may lead to the existence of multiple solutions, a lack of uniqueness, or extreme sensitivity to initial conditions. Unlike linear equation systems, nonlinear problems typically do not have analytical solutions and must rely on numerical iterative methods for approximate solutions. This numerical solving process requires a solid mathematical theoretical foundation as well as exquisite programming implementation skills.
Basic Newton’s Method: A Classic and Efficient Starting Point
The Newton-Raphson method, as the most classic numerical method, has a solid theoretical foundation and satisfactory convergence speed. Its core iterative formula is:
where is the Jacobian matrix of the function at point , containing all first-order partial derivative information.
Let us understand the implementation of Newton’s method through a specific example. Consider solving the equation system:
#include<iostream>
#include<cmath>
#include<Eigen/Dense>
using namespace Eigen;
using namespace std;
struct NewtonSystem {
VectorXd operator()(const VectorXd& x) {
VectorXd f(2);
double x1 = x(0), x2 = x(1);
f(0) = x1*x1 + x2*x2 - 4.0;
f(1) = exp(x1) + x2 - 1.0;
return f;
}
MatrixXd jacobian(const VectorXd& x) {
MatrixXd J(2, 2);
double x1 = x(0);
J(0, 0) = 2*x1;
J(0, 1) = 2*x(1);
J(1, 0) = exp(x1);
J(1, 1) = 1.0;
return J;
}
};
VectorXd newtonSolver(const VectorXd& initialGuess, double tol=1e-6, int max_iter=100) {
NewtonSystem system;
VectorXd x = initialGuess;
for (int i = 0; i < max_iter; ++i) {
VectorXd f = system(x);
double norm = f.norm();
cout << "Iteration " << i << ": norm = " << norm << endl;
if (norm < tol) {
cout << "Converged after " << i << " iterations" << endl;
return x;
}
MatrixXd J = system.jacobian(x);
x -= J.fullPivLu().solve(f);
}
cout << "Warning: Did not converge within max iterations" << endl;
return x;
}
int main() {
VectorXd initialGuess(2);
initialGuess << 1.0, 1.0;
VectorXd solution = newtonSolver(initialGuess);
cout << "Solution: x = " << solution(0) << ", y = " << solution(1) << endl;
NewtonSystem system;
VectorXd residual = system(solution);
cout << "Residual: " << residual.norm() << endl;
return 0;
}This implementation demonstrates the core idea of Newton’s method: approximating the solution to a nonlinear problem through local linearization. The Jacobian matrix provides a local linear approximation of the function at the current point, while the matrix inversion operation is equivalent to solving a linear system to determine the next search direction.
Finite Difference Method: A Practical Choice When Analytical Derivatives Are Unavailable
However, real-world problems are often more complex. Consider a system like this:
The analytical expression of the Jacobian matrix for this system is extremely complex, making manual derivation both difficult and error-prone. In this case, the finite difference method provides an elegant solution.
#include<iostream>
#include<cmath>
#include<functional>
#include<Eigen/Dense>
using namespace Eigen;
using namespace std;
class FiniteDifferenceSolver {
public:
using VectorFunction = function<VectorXd(const VectorXd&)>;
FiniteDifferenceSolver(VectorFunction func, double tol=1e-6, int max_iter=100)
: func(func), tol(tol), max_iter(max_iter) {}
MatrixXd numericalJacobian(const VectorXd& x, double h=1e-6) {
int n = x.size();
VectorXd f0 = func(x);
int m = f0.size();
MatrixXd J(m, n);
for (int j = 0; j < n; ++j) {
VectorXd x_plus = x;
x_plus(j) += h;
VectorXd f_plus = func(x_plus);
J.col(j) = (f_plus - f0) / h;
}
return J;
}
VectorXd solve(const VectorXd& initialGuess) {
VectorXd x = initialGuess;
for (int iter = 0; iter < max_iter; ++iter) {
VectorXd f = func(x);
double norm = f.norm();
cout << "Iteration " << iter << ": norm = " << norm << endl;
if (norm < tol) {
cout << "Converged after " << iter << " iterations" << endl;
return x;
}
MatrixXd J = numericalJacobian(x);
x -= J.colPivHouseholderQr().solve(f);
}
cout << "Warning: Did not converge within max iterations" << endl;
return x;
}
private:
VectorFunction func;
double tol;
int max_iter;
};
VectorXd complexSystem(const VectorXd& x) {
VectorXd f(2);
double x1 = x(0), x2 = x(1);
f(0) = sin(x1*x2) + exp(x1 - x2) - 1.5;
f(1) = log(abs(x1) + 1.0) + x1*x2*x2 - 2.0;
return f;
}
int main() {
FiniteDifferenceSolver solver(complexSystem);
VectorXd initialGuess(2);
initialGuess << 1.0, 1.0;
VectorXd solution = solver.solve(initialGuess);
cout << "Solution: x = " << solution(0) << ", y = " << solution(1) << endl;
VectorXd residual = complexSystem(solution);
cout << "Residual: " << residual.norm() << endl;
return 0;
}The core idea of the finite difference method is based on the definition of derivatives:
By perturbing slightly, we probe the rate of change of the function. This method, while computationally intensive, has excellent versatility.
Broyden’s Method: A Model of Intelligent Update Strategy
When computational efficiency becomes a critical consideration, Broyden’s method demonstrates its unique advantages. This method does not require recalculating the entire Jacobian matrix at each iteration but instead improves the approximate matrix step by step through an intelligent update strategy.
#include<iostream>
#include<cmath>
#include<functional>
#include<Eigen/Dense>
using namespace Eigen;
using namespace std;
class BroydenSolver {
public:
using VectorFunction = function<VectorXd(const VectorXd&)>;
BroydenSolver(VectorFunction func, double tol=1e-6, int max_iter=100)
: func(func), tol(tol), max_iter(max_iter) {}
MatrixXd numericalJacobian(const VectorXd& x, double h=1e-6) {
int n = x.size();
VectorXd f0 = func(x);
int m = f0.size();
MatrixXd J(m, n);
for (int j = 0; j < n; ++j) {
VectorXd x_plus = x;
x_plus(j) += h;
VectorXd f_plus = func(x_plus);
J.col(j) = (f_plus - f0) / h;
}
return J;
}
VectorXd solve(const VectorXd& initialGuess) {
int n = initialGuess.size();
VectorXd x = initialGuess;
VectorXd f = func(x);
MatrixXd B = numericalJacobian(x);
for (int iter = 0; iter < max_iter; ++iter) {
double norm = f.norm();
cout << "Iteration " << iter << ": norm = " << norm << endl;
if (norm < tol) {
cout << "Converged after " << iter << " iterations" << endl;
return x;
}
VectorXd s = B.colPivHouseholderQr().solve(-f);
VectorXd x_new = x + s;
VectorXd f_new = func(x_new);
VectorXd y = f_new - f;
if (s.norm() > 1e-10) {
B += (y - B*s) * s.transpose() / s.squaredNorm();
}
x = x_new;
f = f_new;
}
cout << "Warning: Did not converge within max iterations" << endl;
return x;
}
private:
VectorFunction func;
double tol;
int max_iter;
};
VectorXd testSystem(const VectorXd& x) {
VectorXd f(2);
double x1 = x(0), x2 = x(1);
f(0) = x1*x1 - 2*x1 + x2*x2 - x2 - 8;
f(1) = x1*x1 - 4*x1 + x2*x2 + x2 - 10;
return f;
}
int main() {
BroydenSolver solver(testSystem);
VectorXd initialGuess(2);
initialGuess << 3.0, 2.0;
VectorXd solution = solver.solve(initialGuess);
cout << "Solution: x = " << solution(0) << ", y = " << solution(1) << endl;
VectorXd residual = testSystem(solution);
cout << "Residual: " << residual.norm() << endl;
return 0;
}The update formula of Broyden’s method reflects the intelligent learning mechanism in numerical optimization, improving the understanding of system behavior based on information from the current iteration.
Nelder-Mead Method: A Robust Choice for Derivative-Free Optimization
When derivative information is completely unavailable or the function is noisy, the Nelder-Mead simplex method provides a robust approach based entirely on function value comparisons.
#include<iostream>
#include<cmath>
#include<functional>
#include<vector>
#include<algorithm>
#include<Eigen/Dense>
using namespace Eigen;
using namespace std;
// Disable character encoding warnings from Eigen library
#pragma warning(disable: 4819)
class NelderMeadSolver {
public:
using ScalarFunction = function<double(const VectorXd&)>;
NelderMeadSolver(ScalarFunction func, double tol=1e-6, int max_iter=1000)
: func(func), tol(tol), max_iter(max_iter) {}
VectorXd solve(const VectorXd& initialGuess, double step=0.1) {
// Use Eigen::Index instead of int to avoid type conversion warnings
Eigen::Index n = initialGuess.size();
vector<VectorXd> simplex;
simplex.push_back(initialGuess);
for (Eigen::Index i = 0; i < n; ++i) {
VectorXd point = initialGuess;
point(i) += step;
simplex.push_back(point);
}
for (int iter = 0; iter < max_iter; ++iter) {
vector<pair<double, VectorXd>> evaluated;
for (auto& point : simplex) {
evaluated.emplace_back(func(point), point);
}
// Use traditional comparison function instead of lambda expression
sort(evaluated.begin(), evaluated.end(),
[](const pair<double, VectorXd>& a, const pair<double, VectorXd>& b) {
return a.first < b.first;
});
double range = evaluated.back().first - evaluated[0].first;
if (range < tol) {
cout << "Converged after " << iter << " iterations" << endl;
return evaluated[0].second;
}
VectorXd best = evaluated[0].second;
VectorXd worst = evaluated.back().second;
VectorXd centroid = VectorXd::Zero(n);
for (Eigen::Index i = 0; i < n; ++i) {
centroid += evaluated[i].second;
}
centroid /= n;
VectorXd reflected = centroid + 1.0*(centroid - worst);
double f_reflected = func(reflected);
if (f_reflected < evaluated[0].first) {
VectorXd expanded = centroid + 2.0*(reflected - centroid);
double f_expanded = func(expanded);
if (f_expanded < f_reflected) {
evaluated.back().second = expanded;
} else {
evaluated.back().second = reflected;
}
} else if (f_reflected < evaluated[n-1].first) {
evaluated.back().second = reflected;
} else {
VectorXd contracted = centroid + 0.5*(worst - centroid);
double f_contracted = func(contracted);
if (f_contracted < evaluated.back().first) {
evaluated.back().second = contracted;
} else {
for (size_t i = 1; i < evaluated.size(); ++i) {
evaluated[i].second = best + 0.5*(evaluated[i].second - best);
evaluated[i].first = func(evaluated[i].second);
}
}
}
simplex.clear();
// Replace structured bindings with traditional writing
for (auto& elem : evaluated) {
simplex.push_back(elem.second);
}
}
cout << "Warning: Did not converge within max iterations" << endl;
return simplex[0];
}
private:
ScalarFunction func;
double tol;
int max_iter;
};
double systemNorm(const VectorXd& x) {
VectorXd f(2);
double x1 = x(0), x2 = x(1);
f(0) = x1*x1 + x2*x2 - 4.0;
f(1) = exp(x1) + x2 - 1.0;
return f.squaredNorm();
}
int main() {
NelderMeadSolver solver(systemNorm);
VectorXd initialGuess(2);
initialGuess << 1.0, 1.0;
VectorXd solution = solver.solve(initialGuess);
cout << "Solution: x = " << solution(0) << ", y = " << solution(1) << endl;
double residual = systemNorm(solution);
cout << "Residual norm: " << residual << endl;
return 0;
}The Nelder-Mead method explores the search space through geometric operations such as reflection, expansion, and contraction. Although its convergence speed is slower, it has minimal requirements on the properties of the function and exhibits strong adaptability.
Comprehensive Considerations in Practical Applications
In practical engineering applications, the choice of method needs to consider multiple factors. The scale of the problem, the required precision, the characteristics of the function (smoothness, computational cost, etc.), and the available computational resources all influence the choice of method. Hybrid strategies often provide better robustness—for example, starting with the more robust Nelder-Mead method to obtain a good initial approximation, then switching to the more efficient Newton’s method or Broyden’s method for fine-tuning the solution.
The choice of initial values is crucial for all iterative methods. A reasonable initial guess can significantly improve convergence speed and even avoid getting trapped in local optima. Multi-start strategies, physical insights, problem analysis, or global search methods can be used to obtain better initial points.
Convergence criteria also need to be handled with care. In addition to checking whether the residual norm is less than the tolerance, it is also necessary to consider the change in the solution and the upper limit of iterations to avoid infinite loops or pseudo-convergence.
The journey of solving nonlinear equation systems is both challenging and valuable. From the classic Newton’s method to modern intelligent algorithms, each method has its unique advantages and applicable scenarios. Understanding the mathematical characteristics of the problem, choosing the appropriate algorithm, and implementing and debugging carefully are key to successfully solving nonlinear equation systems. With the continuous development of computational technology, we can expect to solve increasingly complex nonlinear problems more effectively, providing more powerful computational tools for scientific research and engineering applications.
• end • 
Companionship is the longest confession of love
We push the most practical information for you
Scan the QR code to follow us