# Using kinsol to solve a nonlinear pde

**URL:** <https://scilab.discourse.group/t/using-kinsol-to-solve-a-nonlinear-pde/304>\
**Category:** Modeling and simulation\
**Created:** [December 20, 2023, 5:18pm UTC](https://scilab.discourse.group/t/using-kinsol-to-solve-a-nonlinear-pde/304 "2023-12-20T17:18:23Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![mottelet](https://yyz2.discourse-cdn.com/free1/user_avatar/scilab.discourse.group/mottelet/32/13_2.png) [@mottelet](https://scilab.discourse.group/u/mottelet)\
**Post date:** [December 20, 2023, 5:18pm UTC](https://scilab.discourse.group/t/using-kinsol-to-solve-a-nonlinear-pde/304/1 "2023-12-20T17:18:23Z")

</div>

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](https://computing.llnl.gov/projects/sundials)). 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:

(1+u\_x^2)u\_{yy}-2u\_xu\_yu\_{xy}+(1+u\_y^2)u\_{xx}=0,\,\mbox{in}\,\Omega\subset\mathbb{R}^2,

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](https://en.wikipedia.org/wiki/Minimal_surface)).

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:

```scilab
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:

```auto
elapsed time = 21.665255

```

and this nice plot:

 ![Screenshot 2023-12-20 at 15.42.21](https://global.discourse-cdn.com/free1/uploads/utc/original/1X/69b87b29a34749776f78462efd62bec211acdaa7.png)  
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

```auto
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):

 ![Screenshot 2023-12-20 at 17.13.19](https://global.discourse-cdn.com/free1/uploads/utc/original/1X/0a590a1f50f1922e722f9040ad74dbda93bb487d.png)  
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](https://cscapes.cs.purdue.edu/coloringpage/)) 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:

```auto
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:

```auto
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.

---

<div class="post-metadata">

**Author:** ![mottelet](https://yyz2.discourse-cdn.com/free1/user_avatar/scilab.discourse.group/mottelet/32/13_2.png) [@mottelet](https://scilab.discourse.group/u/mottelet)\
**Post date:** [December 21, 2023, 11:21am UTC](https://scilab.discourse.group/t/using-kinsol-to-solve-a-nonlinear-pde/304/2 "2023-12-21T11:21:07Z")

</div>

Hello,

As a comparison, I have solved the same problem (for n=160) with the latest native arm64 release 2023b of Matlab, with fsolve and “trust-region” algorithm:

```scilab
n = 160;
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 = kron(D,speye(n,n));
Dy = kron(speye(n,n),D);
lap = (diag(d1,-1)+diag(d1,1)-2*diag(d0))/h^2;
Dxx = kron(lap,speye(n,n));
Dyy = kron(speye(n,n),lap);

% boundary condition
bdy = find(X(:)==x(1) | X(:)==x(end) | Y(:)==x(1) | Y(:)==x(end));
fBnd = 0.5*sign(sin(2*pi*X).*sin(2*pi*Y));
cnd = fBnd(bdy);

pat = Dxx+Dyy+Dx*Dy;
pat(abs(pat)>0)=1;
opt = optimoptions("fsolve","jacobPattern",pat,"Algorithm","trust-region");
tic
[f1,v,info]=fsolve(@(f) fun(f,Dx,Dy,Dxx,Dyy,bdy,cnd),fBnd(:),opt);
toc
disp("norm of residual = "+string(norm(v)))

mesh(x,x,reshape(f1,n,n))

function out=fun(f,Dx,Dy,Dxx,Dyy,bdy,cnd)
    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
end

```

I got the following output and graph (the script has been launched several times in order to benefit of eventual JIT speedup):

```auto
Equation solved, inaccuracy possible.

fsolve stopped because the vector of function values is near zero, as measured by the value
of the function tolerance. However, the last step was ineffective.
<stopping criteria details>
Elapsed time is 18.810963 seconds.

norm of residual = 2.7544e-07

```

 ![Screenshot 2023-12-21 at 11.48.45](https://global.discourse-cdn.com/free1/uploads/utc/original/1X/e26021bfee5fc4fc9450bc4d572d0e15e7e72fe3.png)

The Scilab program (also for n=160) is the following:

```auto
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 = 160;
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,s3]=kinsol(fun,fBnd(:),jacPattern=Dxx+Dyy+Dx*Dy,jacUpdateFreq=1);
disp("elapsed time = "+string(timer()))
disp("norm of residual = "+string(norm(v)))

mesh(x,x,matrix(f1,n,n))

```

and its output is

```auto
  "elapsed time = 7.0518350"
  "norm of residual = 7.451D-09"

```

 ![Screenshot 2023-12-21 at 12.12.40](https://global.discourse-cdn.com/free1/uploads/utc/original/1X/d992aa74215403692b97979491e88a564d6462b0.png)

Hence Scilab is 2.5x faster than Matlab on this particular problem !