Hello,
As you may have noticed, Scilab 2024.0.0 includes new nonlinear and differential/algebraic equation solvers using the SUNDIALS library (see SUNDIALS | Computing). A common feature of these solvers is the ability to deal with large scale problems by allowing the use of sparse matrices when it is needed. In the following we will show how the new kinsol solver can outperform the legacy fsolve.
In previous Scilab versions, when a nonlinear systems of equations F(y)=0 was solved with fsolve, we could provide the Jacobian F'(y) or let the solver approximate it with finite differences but in both cases the Jacobian had to be given as a dense matrix or approximated as such by the solver. Although this is fine for small problems or when the Jacobian is truly dense, this is a problem for large scale problems like discretized nonlinear partial differential equations.
As an example I will consider the following nonlinear pde:
where u:\Omega\rightarrow\mathbb{R} is constrained on the boundary \Gamma=\partial \Omega by a Dirichlet condition u(x,y)=D(x,y), \forall (x,y)\in\Gamma and u_x, u_y, u_{xy}, u_{xx}, u_{yy} denote the partial derivatives and second derivatives with respect to x and/or y. The solution u of this pde gives the height at each (x,y) of the minimal surface when its boundary is constrained by the given Dirichlet condition (see Minimal surface - Wikipedia).
When \Omega is a rectangle, say [0,1]\times[0,1], the solution of this pde can be approximated by using finite differences quotients to approximate the partial derivatives appearing in the equation and by solving the resulting system of nonlinear equations. In Scilab, we can use the following script:
function out=fun(f)
fx = Dx*f;
fy = Dy*f;
out = (1+fy.*fy).*(Dxx*f)+(1+fx.*fx).*(Dyy*f)-2*fx.*fy.*(Dx*fy);
out(bdy) = f(bdy)-cnd; // boundary condition
endfunction
n = 40;
x = linspace(0,1,n);
[X,Y]=meshgrid(x,x);
// build finite differences operators
h = x(2)-x(1);
d1 = sparse(ones(n-1,1));
d0 = sparse(ones(n,1));
D = (-diag(d1,-1) + diag(d1,1) )/2/h;
// use Kronecker product to build matrices of
// d/dx, d/dy, d^2/dx^2, d^2/dy^2
Dx = D .*. speye(n,n);
Dy = speye(n,n) .*.D;
lap = (diag(d1,-1)+diag(d1,1)-2*diag(d0))/h^2;
Dxx = lap .*. speye(n,n);
Dyy = speye(n,n) .*. lap;
// boundary condition
bdy = find(X(:)==x(1) | X(:)==x($) | Y(:)==x(1) | Y(:)==x($))';
fBnd = 0.5*sign(sin(2*%pi*X).*sin(2*%pi*Y));
cnd = fBnd(bdy);
timer();
[f1,v,info]=fsolve(fBnd(:),fun);
disp("elapsed time = "+string(timer()));
clf
mesh(x,x,matrix(f1,n,n))
We get the following output on a M1 Mac:
elapsed time = 21.665255
and this nice plot:
This time is high because
fsolve approximates the Jacobian of fun as a dense matrix, hence calls fun as many times as the number of columns which is 1600. In fact the Jacobian has the same sparsity pattern as
Dxx+Dyy+Dx*Dy
where at most 9 terms are nonzero on each row, as we can see here (the plot is centered on the (1000,1000) term of the matrix):
If we use
kinsol, giving the pattern of the Jacobian is enough to inform the solver that the Jacobian is sparse and that it has to be approximated by taking account this particular pattern. To this purpose, and as all other SUNDIALS solvers, kinsol uses a particular approximation scheme using the ColPack library (see Sparse Derivative Computation) and in our case the Jacobian is approximated by using only 10 calls to fun().
The call to kinsol to solve the same problem as above is:
timer();
[f2,v,info]=kinsol(fun,fBnd(:), jacPattern=Dxx+Dyy+Dx*Dy, jacUpdateFreq=1);
disp("elapsed time = "+string(timer()));
disp("norm of residual = "+string(norm(v)))
and yields the following output:
kinsol: stopping tolerance on scaled maximum norm of the system function was satisfied.
elapsed time = 0.2178540
On this particular example the kinsol solver is approximately 100 times faster that fsolve. Moreover, since the Jacobian structure does not change with n, the computation time of kinsol scales with the number of unknowns but with a smaller exponent (roughly 2-3 units less) than fsolve (because a dense linear system is solved at each iteration). For example, for n=100, which leads to a system of 10000 (instead of 1600 when n=40) equations, kinsol takes only 1 second to give the solution whereas fsolve takes 8 hours.
So, thanks to SUNDIALS+ColPack, Scilab can now tackle problems which were almost insoluble with the legacy fsolve. But other kinds of problems take advantage of this new infrastructure: I will tell you about the new ode and dae solvers in a forthcoming message.
S.
