First commit

.gitignore

@@ -13,3 +13,5 @@
 
 # Dependency directories (remove the comment below to include it)
 # vendor/
+
+.vscode/*

go.mod

@@ -0,0 +1,3 @@
+module github.com/davidkleiven/nonlin
+
+go 1.13

nonlin/jfsj.go

@@ -0,0 +1,141 @@
+package nonlin
+
+import (
+	"math"
+)
+
+// Problem to be solved
+type Problem struct {
+	// F is a function where we seek the solution of F(x) = 0. The value of the function
+	// should be given to out. out will always be the same length as x. Note that F must
+	// not modify x.
+	F func(out, x []float64)
+}
+
+// Result is a result type returned by the solvers.
+type Result struct {
+	// X is the solution vector at termination
+	X []float64
+
+	// F is the value the function at termination
+	F []float64
+
+	// MaxF is the infinity norm of F
+	MaxF float64
+
+	// Converges is True if the convergence criteria are met, otherwise it is False
+	Converged bool
+}
+
+// Settings is a type that holds parameters required by jfsj algorithm
+type Settings struct {
+	Tol     float64
+	Maxiter int
+}
+
+func defaultSettings() *Settings {
+	return &Settings{
+		Tol:     1e-10,
+		Maxiter: 100000,
+	}
+}
+
+// identify returns the diagonal of the identity matrix
+func identity(n int) []float64 {
+	I := make([]float64, n)
+	for i := range I {
+		I[i] = 1.0
+	}
+	return I
+}
+
+func updateG(G []float64, dF []float64, dx []float64) {
+	dFDotDx := 0.0
+	dFGdF := 0.0
+	dF4 := 0.0
+	for i := range dx {
+		dFDotDx += dF[i] * dx[i]
+		dFGdF += dF[i] * G[i] * dF[i]
+		dF4 += dF[i] * dF[i] * dF[i] * dF[i]
+	}
+
+	for i := range G {
+		G[i] += (dFDotDx - dFGdF) * dF[i] * dF[i] / dF4
+	}
+}
+
+// JFSJ applies the jacobian free singular jacobian algorithm to solve
+// the non-linear system of equations
+func JFSJ(p Problem, x []float64, settings *Settings) Result {
+	G := identity(len(x))
+	if settings == nil {
+		settings = defaultSettings()
+	}
+
+	work := make([]float64, 3*len(x))
+	p.F(work[:len(x)], x)
+
+	var result Result
+	for i := 0; i < settings.Maxiter; i++ {
+		// On even iterations F is placed in work[:len(x)], on odd iterations
+		// F is placed in work[len(x):2*len(x)]
+		F := work[:len(x)]
+		if i%2 == 1 {
+			F = work[len(x) : 2*len(x)]
+		}
+
+		infNormDx := 0.0
+		for j := range x {
+			dx := -G[j] * F[j]
+			x[j] += dx
+			work[2*len(x)+j] = dx
+			if math.Abs(dx) > infNormDx {
+				infNormDx = math.Abs(dx)
+			}
+		}
+
+		dFArray := work[:len(x)]
+		if i%2 == 0 {
+			F = work[len(x) : 2*len(x)]
+			p.F(F, x)
+		} else {
+			F = work[:len(x)]
+			dFArray = work[len(x) : 2*len(x)]
+			p.F(F, x)
+		}
+
+		infNormDF := 0.0
+		infNormF := 0.0
+		for j := range x {
+			dF := math.Abs(work[j] - work[j+len(x)])
+			if dF > infNormDF {
+				infNormDF = dF
+			}
+
+			if math.Abs(F[j]) > infNormF {
+				infNormF = math.Abs(F[j])
+			}
+		}
+
+		if infNormDx+infNormF < settings.Tol {
+			result.X = x
+			result.Converged = true
+			result.F = F
+			return result
+		}
+
+		if infNormDF > settings.Tol {
+			sign := math.Pow(-1.0, float64(i))
+			// Overwrite the part of the work array where F_k is stored with dF_k = F_{k+1} - F_k
+			for j := range x {
+				dFArray[j] = sign * (work[len(x)+j] - work[j])
+			}
+			updateG(G, dFArray, work[2*len(x):])
+		}
+	}
+
+	result.X = x
+	result.Converged = false
+	result.F = work[:len(x)]
+	return result
+}

nonlin/jfsj_test.go

@@ -0,0 +1,78 @@
+package nonlin
+
+import (
+	"math"
+	"testing"
+)
+
+func TestJFSJ(t *testing.T) {
+	for i, test := range []struct {
+		F        func(out []float64, x []float64)
+		Init     []float64
+		Solution []float64
+		Tol      float64
+	}{
+		{
+			F: func(out []float64, x []float64) {
+				out[0] = math.Pow(x[0]-1.0, 2.0) * (x[0] - x[1])
+				out[1] = math.Pow(x[1]-2.0, 5.0) * math.Cos(2.0*x[0]/x[1])
+			},
+			Init:     []float64{0.0, 3.0},
+			Solution: []float64{1.0, 2.0},
+			Tol:      1e-6,
+		},
+		{
+			F: func(out []float64, x []float64) {
+				out[0] = math.Pow(x[0]-1.0, 2.0) * (x[0] - x[1])
+				out[1] = math.Pow(x[1]-2.0, 5.0) * math.Cos(2.0*x[0]/x[1])
+			},
+			Init:     []float64{0.5, 2.0},
+			Solution: []float64{1.0, 2.0},
+			Tol:      1e-6,
+		},
+		{
+			F: func(out []float64, x []float64) {
+				out[0] = 4.0*x[0] - 2.0*x[1] + x[0]*x[0] - 3.0
+				out[1] = -x[0] + 4*x[1] - x[2] + x[1]*x[1] - 3.0
+				out[2] = -2.0*x[1] + 4.0*x[2] + x[2]*x[2] - 3.0
+			},
+			Init:     []float64{-1.5, 0.0, -1.5},
+			Solution: []float64{1.0, 1.0, 1.0},
+			Tol:      1e-6,
+		},
+		{
+			F: func(out []float64, x []float64) {
+				out[0] = 2.0/(1.0+x[0]*x[0]) + math.Sin(x[1]-1.0) - 1.0
+				out[1] = math.Sin(x[0]-1.0) + 2.0/(1.0+x[1]*x[1]) - 1.0
+			},
+			Init:     []float64{0.7, 0.7},
+			Solution: []float64{1.0, 1.0},
+			Tol:      1e-6,
+		},
+		{
+			F: func(out []float64, x []float64) {
+				out[0] = math.Exp(x[0]) + x[1] - 1.0
+				out[1] = math.Exp(x[1]) + x[0] - 1.0
+			},
+			Init:     []float64{-0.5, -0.5},
+			Solution: []float64{0.0, 0.0},
+			Tol:      1e-6,
+		},
+	} {
+		p := Problem{
+			F: test.F,
+		}
+
+		res := JFSJ(p, test.Init, nil)
+		if !res.Converged {
+			t.Errorf("Test #%d: Did not converge\n", i)
+		}
+
+		for j := range res.X {
+			if math.Abs(res.X[j]-test.Solution[j]) > test.Tol {
+				t.Errorf("Test #%d:\nExpected\n%v\nGot\n%v\n", i, test.Solution, res.X)
+				break
+			}
+		}
+	}
+}