Multidimensional Newton-Raphson method: Difference between revisions

(code - begin)
Line 8:
For matrix inversion and matrix and vector definitions - see C# source from [[Gauss-Jordan matrix inversion]]
<lang C#>
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);
D = J.Inv() * F;
X -= D;
//need weight vector because different coordinates can diffs by order of magnitudes
} while (D.norm(fun.weights())>1e-12);
return X;
}
}
}
</lang>
<lang C#>
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();
}
}
Line 23 ⟶ 101:
</lang>
{{out}}<pre>
2.54258545959024
1.06149981539336
</pre>
Anonymous user