## PJFNK Matlab

Version 1.0.0 (37.2 KB) by
This repository is for an implementation of the Preconditioned Jacobian-Free Newton-Krylov method in Matlab.

Updated 19 Mar 2023

From GitHub

# PJFNK MATLAB

This repository is for an implementation of the Preconditioned Jacobian-Free Newton-Krylov method in Matlab, as well as methods for standard full Newton method and formulation of full and sparse Numerical Jacobian matrices. With this method, users can choose a Krylov method from any of the following:

• 'gmres' - Generalized Minimum RESidual

• (default) Most flexible method, but may require restarts to be practical
• 'bicgstab' - BiConjugate Gradients Stabilized

• Enhanced BiCG method with a GMRES(1) step for stabilization
• 'bicgstabl'

• Enhanced BiCGStab method with a GMRES(2) step for stabilization
• 'pcg'

• Optimal method for Symmetric Positive Definite Jacobians
• 'minres'

• Jacobians must be Symmetric, but does not need to be Positive Definite
• 'symmlq'

• Jacobians must be Symmetric, but does not need to be Positive Definite
• 'cgs'

• Similar to BiCG, but does not require transposition of matrix
• 'tfqmr'

• Optimal for Symmetric Jacobians that are indefinite

# Requirements

• MATLAB
• [Optional] Optimization Toolbox

NOTE: The Optimization Toolbox in MATLAB is not required to run these codes, but is used for validation checks in some of the unit tests.

# Basic Usage

To use this function, you must provide (at a minimum) a:

• Function Handle to evaluate a set of equations (F = @(x) ...)

• An initial guess vector (x0)

The size of the vector (x0) and the size of the evaluated function handle (F(x)) must both be Nx1.

``````% Example 1 - Basics

F = @(x) [2*x(1)*x(2) - x(2)^2;
-x(1)+2*x(2)-exp(-x(3));
-x(2)+2*x(3)^3];

x0 = [1; 1; 1];

[x,stats,opts] = JacobianFreeNewtonKrylov(F,x0);

``````

The function will return the solution to the system (x) and 2 struct() objects containing: (i) the status information of the solve (stats) and (ii) a set of options used in the solve (opts). Several default options are provided.

You can also provide a struct() to the input arguments of the function to change the behavior of the solve as needed.

``````% Example 2 - Solver Options
options = struct();
options.maxiter = 20;  % Set the maximum number of non-linear steps
options.ftol = 1e-6;   % Set the solver tolerance
options.krylov_solver = 'gmres';        % Pick from a selection of Krylov methods
options.linesearch = 'backtracking';    % Pick from a selection of linesearch methods

% [Optional] Provide a function handle for a known Jacobian
options.jacfun = J = @(x)    [2*x(2), 2*x(1)-2*x(2), 0;
-1, 2, exp(-x(3));
0, -1, 6*x(3)^2];
``````

This is a Jacobian-Free method, so you technically never need to supply the Jacobian. However, you can actually choose to use the full (or sparse) Jacobian (if you desire).

``````options.use_matrix = true;   % Specifies to use the Jacobian for the solve
``````

By default, the method will set `use_matrix` to false, and just use the `JacobianOperator` function in conjunction with a Krylov solver. However, if you tell the solver to use the matrix form, it will do that instead.

For the Krylov solvers, you can also control how those methods are solved by adding to the `options` data structure.

``````% Example 3 - Setting Krylov options for 'gmres'
options.krylov_opts = struct();
options.krylov_opts.restart = 50;  % Iterations before restart
options.krylov_opts.tol = 1e-6;    % Linear solver tolerance
``````

NOTE: Each Krylov solver will have its own set of options. The options available (including preconditioners) follows exactly with the Matlab standard set of options for each iterative solver (see here).

Lastly, there is also a full `NewtonMethod` implementation. You can use this method if you do not want to use Krylov solvers underneath the non-linear problem. In this case, the search direction is resolved using the Matlab `mldivide` () function. This method also has almost all the same options available to it as the `JacobianFreeNewtonKrylov` method.

``````% Example 4 - Using full Newton Method
F = @(x) [2*x(1) - x(2);
-x(1)+2*x(2)-x(3);
-x(2)+2*x(3) + 1];

x0 = [1; 0; 0];

% Exact Jacobian
J = @(x) [2, -1, 0;
-1, 2, -1;
0, -1, 2];

% Store the Jacobian function handle in the data structure
options = struct();
options.jacfunc = J;

[x,stats,opts] = NewtonMethod(F,x0,options);
``````

For more usage examples, see the `tests` folder of this repository.

# Functions

• `JacobianOperator`

Performs an action of a Jacobian on any vector v without forming the Jacobian explicitly.

``````% Pseudo Code
Jv = (F(x+epsilon*v) - F(x))/epsilon;
``````

where F(x) evaluates the non-linear system at state x.

where v is the vector being multiplied by the Jacobian.

where epsilon is a small perturbation value.

• `NumericalJacobianMatrix`

Approximates the full Jacobian with finite-differences approach. Can return either a 'dense' or 'sparse' Jacobian via user request. Default is to return a 'dense' Jacobian.

``````% Pseudo Code
for i=1:N
dx = x;
dx(i,1) =  x(i,1) + epsilon;
J(:, i) = (F(dx) - F(x))/epsilon;
end
``````

where F(x) evaluates the non-linear system at state x.

where epsilon is a small perturbation value.

• `NewtonMethod`

Performs a basic Newton's method with direct inversion of the Jacobian matrix. The user may provide an exact Jacobian function as a function handle. If they do not, then the method will use the `NumericalJacobianMatrix` function to approximate the Jacobian.

``````% Pseudo Code
while (norm(F(x)) > ftol)
% Perform a Newton Step
%   J*s = -F --> s = -J\F
x = x0 - J(x)\F(x);

% Update x
x0=x;
end
``````

where F(x) evaluates the non-linear system at state x.

where J(x) evaluates the Jacobian matrix at state x.

where ftol is the error or residual tolerance for convergence.

• `JacobianFreeNewtonKrylov`

Performs a Jacobian-Free Newton-Krylov solve with associated user options. The user may choose any valid Krylov method to solve the Newton step.

``````% Pseudo Code
while (norm(F(x)) > ftol)
% Perform a Newton Step
%   J*s = -F --> s = -J\F
x = x0 + KrylovSolve( @JacOp( @F(x0),x0,F), -F );

% Update x
x0=x;
end
``````

where F(x) evaluates the non-linear system at state x.

where JacOp evaluates the Jacobian operator at state x.

where KrylovSolve evaluates the solution to the linear system.

where ftol is the error or residual tolerance for convergence.

# Linesearch Methods

• `StandardNewtonStep`

Takes a standard Newton step without any linesearching.

``````% Invoked via the 'options' struct() for Newton methods
options.linesearch = 'none';
``````

• `BacktrackingLinesearch`

Attempts to take a standard Newton step, then if the step fails to show sufficient reduction in the norm, it reduces the step size by the 'contraction_rate' factor. Reduction of step size stops when it reaches the 'min_step' tolerance, regardless of whether or not sufficient norm reduction was achieved.

Guarantees that the residual reduction is monotonic, but may result in step sizes too small to get a solution in reasonable timeframe.

``````% Invoked via the 'options' struct() for Newton methods
options.linesearch = 'backtracking';
options.linesearch_opts = struct('min_step',1e-3,'contraction_rate',0.5);
``````

# Preconditioning

Users may provide any custom preconditioning matrices or function handles to the `JacobianFreeNewtonKrylov` method. Those preconditioning functions must have the following format

``````% Custom preconditioner function used to help establish the search direction
%
%       Consider this function to be a Pseudo-Inverse for the following:
%           s = J(x)\b
%               where s = a vector to solve for, J(x) is the current
%               Jacobian at the current non-linear state, and b is the
%               known right-hand side vector
%
%   @param b = Nx1 vector representing the right-hand size of J*s = b
%   @param Jacfun = Function handle for the Jacobian (J = Jacfun(x))
%   @param x = Nx1 vector for the current non-linear state
%   @param options = An optional user defined struct() for any other information
%               you may need to perform this action.
function s = custom_precon(b,Jacfun,x,options)
...
end
``````

You provide this preconditioner to the `JacobianFreeNewtonKrylov` function through the `options` struct().

``````options.krylov_opts = struct();
options.krylov_opts.M1 = @(b,Jacfun,x,options) custom_precon(b,Jacfun,x,options);
``````

NOTE: We use the same preconditioning syntax from the standard Matlab library of Iterative Krylov Methods (see here). This allows users to provide either a single preconditioner Matrix/Function [M1], or a pair of preconditioners [M1 and M2]. If you provide both preconditioners, note that they are resolve in the following order: M2(M1\b).

# Automatic Matrix Reordering for JFNK

When using Krylov methods to iteratively solve the linear sub-problems, often times the condition number of the matrix will make solutions difficult to accomplish (even when preconditioning!). Thus, the `JacobianFreeNewtonKrylov` implementation here also allows the user to specify a reordering scheme.

``````% Example - Specifying a reordering scheme
options.use_matrix = true;    % Tells the method to use either the user
% provided Jacobian function or a numerical
% Jacobian matrix to apply ordering to

options.equilibrate = true;   % Tells the method that you will request an
% automatic reordering of the Jacobian

options.reordering_method = 'dissect';  % Tells the method to use the nested
% dissection method for reordering

``````

This option is only available to you if you specify to `use_matrix` for the solve. Otherwise, if you do not set `use_matrix` to `true`, the `equilibrate` option will just be ignored.

There are 3 supported reordering methods:

IMPORTANT NOTE: You MUST supply all preconditioners to the `JacobianFreeNewtonKrylov` function as function handles, AND you must use the newly calculated Jacobian from `Jacfun(x)` inside your preconditioning function to formulate the preconditioner result. This is because when doing reordering, the shape/form of the original problem has changed, and may have changed in an unknown way. Thus, using this option means that you may not know what your Jacobian looks like when it gets to your preconditioning function.

# References

** JacobianOperator and JacobianFreeNewtonKrylov**

• Knoll D.A., Keyes D.E. (2003) "Jacobian-Free Newton-Krylov methods: a survey of approaches and applications", Journal of Computational Physics.

** BacktrackLinesearch **

• Armijo, Larry (1966). "Minimization of functions having Lipschitz continuous first partial derivatives". Pacific J. Math. 16 (1): 1–3.

** YBusFormation for Power Flow **

# Citation

Ladshaw, A.P., "PJFNK MATLAB: A MATLAB implementation for Jacobian-Free Newton-Krlov solvers for non-linear systems," https://github.com/aladshaw3/pjfnk_matlab, Accessed (Month) (Day), (Year).

### Cite As

Ladshaw, A.P., "PJFNK MATLAB: A MATLAB implementation for Jacobian-Free Newton-Krlov solvers for non-linear systems," https://github.com/aladshaw3/pjfnk_matlab, Accessed (Month) (Day), (Year).

##### MATLAB Release Compatibility
Created with R2023a
Compatible with any release
##### Platform Compatibility
Windows macOS Linux

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

#### tests

Version Published Release Notes
1.0.0

To view or report issues in this GitHub add-on, visit the GitHub Repository.
To view or report issues in this GitHub add-on, visit the GitHub Repository.