Multidimensional Newton-Raphson method
- Task
Create a program that finds and outputs the root of a system of nonlinear equations
using Newton-Raphson metod.
C#
For matrix inversion and matrix and vector definitions - see C# source from Gaussian elimination <lang csharp> 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; } }
} </lang> <lang csharp> 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(); } }
} </lang>
- Output:
2.54258545959024 1.06149981539336
zkl
This doesn't use Newton-Raphson (with derivatives) but a hybrid algorithm. <lang zkl>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
xy:=v.toList(); // Vector to List fs.apply('wrap(f){ f(xy.xplode()) }).println(); // deltas from zero</lang>
- Output:
2.59807621, 1.06365371 L(2.13651e-09,2.94321e-10)
A condensed solver (for a different set of functions): <lang zkl>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)</lang>
- Output:
1.00,1.00