Multidimensional Newton-Raphson method
Multidimensional Newton-Raphson method is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.
- Task
Create a program that finds and outputs the root of a system of nonlinear equations
using Newton-Raphson method.
C#
For matrix inversion and matrix and vector definitions - see C# source from Gaussian elimination
using System;
namespace Rosetta
{
internal interface IFun
{
double F(int index, Vector x);
double df(int index, int derivative, Vector x);
double[] weights();
}
class Newton
{
internal Vector Do(int size, IFun fun, Vector start)
{
Vector X = start.Clone();
Vector F = new Vector(size);
Matrix J = new Matrix(size, size);
Vector D;
do
{
for (int i = 0; i < size; i++)
F[i] = fun.F(i, X);
for (int i = 0; i < size; i++)
for (int j = 0; j < size; j++)
J[i, j] = fun.df(i, j, X);
J.ElimPartial(F);
X -= F;
//need weight vector because different coordinates can diffs by order of magnitudes
} while (F.norm(fun.weights()) > 1e-12);
return X;
}
}
}
using System;
//example from https://eti.pg.edu.pl/documents/176593/26763380/Wykl_AlgorOblicz_7.pdf
namespace Rosetta
{
class Program
{
class Fun: IFun
{
private double[] w = new double[] { 1,1 };
public double F(int index, Vector x)
{
switch (index)
{
case 0: return Math.Atan(x[0]) - x[1] * x[1] * x[1];
case 1: return 4 * x[0] * x[0] + 9 * x[1] * x[1] - 36;
}
throw new Exception("bad index");
}
public double df(int index, int derivative, Vector x)
{
switch (index)
{
case 0:
switch (derivative)
{
case 0: return 1 / (1 + x[0] * x[0]);
case 1: return -3*x[1]*x[1];
}
break;
case 1:
switch (derivative)
{
case 0: return 8 * x[0];
case 1: return 18 * x[1];
}
break;
}
throw new Exception("bad index");
}
public double[] weights() { return w; }
}
static void Main(string[] args)
{
Fun fun = new Fun();
Newton newton = new Newton();
Vector start = new Vector(new double[] { 2.75, 1.25 });
Vector X = newton.Do(2, fun, start);
X.print();
}
}
}
- Output:
2.54258545959024 1.06149981539336
C++
#include <cmath>
#include <cstdint>
#include <iostream>
#include <functional>
#include <stdexcept>
#include <vector>
std::vector<std::vector<double>> to_reduced_row_echelon_form(std::vector<std::vector<double>> vec) {
const uint32_t row_count = vec.size();
const uint32_t col_count = vec[0].size();
uint32_t lead = 0;
for ( uint32_t row = 0; row < row_count; ++row ) {
if ( col_count <= lead ) { return vec; }
uint32_t i = row;
while ( vec[i][lead] == 0 ) {
i += 1;
if ( row_count == i ) {
i = row;
lead += 1;
if ( col_count == lead ) { return vec; }
}
}
const std::vector<double> temp = vec[i]; vec[i] = vec[row]; vec[row] = temp;
if ( vec[row][lead] != 0 ) {
const double divisor = vec[row][lead];
for ( uint32_t col = 0; col < col_count; ++col ) {
vec[row][col] /= divisor;
}
}
for ( uint32_t k = 0; k < row_count; ++k ) {
if ( k != row ) {
const double multiplier = vec[k][lead];
for ( uint32_t j = 0; j < col_count; ++j ) {
vec[k][j] -= vec[row][j] * multiplier;
}
}
}
lead += 1;
}
return vec;
}
std::vector<std::vector<double>> inverse(const std::vector<std::vector<double>>& vec) {
if ( vec.size() != vec[0].size() ) {
throw std::invalid_argument("Not a square vector");
}
const uint32_t size = vec.size();
std::vector<std::vector<double>> augmented(size, std::vector<double>(2 * size, 0.0));
for ( uint32_t row = 0; row < size; ++row ) {
for ( uint32_t col = 0; col < size; ++col ) {
augmented[row][col] = vec[row][col];
}
// Copy identity matrix to the right hand side of augmented matrix
augmented[row][row + size] = 1.0;
}
augmented = to_reduced_row_echelon_form(augmented);
std::vector<std::vector<double>> result(size, std::vector<double>(size, 0.0));
// Copy inverse matrix from right hand side of augmented matrix
for ( uint32_t row = 0; row < size; ++row ) {
for ( uint32_t col = 0; col < size; ++col ) {
result[row][col] = augmented[row][col + size];
}
}
return result;
}
std::vector<std::vector<double>> multiply(const std::vector<std::vector<double>>& one,
const std::vector<std::vector<double>>& two) {
if ( one[0].size() != two.size() ) {
throw std::invalid_argument("Incompatible vector sizes");
}
const uint32_t row_count = one.size();
const uint32_t col_count = two[0].size();
std::vector<std::vector<double>> result(row_count, std::vector<double>(col_count, 0.0));
for ( uint32_t row = 0; row < row_count; ++row ) {
for ( uint32_t col = 0; col < col_count; ++col ) {
for ( uint32_t k = 0; k < row_count; ++k ) {
result[row][col] += one[row][k] * two[k][col];
}
}
}
return result;
}
std::vector<double> solve(
const std::vector<std::function<double(std::vector<double>)>>& functions,
const std::vector<std::vector<std::function<double(std::vector<double>)>>>& jacobian,
std::vector<double> values) {
const uint32_t size = functions.size();
const double epsilon = 0.000'000'08;
const uint32_t maxIterations = 4;
uint32_t iteration = 0;
double maxChange = 0.0;
while ( iteration < maxIterations || maxChange < epsilon ) {
std::vector<std::vector<double>> column(size, std::vector<double>(1, 0.0));
for ( uint32_t i = 0; i < size; ++i ) {
column[i][0] = functions[i](values);
}
std::vector<std::vector<double>> jac(size, std::vector<double>(values.size(), 0.0));
for ( uint32_t i = 0; i < size; ++i ) {
for ( uint32_t j = 0; j < size; ++j ) {
jac[i][j] = jacobian[i][j](values);
}
}
const std::vector<std::vector<double>> jacInverse = inverse(jac);
std::vector<std::vector<double>> delta = multiply(jacInverse, column);
for ( uint32_t i = 0; i < size; ++i ) {
values[i] -= delta[i][0];
if ( std::abs(delta[i][0]) > maxChange ) {
maxChange = std::abs(delta[i][0]);
}
}
iteration += 1;
}
return values;
}
int main() {
/**
* Solve the two non-linear equations:
* y + x^2 - x - 0.5 = 0
* y + 5xy - x^2 = 0
* with initial values of x = 1.2, y = 1.2
*/
std::vector<std::function<double(std::vector<double>)>> functions = {
[](const std::vector<double>& a) { return a[1] + a[0] * a[0] - a[0] - 0.5; },
[](const std::vector<double>& a) { return a[1] + 5.0 * a[0] * a[1] - a[0] * a[0]; }
};
std::vector<std::vector<std::function<double(std::vector<double>)>>> jacobian = {
{ [](const std::vector<double>& a) { return 2.0 * a[0] - 1.0; },
[](const std::vector<double>& a) { return 1.0; }
},
{ [](const std::vector<double>& a) { return 5.0 * a[1] - 2.0 * a[0]; },
[](const std::vector<double>& a) { return 1.0 + 5.0 * a[0]; }
}
};
std::vector<double> initial_values = { 1.2, 1.2 };
std::vector<double> result = solve(functions, jacobian, initial_values);
std::cout << "x = " << result[0] << ", y = " << result[1] << std::endl;
/**
* Solve the three non-linear equations:
* 9x^2 + 36y^2 + 4z^2 - 36 = 0
* x^2 - 2y^2 - 20z = 0
* x^2 - y^2 + z^2 = 0
* with initial values of x = 1.0, y = 1.0 and z = 0.0
*/
functions = {
[](const std::vector<double>& a)
{ return 9.0 * a[0] * a[0] + 36 * a[1] * a[1] + 4 * a[2] * a[2] - 36.0; },
[](const std::vector<double>& a) { return a[0] * a[0] - 2.0 * a[1] * a[1] - 20.0 * a[2]; },
[](const std::vector<double>& a) { return a[0] * a[0] - a[1] * a[1] + a[2] * a[2] ; }
};
jacobian = {
{ [](const std::vector<double>& a) { return 18.0 * a[0]; },
[](const std::vector<double>& a) { return 72.0 * a[1]; },
[](const std::vector<double>& a) { return 8.0 * a[2]; }
},
{ [](const std::vector<double>& a) { return 2.0 * a[0]; },
[](const std::vector<double>& a) { return -4.0 * a[1]; },
[](const std::vector<double>& a) { return -20.0; }
},
{ [](const std::vector<double>& a) { return 2.0 * a[0]; },
[](const std::vector<double>& a) { return -2.0 * a[1]; },
[](const std::vector<double>& a) { return 2.0 * a[2]; }
}
};
initial_values = { 1.0, 1.0, 0.0 };
result = solve(functions, jacobian, initial_values);
std::cout << "x = " << result[0] << ", y = " << result[1] << ", z = " << result[2] <<std::endl;
}
- Output:
x = 1.23332, y = 0.212245 x = 0.893628, y = 0.894527, z = -0.0400893
FreeBASIC
Type Matrix
Dim As Double m(Any, Any)
End Type
Dim Shared As Matrix gu1
Sub rowswap(Byref M As Matrix, i As Uinteger, j As Uinteger)
Dim As Integer k
For k = 0 To Ubound(M.m, 2)
Swap M.m(j, k), M.m(i, k)
Next k
End Sub
Function rowech(Byref M As Matrix) As Matrix
Dim As Uinteger lead = 0, rowCount = 1 + Ubound(M.m, 1), colCount = 1 + Ubound(M.m, 2)
Dim As Uinteger r, i, j
Dim As Double K
For r = 0 To rowCount - 1
If lead >= colCount Then Exit For
i = r
While M.m(i, lead) = 0
i += 1
If i = rowCount Then
i = r
lead += 1
If lead = colCount Then Exit For
End If
Wend
rowswap M, r, i
K = M.m(r, lead)
If K <> 0 Then
For j = 0 To colCount - 1
M.m(r, j) /= K
Next j
End If
For i = 0 To rowCount - 1
If i <> r Then
K = M.m(i, lead)
For j = 0 To colCount - 1
M.m(i, j) -= M.m(r, j) * K
Next j
End If
Next i
lead += 1
Next r
Return M
End Function
Function MatrixTranspose(Byref A As Matrix) As Matrix
Dim mtranspuesta As Matrix
Redim mtranspuesta.m(Lbound(A.m, 2) To Ubound(A.m, 2), Lbound(A.m, 1) To Ubound(A.m, 1))
Dim As Integer fila, columna
For fila = Lbound(A.m, 1) To Ubound(A.m, 1)
For columna = Lbound(A.m, 2) To Ubound(A.m, 2)
mtranspuesta.m(columna, fila) = A.m(fila, columna)
Next columna
Next fila
Return mtranspuesta
End Function
Function MatrixInverse(Byref A As Matrix) As Matrix
Dim ret As Matrix, working As Matrix
Dim As Uinteger i, j
If Ubound(A.m, 1) <> Ubound(A.m, 2) Then Return ret
Dim As Integer n = Ubound(A.m, 1)
Redim ret.m(n, n)
Redim working.m(n, 2 * (n + 1))
For i = 0 To n
For j = 0 To n
working.m(i, j) = A.m(i, j)
Next j
working.m(i, i + n + 1) = 1
Next i
working = rowech(working)
For i = 0 To n
For j = 0 To n
ret.m(i, j) = working.m(i, j + n + 1)
Next j
Next i
Return ret
End Function
Function MatrixMultiply(Byref A As Matrix, Byref B As Matrix) As Matrix
Dim As Matrix result
Redim result.m(Ubound(A.m, 1), Ubound(B.m, 2))
For i As Integer = 0 To Ubound(A.m, 1)
For j As Integer = 0 To Ubound(B.m, 2)
result.m(i, j) = 0
For k As Integer = 0 To Ubound(A.m, 2)
result.m(i, j) += A.m(i, k) * B.m(k, j)
Next k
Next j
Next i
Return result
End Function
Function MatrixSubtract(Byref A As Matrix, Byref B As Matrix) As Matrix
Dim As Matrix result
Redim result.m(Ubound(A.m, 1), Ubound(A.m, 2))
For i As Integer = 0 To Ubound(A.m, 1)
For j As Integer = 0 To Ubound(A.m, 2)
result.m(i, j) = A.m(i, j) - B.m(i, j)
Next j
Next i
Return result
End Function
Type FuncType As Function(x() As Double) As Double
Type JacobianType As Function(x() As Double, i As Integer, j As Integer) As Double
Function solve(funcs() As FuncType, Byref jacobian As JacobianType, guesses() As Double) As Double Ptr
Dim As Integer size = Ubound(funcs)
Dim As Double Ptr gu2 = Callocate((size + 1) * Sizeof(Double))
Dim As Matrix jac
Redim jac.m(size, size)
Dim As Double tol = 1e-8
Dim As Integer maxIter = 12
Dim As Integer iter = 0
Dim As Integer i, j, k
Dim As Double temp_x(size)
For i = 0 To size
gu2[i] = guesses(i)
Next
Do
Dim As Matrix f
Redim f.m(size, 0)
For i = 0 To size
For j = 0 To size
temp_x(j) = gu2[j]
Next j
f.m(i, 0) = funcs(i)(temp_x())
Next
For i = 0 To size
For j = 0 To size
For k = 0 To size
temp_x(k) = gu2[k]
Next k
jac.m(i, j) = jacobian(temp_x(), i, j)
Next j
Next i
Dim As Matrix jacInv = MatrixInverse(jac)
Dim As Matrix delta = MatrixMultiply(jacInv, f)
Dim maxChange As Double = 0
For i = 0 To size
gu2[i] -= delta.m(i, 0)
If Abs(delta.m(i, 0)) > maxChange Then
maxChange = Abs(delta.m(i, 0))
End If
Next
iter += 1
If iter >= maxIter Or maxChange < tol Then Exit Do
Loop
Return gu2
End Function
' Funciones para el primer sistema de ecuaciones
Function f1(x() As Double) As Double
Return -x(0) * x(0) + x(0) + 0.5 - x(1)
End Function
Function f2(x() As Double) As Double
Return x(1) + 5 * x(0) * x(1) - x(0) * x(0)
End Function
Function jacobian1(x() As Double, i As Integer, j As Integer) As Double
Select Case As Const i
Case 0
If j = 0 Then Return -2 * x(0) + 1
If j = 1 Then Return -1
Case 1
If j = 0 Then Return 5 * x(1) - 2 * x(0)
If j = 1 Then Return 1 + 5 * x(0)
End Select
Return 0
End Function
Function f3(x() As Double) As Double
Return 9 * x(0) * x(0) + 36 * x(1) * x(1) + 4 * x(2) * x(2) - 36
End Function
Function f4(x() As Double) As Double
Return x(0) * x(0) - 2 * x(1) * x(1) - 20 * x(2)
End Function
Function f5(x() As Double) As Double
Return x(0) * x(0) - x(1) * x(1) + x(2) * x(2)
End Function
Function jacobian2(x() As Double, i As Integer, j As Integer) As Double
Select Case As Const i
Case 0
If j = 0 Then Return 18 * x(0)
If j = 1 Then Return 72 * x(1)
If j = 2 Then Return 8 * x(2)
Case 1
If j = 0 Then Return 2 * x(0)
If j = 1 Then Return -4 * x(1)
If j = 2 Then Return -20
Case 2
If j = 0 Then Return 2 * x(0)
If j = 1 Then Return -2 * x(1)
If j = 2 Then Return 2 * x(2)
End Select
Return 0
End Function
' Main program
Dim As FuncType funcs(1)
funcs(0) = @f1
funcs(1) = @f2
Dim As Double guesses(1) = {1.2, 1.2}
Dim As Double Ptr sols = solve(funcs(), @jacobian1, guesses())
Print Using "Approximate solutions are x = #.#######, y = #.#######"; sols[0]; sols[1]
Dim As FuncType funcs2(2)
funcs2(0) = @f3
funcs2(1) = @f4
funcs2(2) = @f5
Dim As Double guesses2(2) = {1.0, 1.0, 0.0}
sols = solve(funcs2(), @jacobian2, guesses2())
Print Using !"\nApproximate solutions are x = #.#######, y = #.#######, z = #.#######"; sols[0]; sols[1]; sols[2]
Deallocate(sols)
Sleep
- Output:
Approximate solutions are x = 1.2333178, y = 0.2122450 Approximate solutions are x = 0.8936282, y = 0.8945270, z = -.0400893
Go
We follow the Kotlin example of coding our own matrix methods rather than using a third party library.
package main
import (
"fmt"
"math"
)
type vector = []float64
type matrix []vector
type fun = func(vector) float64
type funs = []fun
type jacobian = []funs
func (m1 matrix) mul(m2 matrix) matrix {
rows1, cols1 := len(m1), len(m1[0])
rows2, cols2 := len(m2), len(m2[0])
if cols1 != rows2 {
panic("Matrices cannot be multiplied.")
}
result := make(matrix, rows1)
for i := 0; i < rows1; i++ {
result[i] = make(vector, cols2)
for j := 0; j < cols2; j++ {
for k := 0; k < rows2; k++ {
result[i][j] += m1[i][k] * m2[k][j]
}
}
}
return result
}
func (m1 matrix) sub(m2 matrix) matrix {
rows, cols := len(m1), len(m1[0])
if rows != len(m2) || cols != len(m2[0]) {
panic("Matrices cannot be subtracted.")
}
result := make(matrix, rows)
for i := 0; i < rows; i++ {
result[i] = make(vector, cols)
for j := 0; j < cols; j++ {
result[i][j] = m1[i][j] - m2[i][j]
}
}
return result
}
func (m matrix) transpose() matrix {
rows, cols := len(m), len(m[0])
trans := make(matrix, cols)
for i := 0; i < cols; i++ {
trans[i] = make(vector, rows)
for j := 0; j < rows; j++ {
trans[i][j] = m[j][i]
}
}
return trans
}
func (m matrix) inverse() matrix {
le := len(m)
for _, v := range m {
if len(v) != le {
panic("Not a square matrix")
}
}
aug := make(matrix, le)
for i := 0; i < le; i++ {
aug[i] = make(vector, 2*le)
copy(aug[i], m[i])
// augment by identity matrix to right
aug[i][i+le] = 1
}
aug.toReducedRowEchelonForm()
inv := make(matrix, le)
// remove identity matrix to left
for i := 0; i < le; i++ {
inv[i] = make(vector, le)
copy(inv[i], aug[i][le:])
}
return inv
}
// note: this mutates the matrix in place
func (m matrix) toReducedRowEchelonForm() {
lead := 0
rowCount, colCount := len(m), len(m[0])
for r := 0; r < rowCount; r++ {
if colCount <= lead {
return
}
i := r
for m[i][lead] == 0 {
i++
if rowCount == i {
i = r
lead++
if colCount == lead {
return
}
}
}
m[i], m[r] = m[r], m[i]
if div := m[r][lead]; div != 0 {
for j := 0; j < colCount; j++ {
m[r][j] /= div
}
}
for k := 0; k < rowCount; k++ {
if k != r {
mult := m[k][lead]
for j := 0; j < colCount; j++ {
m[k][j] -= m[r][j] * mult
}
}
}
lead++
}
}
func solve(fs funs, jacob jacobian, guesses vector) vector {
size := len(fs)
var gu1 vector
gu2 := make(vector, len(guesses))
copy(gu2, guesses)
jac := make(matrix, size)
for i := 0; i < size; i++ {
jac[i] = make(vector, size)
}
tol := 1e-8
maxIter := 12
iter := 0
for {
gu1 = gu2
g := matrix{gu1}.transpose()
t := make(vector, size)
for i := 0; i < size; i++ {
t[i] = fs[i](gu1)
}
f := matrix{t}.transpose()
for i := 0; i < size; i++ {
for j := 0; j < size; j++ {
jac[i][j] = jacob[i][j](gu1)
}
}
g1 := g.sub(jac.inverse().mul(f))
gu2 = make(vector, size)
for i := 0; i < size; i++ {
gu2[i] = g1[i][0]
}
iter++
any := false
for i, v := range gu2 {
if math.Abs(v)-gu1[i] > tol {
any = true
break
}
}
if !any || iter >= maxIter {
break
}
}
return gu2
}
func main() {
/*
solve the two non-linear equations:
y = -x^2 + x + 0.5
y + 5xy = x^2
given initial guesses of x = y = 1.2
Example taken from:
http://www.fixoncloud.com/Home/LoginValidate/OneProblemComplete_Detailed.php?problemid=286
Expected results: x = 1.23332, y = 0.2122
*/
f1 := func(x vector) float64 { return -x[0]*x[0] + x[0] + 0.5 - x[1] }
f2 := func(x vector) float64 { return x[1] + 5*x[0]*x[1] - x[0]*x[0] }
fs := funs{f1, f2}
jacob := jacobian{
funs{
func(x vector) float64 { return -2*x[0] + 1 },
func(x vector) float64 { return -1 },
},
funs{
func(x vector) float64 { return 5*x[1] - 2*x[0] },
func(x vector) float64 { return 1 + 5*x[0] },
},
}
guesses := vector{1.2, 1.2}
sol := solve(fs, jacob, guesses)
fmt.Printf("Approximate solutions are x = %.7f, y = %.7f\n", sol[0], sol[1])
/*
solve the three non-linear equations:
9x^2 + 36y^2 + 4z^2 - 36 = 0
x^2 - 2y^2 - 20z = 0
x^2 - y^2 + z^2 = 0
given initial guesses of x = y = 1.0 and z = 0.0
Example taken from:
http://mathfaculty.fullerton.edu/mathews/n2003/FixPointNewtonMod.html (exercise 5)
Expected results: x = 0.893628, y = 0.894527, z = -0.0400893
*/
fmt.Println()
f3 := func(x vector) float64 { return 9*x[0]*x[0] + 36*x[1]*x[1] + 4*x[2]*x[2] - 36 }
f4 := func(x vector) float64 { return x[0]*x[0] - 2*x[1]*x[1] - 20*x[2] }
f5 := func(x vector) float64 { return x[0]*x[0] - x[1]*x[1] + x[2]*x[2] }
fs = funs{f3, f4, f5}
jacob = jacobian{
funs{
func(x vector) float64 { return 18 * x[0] },
func(x vector) float64 { return 72 * x[1] },
func(x vector) float64 { return 8 * x[2] },
},
funs{
func(x vector) float64 { return 2 * x[0] },
func(x vector) float64 { return -4 * x[1] },
func(x vector) float64 { return -20 },
},
funs{
func(x vector) float64 { return 2 * x[0] },
func(x vector) float64 { return -2 * x[1] },
func(x vector) float64 { return 2 * x[2] },
},
}
guesses = vector{1, 1, 0}
sol = solve(fs, jacob, guesses)
fmt.Printf("Approximate solutions are x = %.7f, y = %.7f, z = %.7f\n", sol[0], sol[1], sol[2])
}
- Output:
Approximate solutions are x = 1.2333178, y = 0.2122450 Approximate solutions are x = 0.8936282, y = 0.8945270, z = -0.0400893
Java
import java.util.Arrays;
import java.util.List;
import java.util.function.Function;
public final class MultidimensionalNewtonRaphsonMethod {
public static void main(String[] args) {
/**
* Solve the two non-linear equations:
* y + x^2 - x - 0.5 = 0
* y + 5xy - x^2 = 0
* with initial values of x = 1.2, y = 1.2
*/
List<Function<Double[], Double>> functions = List.of(
a -> a[1] + a[0] * a[0] - a[0] - 0.5,
a -> a[1] + 5.0 * a[0] * a[1] - a[0] * a[0]
);
List<List<Function<Double[], Double>>> jacobian = List.of(
List.of( a -> 2.0 * a[0] - 1.0, a -> 1.0 ),
List.of( a -> 5.0 * a[1] - 2.0 * a[0], a -> 1.0 + 5.0 * a[0] )
);
Double[] initialValues = new Double[] { 1.2, 1.2 };
Double[] result = solve(functions, jacobian, initialValues);
System.out.println(Arrays.toString(result));
/**
* Solve the three non-linear equations:
* 9x^2 + 36y^2 + 4z^2 - 36 = 0
* x^2 - 2y^2 - 20z = 0
* x^2 - y^2 + z^2 = 0
* with initial values of x = 1.0, y = 1.0 and z = 0.0
*/
functions = List.of(
a -> 9.0 * a[0] * a[0] + 36.0 * a[1] * a[1] + 4 * a[2] * a[2] - 36.0,
a -> a[0] * a[0] - 2.0 * a[1] * a[1] - 20.0 * a[2],
a -> a[0] * a[0] - a[1] * a[1] + a[2] * a[2]
);
jacobian = List.of(
List.of( a -> 18.0 * a[0], a -> 72.0 * a[1], a -> 8.0 * a[2] ),
List.of( a -> 2.0 * a[0], a -> -4.0 * a[1], a -> -20.0 ),
List.of( a -> 2.0 * a[0], a -> -2.0 * a[1], a -> 2.0 * a[2] )
);
initialValues = new Double[] { 1.0, 1.0, 0.0 };
result = solve(functions, jacobian, initialValues);
System.out.println(Arrays.toString(result));
}
private static Double[] solve(List<Function<Double[], Double>> functions,
List<List<Function<Double[], Double>>> jacobian,
Double[] values) {
final int size = functions.size();
final double epsilon = 0.000_000_08;
final int maxIterations = 4;
int iteration = 0;
double maxChange = 0.0;
while ( iteration < maxIterations || maxChange < epsilon ) {
double[][] column = new double[size][1];
for ( int i = 0; i < size; i++ ) {
column[i][0] = functions.get(i).apply(values);
}
double[][] jac = new double[size][values.length];
for ( int i = 0; i < size; i++ ) {
for ( int j = 0; j < size; j++ ) {
jac[i][j] = jacobian.get(i).get(j).apply(values);
}
}
double[][] jacInverse = inverse(jac);
double[][] delta = multiply(jacInverse, column);
for ( int i = 0; i < size; i++ ) {
values[i] -= delta[i][0];
if ( Math.abs(delta[i][0]) > maxChange ) {
maxChange = Math.abs(delta[i][0]);
}
}
iteration += 1;
}
return values;
}
private static double[][] multiply(double[][] one, double[][] two) {
if ( one[0].length != two.length ) {
throw new AssertionError("Incompatible array sizes");
}
final int rowCount = one.length;
final int colCount = two[0].length;
double[][] result = new double[rowCount][colCount];
for ( int row = 0; row < rowCount; row++ ) {
for ( int col = 0; col < colCount; col++ ) {
for ( int k = 0; k < rowCount; k++ ) {
result[row][col] += one[row][k] * two[k][col];
}
}
}
return result;
}
private static double[][] inverse(double[][] array) {
if ( array.length != array[0].length ) {
throw new AssertionError("Not a square array");
}
final int size = array.length;
double[][] augmented = new double[size][2 * size];
for ( int row = 0; row < size; row++ ) {
for ( int col = 0; col < size; col++ ) {
augmented[row][col] = array[row][col];
}
// Copy identity matrix to the right hand side of augmented matrix
augmented[row][row + size] = 1.0;
}
augmented = toReducedRowEchelonForm(augmented);
double[][] result = new double[size][size];
// Copy inverse matrix from right hand side of augmented matrix
for ( int row = 0; row < size; row++ ) {
for ( int col = 0; col < size; col++ ) {
result[row][col] = augmented[row][col + size];
}
}
return result;
}
private static double[][] toReducedRowEchelonForm(double[][] array) {
final int rowCount = array.length;
final int colCount = array[0].length;
int lead = 0;
for ( int row = 0; row < rowCount; row++ ) {
if ( colCount <= lead ) { return array; }
int i = row;
while ( array[i][lead] == 0 ) {
i += 1;
if ( rowCount == i ) {
i = row;
lead += 1;
if ( colCount == lead ) { return array; }
}
}
final double[] temp = array[i]; array[i] = array[row]; array[row] = temp;
if ( array[row][lead] != 0 ) {
final double divisor = array[row][lead];
for ( int col = 0; col < colCount; col++ ) {
array[row][col] /= divisor;
}
}
for ( int k = 0; k < rowCount; k++ ) {
if ( k != row ) {
final double multiplier = array[k][lead];
for ( int j = 0; j < colCount; j++ ) {
array[k][j] -= array[row][j] * multiplier;
}
}
}
lead += 1;
}
return array;
}
}
- Output:
[1.2333177930042694, 0.212245014464034] [0.8936282344764825, 0.8945270103905782, -0.040089286159152804]
Julia
NLsolve is a Julia package for nonlinear systems of equations, with the Newton-Raphson method one of the choices for solvers.
# from the NLSolve documentation: to solve
# (x, y) -> ((x+3)*(y^3-7)+18, sin(y*exp(x)-1))
using NLsolve
function f!(F, x)
F[1] = (x[1]+3)*(x[2]^3-7)+18
F[2] = sin(x[2]*exp(x[1])-1)
end
function j!(J, x)
J[1, 1] = x[2]^3-7
J[1, 2] = 3*x[2]^2*(x[1]+3)
u = exp(x[1])*cos(x[2]*exp(x[1])-1)
J[2, 1] = x[2]*u
J[2, 2] = u
end
println(nlsolve(f!, j!, [ 0.1; 1.2], method = :newton))
- Output:
Results of Nonlinear Solver Algorithm * Algorithm: Newton with line-search * Starting Point: [0.1, 1.2] * Zero: [-3.7818e-16, 1.0] * Inf-norm of residuals: 0.000000 * Iterations: 4 * Convergence: true * |x - x'| < 0.0e+00: false * |f(x)| < 1.0e-08: true * Function Calls (f): 5 * Jacobian Calls (df/dx): 4
Kotlin
A straightforward approach multiplying by the inverse of the Jacobian, rather than dividing by f'(x) as one would do in the single dimensional case, which is quick enough here.
As neither the JDK nor the Kotlin Standard Library have matrix functions built in, most of the functions used have been taken from other tasks.
// Version 1.2.31
import kotlin.math.abs
typealias Vector = DoubleArray
typealias Matrix = Array<Vector>
typealias Func = (Vector) -> Double
typealias Funcs = List<Func>
typealias Jacobian = List<Funcs>
operator fun Matrix.times(other: Matrix): Matrix {
val rows1 = this.size
val cols1 = this[0].size
val rows2 = other.size
val cols2 = other[0].size
require(cols1 == rows2)
val result = Matrix(rows1) { Vector(cols2) }
for (i in 0 until rows1) {
for (j in 0 until cols2) {
for (k in 0 until rows2) {
result[i][j] += this[i][k] * other[k][j]
}
}
}
return result
}
operator fun Matrix.minus(other: Matrix): Matrix {
val rows = this.size
val cols = this[0].size
require(rows == other.size && cols == other[0].size)
val result = Matrix(rows) { Vector(cols) }
for (i in 0 until rows) {
for (j in 0 until cols) {
result[i][j] = this[i][j] - other[i][j]
}
}
return result
}
fun Matrix.transpose(): Matrix {
val rows = this.size
val cols = this[0].size
val trans = Matrix(cols) { Vector(rows) }
for (i in 0 until cols) {
for (j in 0 until rows) trans[i][j] = this[j][i]
}
return trans
}
fun Matrix.inverse(): Matrix {
val len = this.size
require(this.all { it.size == len }) { "Not a square matrix" }
val aug = Array(len) { DoubleArray(2 * len) }
for (i in 0 until len) {
for (j in 0 until len) aug[i][j] = this[i][j]
// augment by identity matrix to right
aug[i][i + len] = 1.0
}
aug.toReducedRowEchelonForm()
val inv = Array(len) { DoubleArray(len) }
// remove identity matrix to left
for (i in 0 until len) {
for (j in len until 2 * len) inv[i][j - len] = aug[i][j]
}
return inv
}
fun Matrix.toReducedRowEchelonForm() {
var lead = 0
val rowCount = this.size
val colCount = this[0].size
for (r in 0 until rowCount) {
if (colCount <= lead) return
var i = r
while (this[i][lead] == 0.0) {
i++
if (rowCount == i) {
i = r
lead++
if (colCount == lead) return
}
}
val temp = this[i]
this[i] = this[r]
this[r] = temp
if (this[r][lead] != 0.0) {
val div = this[r][lead]
for (j in 0 until colCount) this[r][j] /= div
}
for (k in 0 until rowCount) {
if (k != r) {
val mult = this[k][lead]
for (j in 0 until colCount) this[k][j] -= this[r][j] * mult
}
}
lead++
}
}
fun solve(funcs: Funcs, jacobian: Jacobian, guesses: Vector): Vector {
val size = funcs.size
var gu1: Vector
var gu2 = guesses.copyOf()
val jac = Matrix(size) { Vector(size) }
val tol = 1.0e-8
val maxIter = 12
var iter = 0
do {
gu1 = gu2
val g = arrayOf(gu1).transpose()
val f = arrayOf(Vector(size) { funcs[it](gu1) }).transpose()
for (i in 0 until size) {
for (j in 0 until size) {
jac[i][j] = jacobian[i][j](gu1)
}
}
val g1 = g - jac.inverse() * f
gu2 = Vector(size) { g1[it][0] }
iter++
}
while (gu2.withIndex().any { iv -> abs(iv.value - gu1[iv.index]) > tol } && iter < maxIter)
return gu2
}
fun main(args: Array<String>) {
/* solve the two non-linear equations:
y = -x^2 + x + 0.5
y + 5xy = x^2
given initial guesses of x = y = 1.2
Example taken from:
http://www.fixoncloud.com/Home/LoginValidate/OneProblemComplete_Detailed.php?problemid=286
Expected results: x = 1.23332, y = 0.2122
*/
val f1: Func = { x -> -x[0] * x[0] + x[0] + 0.5 - x[1] }
val f2: Func = { x -> x[1] + 5 * x[0] * x[1] - x[0] * x[0] }
val funcs = listOf(f1, f2)
val jacobian = listOf(
listOf<Func>({ x -> - 2.0 * x[0] + 1.0 }, { _ -> -1.0 }),
listOf<Func>({ x -> 5.0 * x[1] - 2.0 * x[0] }, { x -> 1.0 + 5.0 * x[0] })
)
val guesses = doubleArrayOf(1.2, 1.2)
val (xx, yy) = solve(funcs, jacobian, guesses)
System.out.printf("Approximate solutions are x = %.7f, y = %.7f\n", xx, yy)
/* solve the three non-linear equations:
9x^2 + 36y^2 + 4z^2 - 36 = 0
x^2 - 2y^2 - 20z = 0
x^2 - y^2 + z^2 = 0
given initial guesses of x = y = 1.0 and z = 0.0
Example taken from:
http://mathfaculty.fullerton.edu/mathews/n2003/FixPointNewtonMod.html (exercise 5)
Expected results: x = 0.893628, y = 0.894527, z = -0.0400893
*/
println()
val f3: Func = { x -> 9.0 * x[0] * x[0] + 36.0 * x[1] * x[1] + 4.0 * x[2] * x[2] - 36.0 }
val f4: Func = { x -> x[0] * x[0] - 2.0 * x[1] * x[1] - 20.0 * x[2] }
val f5: Func = { x -> x[0] * x[0] - x[1] * x[1] + x[2] * x[2] }
val funcs2 = listOf(f3, f4, f5)
val jacobian2 = listOf(
listOf<Func>({ x -> 18.0 * x[0] }, { x -> 72.0 * x[1] }, { x -> 8.0 * x[2] }),
listOf<Func>({ x -> 2.0 * x[0] }, { x -> -4.0 * x[1] }, { _ -> -20.0 }),
listOf<Func>({ x -> 2.0 * x[0] }, { x -> -2.0 * x[1] }, { x -> 2.0 * x[2] })
)
val guesses2 = doubleArrayOf(1.0, 1.0, 0.0)
val (xx2, yy2, zz2) = solve(funcs2, jacobian2, guesses2)
System.out.printf("Approximate solutions are x = %.7f, y = %.7f, z = %.7f\n", xx2, yy2, zz2)
}
- Output:
Approximate solutions are x = 1.2333178, y = 0.2122450 Approximate solutions are x = 0.8936282, y = 0.8945270, z = -0.0400893
Nim
import sequtils, strformat, sugar
type
Vector = seq[float]
Matrix = seq[Vector]
Func = (Vector) -> float
Funcs = seq[Func]
Jacobian = seq[Funcs]
func `*`(m1, m2: Matrix): Matrix =
let
rows1 = m1.len
cols1 = m1[0].len
rows2 = m2.len
cols2 = m2[0].len
doAssert cols1 == rows2
result = newSeqWith(rows1, newSeq[float](cols2))
for i in 0..<rows1:
for j in 0..<cols2:
for k in 0..<rows2:
result[i][j] += m1[i][k] * m2[k][j]
func `-`(m1, m2: Matrix): Matrix =
let
rows = m1.len
cols = m1[0].len
doAssert m2.len == rows and m2[0].len == cols
result = newSeqWith(rows, newSeq[float](cols))
for i in 0..<rows:
for j in 0..<cols:
result[i][j] = m1[i][j] - m2[i][j]
func transposed(m: Matrix): Matrix =
let
rows = m.len
cols = m[0].len
result = newSeqWith(cols, newSeq[float](rows))
for i in 0..<cols:
for j in 0..<rows:
result[i][j] = m[j][i]
func toReducedRowEchelonForm(m: var Matrix) =
var lead = 0
let rowCount = m.len
let colCount = m[0].len
for r in 0..<rowCount:
if colCount <= lead: return
var i = r
while m[i][lead] == 0:
inc i
if rowCount == i:
i = r
inc lead
if colCount == lead: return
swap m[i], m[r]
if m[r][lead] != 0:
let divisor = m[r][lead]
for j in 0..<colCount: m[r][j] /= divisor
for k in 0..<rowCount:
if k != r:
let mult = m[k][lead]
for j in 0..<colCount: m[k][j] -= m[r][j] * mult
inc lead
func inverse(m: Matrix): Matrix =
let size = m.len
doAssert m.allIt(it.len == size), "not a square matrix."
var aug = newSeqWith(size, newSeq[float](2 * size))
for i in 0..<size:
for j in 0..<size: aug[i][j] = m[i][j]
# Augment by identity matrix to right.
aug[i][i + size] = 1
aug.toReducedRowEchelonForm()
result = newSeqWith(size, newSeq[float](size))
# Remove identity matrix to left.
for i in 0..<size:
for j in 0..<size: result[i][j] = aug[i][j + size]
proc solve(funcs: Funcs; jacobian: Jacobian; guesses: Vector): Vector =
let size = funcs.len
result = guesses
var jac = newSeqWith(size, newSeq[float](size))
const Tol = 1e-8
let MaxIter = 12
var iter = 1
while true:
let gu = move(result)
let g = transposed(@[gu])
let f = transposed(@[funcs.mapIt(it(gu))])
for i in 0..<size:
for j in 0..<size:
jac[i][j] = jacobian[i][j](gu)
let g1 = g - inverse(jac) * f
result = g1.mapIt(it[0])
inc iter
if iter > MaxIter: break
var exit = true
for idx, val in result:
if abs(val - gu[idx]) > Tol:
exit = false
break
if exit: break
when isMainModule:
#[ Solve the two non-linear equations:
y = -x^2 + x + 0.5
y + 5xy = x^2
given initial guesses of x = y = 1.2
Example taken from:
http://www.fixoncloud.com/Home/LoginValidate/OneProblemComplete_Detailed.php?problemid=286
Expected results: x = 1.23332, y = 0.2122
]#
let
f1: Func = (x: Vector) => -x[0] * x[0] + x[0] + 0.5 - x[1]
f2: Func = (x: Vector) => x[1] + 5 * x[0] * x[1] - x[0] * x[0]
funcs1: Funcs = @[f1, f2]
jacobian1: Jacobian = @[@[Func((x: Vector) => - 2 * x[0] + 1),
Func((x: Vector) => -1.0)],
@[Func((x: Vector) => 5 * x[1] - 2 * x[0]),
Func((x: Vector) => 1 + 5 * x[0])]]
guesses1 = @[1.2, 1.2]
sol1 = solve(funcs1, jacobian1, guesses1)
echo &"Approximate solutions are x = {sol1[0]:.7f}, y = {sol1[1]:.7f}"
#[ Solve the three non-linear equations:
9x^2 + 36y^2 + 4z^2 - 36 = 0
x^2 - 2y^2 - 20z = 0
x^2 - y^2 + z^2 = 0
given initial guesses of x = y = 1.0 and z = 0.0
Example taken from:
http://mathfaculty.fullerton.edu/mathews/n2003/FixPointNewtonMod.html (exercise 5)
Expected results: x = 0.893628, y = 0.894527, z = -0.0400893
]#
echo()
let
f3: Func = (x: Vector) => 9 * x[0] * x[0] + 36 * x[1] * x[1] + 4 * x[2] * x[2] - 36
f4: Func = (x: Vector) => x[0] * x[0] - 2 * x[1] * x[1] - 20 * x[2]
f5: Func = (x: Vector) => x[0] * x[0] - x[1] * x[1] + x[2] * x[2]
funcs2: Funcs = @[f3, f4, f5]
jacobian2: Jacobian = @[@[Func((x: Vector) => 18 * x[0]),
Func((x: Vector) => 72 * x[1]),
Func((x: Vector) => 8 * x[2])],
@[Func((x: Vector) => 2 * x[0]),
Func((x: Vector) => -4 * x[1]),
Func((x: Vector) => -20.0)],
@[Func((x: Vector) => 2 * x[0]),
Func((x: Vector) => -2 * x[1]),
Func((x: Vector) => 2 * x[2])]]
guesses2 = @[1.0, 1.0, 0.0]
sol2 = solve(funcs2, jacobian2, guesses2)
echo &"Approximate solutions are x = {sol2[0]:.7f}, y = {sol2[1]:.7f}, z = {sol2[2]:.7f}"
- Output:
Approximate solutions are x = 1.2333178, y = 0.2122450 Approximate solutions are x = 0.8936282, y = 0.8945270, z = -0.0400893
Phix
Uses code from Reduced_row_echelon_form#Phix,
Gauss-Jordan_matrix_inversion#Phix,
Matrix_transposition#Phix, and
Matrix_multiplication#Phix
See std distro for a complete runnable version.
-- demo\rosetta\Multidimensional_Newton-Raphson_method.exw with javascript_semantics function solve(sequence fs, jacob, guesses) integer size = length(fs), maxIter = 12, iter = 0 sequence gu1, g, t, f, g1, gu2 = guesses, jac = repeat(repeat(0,size),size) atom tol = 1e-8 while true do gu1 := gu2 g := matrix_transpose({gu1}) t := repeat(0, size) for i=1 to size do t[i] := call_func(fs[i],{gu1}) end for f := matrix_transpose({t}) for i=1 to size do for j=1 to size do jac[i][j] := call_func(jacob[i][j],{gu1}) end for end for g1 := sq_sub(g,matrix_mul(inverse(jac),f)) gu2 := vslice(g1,1) iter += 1 bool any := find(true,sq_gt(sq_sub(sq_abs(gu2),gu1),tol))!=0 if not any or iter >= maxIter then exit end if end while return gu2 end function function f1(sequence v) atom {x,y} = v return -x*x+x+0.5-y end function function f2(sequence v) atom {x,y} = v return y+5*x*y-x*x end function function f3(sequence v) atom {x,y,z} = v return 9*x*x+36*y*y+4*z*z-36 end function function f4(sequence v) atom {x,y,z} = v return x*x-2*y*y-20*z end function function f5(sequence v) atom {x,y,z} = v return x*x-y*y+z*z end function function j1(sequence v) atom {x} = v return -2*x+1 end function function j2(sequence /*v*/) return -1 end function function j3(sequence v) atom {x,y} = v return 5*y-2*x end function function j4(sequence v) atom {x} = v return 1+5*x end function function j11(sequence v) atom {x} = v return 18*x end function function j12(sequence v) atom {?,y} = v return 72*y end function function j13(sequence v) atom {?,?,z} = v return 8*z end function function j21(sequence v) atom {x} = v return 2*x end function function j22(sequence v) atom {?,y} = v return -4*y end function function j23(sequence /*v*/) return -20 end function function j31(sequence v) atom {x} = v return 2*x end function function j32(sequence v) atom {?,y} = v return -2*y end function function j33(sequence v) atom {?,?,z} = v return 2*z end function procedure main() sequence fs, jacob, guesses /* solve the two non-linear equations: y = -x^2 + x + 0.5 y + 5xy = x^2 given initial guesses of x = y = 1.2 Example taken from: http://www.fixoncloud.com/Home/LoginValidate/OneProblemComplete_Detailed.php?problemid=286 Expected results: x = 1.23332, y = 0.2122 */ fs = {f1,f2} jacob = {{j1,j2}, {j3,j4}} guesses := {1.2, 1.2} printf(1,"Approximate solutions are x = %.7f, y = %.7f\n\n", solve(fs, jacob, guesses)) /* solve the three non-linear equations: 9x^2 + 36y^2 + 4z^2 - 36 = 0 x^2 - 2y^2 - 20z = 0 x^2 - y^2 + z^2 = 0 given initial guesses of x = y = 1.0 and z = 0.0 Example taken from: http://mathfaculty.fullerton.edu/mathews/n2003/FixPointNewtonMod.html (exercise 5) Expected results: x = 0.893628, y = 0.894527, z = -0.0400893 */ fs = {f3, f4, f5} jacob = {{j11,j12,j13}, {j21,j22,j23}, {j31,j32,j33}} guesses = {1, 1, 0} printf(1,"Approximate solutions are x = %.7f, y = %.7f, z = %.7f\n", solve(fs, jacob, guesses)) end procedure main()
- Output:
Approximate solutions are x = 1.2333178, y = 0.2122450 Approximate solutions are x = 0.8936282, y = 0.8945270, z = -0.04008929
Raku
(formerly Perl 6)
# Reference:
# https://github.com/pierre-vigier/Perl6-Math-Matrix
# Mastering Algorithms with Perl
# By Jarkko Hietaniemi, John Macdonald, Jon Orwant
# Publisher: O'Reilly Media, ISBN-10: 1565923987
# https://resources.oreilly.com/examples/9781565923980/blob/master/ch16/solve
use v6;
sub solve_funcs ($funcs, @guesses, $iterations, $epsilon) {
my ($error_value, @values, @delta, @jacobian); my \ε = $epsilon;
for 1 .. $iterations {
for ^+$funcs { @values[$^i] = $funcs[$^i](|@guesses); }
$error_value = 0;
for ^+$funcs { $error_value += @values[$^i].abs }
return @guesses if $error_value ≤ ε;
for ^+$funcs { @delta[$^i] = -@values[$^i] }
@jacobian = jacobian $funcs, @guesses, ε;
@delta = solve_matrix @jacobian, @delta;
loop (my $j = 0, $error_value = 0; $j < +$funcs; $j++) {
$error_value += @delta[$j].abs ;
@guesses[$j] += @delta[$j];
}
return @guesses if $error_value ≤ ε;
}
return @guesses;
}
sub jacobian ($funcs is copy, @points is copy, $epsilon is copy) {
my ($Δ, @P, @M, @plusΔ, @minusΔ);
my Array @jacobian; my \ε = $epsilon;
for ^+@points -> $i {
@plusΔ = @minusΔ = @points;
$Δ = (ε * @points[$i].abs) || ε;
@plusΔ[$i] = @points[$i] + $Δ;
@minusΔ[$i] = @points[$i] - $Δ;
for ^+$funcs { @P[$^k] = $funcs[$^k](|@plusΔ); }
for ^+$funcs { @M[$^k] = $funcs[$^k](|@minusΔ); }
for ^+$funcs -> $j {
@jacobian[$j][$i] = (@P[$j] - @M[$j]) / (2 * $Δ);
}
}
return @jacobian;
}
sub solve_matrix (@matrix_array is copy, @delta is copy) {
# https://github.com/pierre-vigier/Perl6-Math-Matrix/issues/56
{ use Math::Matrix;
my $matrix = Math::Matrix.new(@matrix_array);
my $vector = Math::Matrix.new(@delta.map({.list}));
die "Matrix is not invertible" unless $matrix.is-invertible;
my @result = ( $matrix.inverted dot $vector ).transposed;
return @result.split(" ");
}
}
my $funcs = [
{ 9*$^x² + 36*$^y² + 4*$^z² - 36 }
{ $^x² - 2*$^y² - 20*$^z }
{ $^x² - $^y² + $^z² }
];
my @guesses = (1,1,0);
my @solution = solve_funcs $funcs, @guesses, 20, 1e-8;
say "Solution: ", @solution;
- Output:
Solution: [0.8936282344764825 0.8945270103905782 -0.04008928615915281]
Wren
import "./matrix" for Matrix
import "./fmt" for Fmt
var solve = Fn.new { |funcs, jacobian, guesses|
var size = funcs.count
var gu1 = []
var gu2 = guesses.toList
var jac = Matrix.new(size, size)
var tol = 1e-8
var maxIter = 12
var iter = 0
while (true) {
gu1 = gu2
var g = Matrix.new([gu1]).transpose
var v = List.filled(size, 0)
for (i in 0...size) v[i] = funcs[i].call(gu1)
var f = Matrix.new([v]).transpose
for (i in 0...size) {
for (j in 0...size) jac[i, j] = jacobian[i][j].call(gu1)
}
var g1 = g - jac.inverse * f
gu2 = List.filled(size, 0)
for (i in 0...size) gu2[i] = g1[i][0]
iter = iter + 1
if (iter == maxIter) break
if ((0...size).all { |i| (gu2[i] - gu1[i]).abs <= tol }) break
}
return gu2
}
/* solve the two non-linear equations:
y = -x^2 + x + 0.5
y + 5xy = x^2
given initial guesses of x = y = 1.2
Example taken from:
http://www.fixoncloud.com/Home/LoginValidate/OneProblemComplete_Detailed.php?problemid=286
Expected results: x = 1.23332, y = 0.2122
*/
var f1 = Fn.new { |x| -x[0] * x[0] + x[0] + 0.5 - x[1] }
var f2 = Fn.new { |x| x[1] + 5 * x[0] * x[1] - x[0] * x[0] }
var funcs = [f1, f2]
var jacobian = [
[ Fn.new { |x| - 2 * x[0] + 1 }, Fn.new { |x| -1 } ],
[ Fn.new { |x| 5 * x[1] - 2 * x[0] }, Fn.new { |x| 1 + 5 * x[0] } ]
]
var guesses = [1.2, 1.2]
var sols = solve.call(funcs, jacobian, guesses)
Fmt.print("Approximate solutions are x = $.7f, y = $.7f", sols[0], sols[1])
/* solve the three non-linear equations:
9x^2 + 36y^2 + 4z^2 - 36 = 0
x^2 - 2y^2 - 20z = 0
x^2 - y^2 + z^2 = 0
given initial guesses of x = y = 1.0 and z = 0.0
Example taken from:
http://mathfaculty.fullerton.edu/mathews/n2003/FixPointNewtonMod.html (exercise 5)
Expected results: x = 0.893628, y = 0.894527, z = -0.0400893
*/
System.print()
var f3 = Fn.new { |x| 9 * x[0] * x[0] + 36 * x[1] *x[1] + 4 * x[2] * x[2] - 36 }
var f4 = Fn.new { |x| x[0] * x[0] - 2 * x[1] * x[1] - 20 * x[2] }
var f5 = Fn.new { |x| x[0] * x[0] - x[1] * x[1] + x[2] * x[2] }
funcs = [f3, f4, f5]
jacobian = [
[ Fn.new { |x| 18 * x[0] }, Fn.new { |x| 72 * x[1] }, Fn.new { |x| 8 * x[2] } ],
[ Fn.new { |x| 2 * x[0] }, Fn.new { |x| -4 * x[1] }, Fn.new { |x| -20 } ],
[ Fn.new { |x| 2 * x[0] }, Fn.new { |x| -2 * x[1] }, Fn.new { |x| 2 * x[2] } ]
]
guesses = [1, 1, 0]
sols = solve.call(funcs, jacobian, guesses)
Fmt.print("Approximate solutions are x = $.7f, y = $.7f, z = $.7f", sols[0], sols[1], sols[2])
- Output:
Approximate solutions are x = 1.2333178, y = 0.2122450 Approximate solutions are x = 0.8936282, y = 0.8945270, z = -0.0400893
zkl
This doesn't use Newton-Raphson (with derivatives) but a hybrid algorithm.
var [const] GSL=Import.lib("zklGSL"); // libGSL (GNU Scientific Library)
// two functions of two variables: f(x,y)=0
fs:=T(fcn(x,y){ x.atan() - y*y*y }, fcn(x,y){ 4.0*x*x + 9*y*y - 36 });
v=GSL.VectorFromData(2.75, 1.25); // an initial guess at the solution
GSL.multiroot_fsolver(fs,v);
v.format(11,8).println(); // answer overwrites initial guess
fs.run(True,v.toList().xplode()).println(); // deltas from zero- Output:
2.59807621, 1.06365371 L(2.13651e-09,2.94321e-10)
A condensed solver (for a different set of functions):
v:=GSL.VectorFromData(-10.0, -15.0);
GSL.multiroot_fsolver(T( fcn(x,y){ 1.0 - x }, fcn(x,y){ 10.0*(y - x*x) }),v)
.format().println(); // --> (1,1)- Output:
1.00,1.00
Another example:
v:=GSL.VectorFromData(1.0, 1.0, 0.0); // initial guess
fxyzs:=T(
fcn(x,y,z){ x*x*9 + y*y*36 + z*z*4 - 36 }, // 9x^2 + 36y^2 + 4z^2 - 36 = 0
fcn(x,y,z){ x*x - y*y*2 - z*20 }, // x^2 - 2y^2 - 20z = 0
fcn(x,y,z){ x*x - y*y + z*z }); // x^2 - y^2 + z^2 = 0
(v=GSL.multiroot_fsolver(fxyzs,v)).format(12,8).println();
fxyzs.run(True,v.toList().xplode()).println(); // deltas from zero- Output:
0.89362824, 0.89452701, -0.04008929 L(6.00672e-08,1.0472e-08,9.84017e-09)