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 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