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