using System; namespace Science.Mathematics.VectorCalculus { public class NonlinearFunction { public NonlinearFunction(Function.ToLastType function) { func = function; } private double[] ans; private Function.ToLastType func; public void FindRoot(double[] initialPoint) { ans = new double[initialPoint.Length]; for (int j = 0; j < initialPoint.Length; j++) ans[j] = initialPoint[j]; Newt dd = new Newt(); dd.newt(ans, func); } public double[] Root { get { return ans; } } } class Lnsrch { private double[] x; public double[] X // output { get { return x; } } bool check; public bool Check // output { get { return check; } } double f; public double F // output { get { return f; } } public Lnsrch() { } public void lnsrch(double[] xold, double fold, double[] g, double[] p, double stpmax, Function.ToLastType func) { double ALF = 1.0e-4; double TOLX = 1.0e-15; double a, alam, alam2 = 0.0, alamin, b, disc, f2 = 0.0; double rhs1, rhs2, slope = 0.0, sum = 0.0, temp, test, tmplam; int i, n = xold.Length; x = new double[n]; check = false; for (i = 0; i < n; i++) sum += p[i] * p[i]; sum = Math.Sqrt(sum); if (sum > stpmax) for (i = 0; i < n; i++) p[i] *= stpmax / sum; for (i = 0; i < n; i++) slope += g[i] * p[i]; if (slope >= 0.0) try { throw new Science.Error(); } catch (Science.Error e) { e.Write("Roundoff problem in lnsrch."); } test = 0.0; for (i = 0; i < n; i++) { temp = Math.Abs(p[i]) / Math.Max(Math.Abs(xold[i]), 1.0); if (temp > test) test = temp; } alamin = TOLX / test; alam = 1.0; for (; ; ) { for (i = 0; i < n; i++) x[i] = xold[i] + alam * p[i]; f = func(x); if (alam < alamin) { for (i = 0; i < n; i++) x[i] = xold[i]; check = true; return; } else if (f <= fold + ALF * alam * slope) return; else { if (alam == 1.0) tmplam = -slope / (2.0 * (f - fold - slope)); else { rhs1 = f - fold - alam * slope; rhs2 = f2 - fold - alam2 * slope; a = (rhs1 / (alam * alam) - rhs2 / (alam2 * alam2)) / (alam - alam2); b = (-alam2 * rhs1 / (alam * alam) + alam * rhs2 / (alam2 * alam2)) / (alam - alam2); if (a == 0.0) tmplam = -slope / (2.0 * b); else { disc = b * b - 3.0 * a * slope; if (disc < 0.0) tmplam = 0.5 * alam; else if (b <= 0.0) tmplam = (-b + Math.Sqrt(disc)) / (3.0 * a); else tmplam = -slope / (b + Math.Sqrt(disc)); } if (tmplam > 0.5 * alam) tmplam = 0.5 * alam; } } alam2 = alam; f2 = f; alam = Math.Max(tmplam, 0.1 * alam); } } } class Newt { private bool check; public bool Check // output { get { return check; } } public Newt() { } public void newt(double[] x, Function.ToLastType vecfunc) // x is input and output { int MAXITS = 200; double TOLF = 1.0e-8; double TOLMIN = 1.0e-12; double TOLX = 10e-15; double STPMX = 100.0; int i, its, j, n = x.Length; double den, f, fold, stpmax, sum, temp, test; double[,] fjac = new Double[n, n]; double[] g = new Double[n]; double[] p = new Double[n]; double[] xold = new Double[n]; NRfdjac fd = new NRfdjac(vecfunc); NRfmin fm = new NRfmin(vecfunc); f = fm.Evaluate(x); double[] fvec = fm.Fvec; test = 0.0; for (i = 0; i < n; i++) if (Math.Abs(fvec[i]) > test) test = Math.Abs(fvec[i]); if (test < 0.01 * TOLF) { check = false; return; } sum = 0.0; for (i = 0; i < n; i++) sum += x[i] * x[i]; stpmax = STPMX * Math.Max(Math.Sqrt(sum), (double)n); for (its = 0; its < MAXITS; its++) { fvec = vecfunc(x); fjac = fd.Evaluate(x, fvec); for (i = 0; i < n; i++) { sum = 0.0; for (j = 0; j < n; j++) sum += fjac[j, i] * fvec[j]; g[i] = sum; } for (i = 0; i < n; i++) xold[i] = x[i]; fold = f; for (i = 0; i < n; i++) p[i] = -fvec[i]; LUdcmp obj = new LUdcmp(fjac); obj.solve(p, p); Function.ToLastType func = new Function.ToLastType(fm.Evaluate); Lnsrch ln = new Lnsrch(); ln.lnsrch(xold, fold, g, p, stpmax, func); check = ln.Check; for (i = 0; i < n; i++) x[i] = ln.X[i]; f = ln.F; test = 0.0; for (i = 0; i < n; i++) if (Math.Abs(fvec[i]) > test) test = Math.Abs(fvec[i]); if (test < TOLF) { check = false; return; } if (check) { test = 0.0; den = Math.Max(f, 0.5 * n); for (i = 0; i < n; i++) { temp = Math.Abs(g[i]) * Math.Max(Math.Abs(x[i]), 1.0) / den; if (temp > test) test = temp; } check = (test < TOLMIN ? true : false); return; } test = 0.0; for (i = 0; i < n; i++) { temp = (Math.Abs(x[i] - xold[i])) / Math.Max(Math.Abs(x[i]), 1.0); if (temp > test) test = temp; } if (test < TOLX) return; } try { throw new Science.Error(); } catch (Science.Error e) { e.Write("MAXITS exceeded in newt"); } } } class NRfdjac { private Function.ToLastType func; public NRfdjac(Function.ToLastType funcc) { func = funcc; } public double[,] Evaluate(double[] x, double[] fvec) { double EPS = 1.0e-8; int n = x.Length; double[,] df = new double[n, n]; double[] xh = x; for (int j = 0; j < n; j++) { double temp = xh[j]; double h = EPS * Math.Abs(temp); if (h == 0.0) h = EPS; xh[j] = temp + h; h = xh[j] - temp; double[] f = func(xh); xh[j] = temp; for (int i = 0; i < n; i++) df[i, j] = (f[i] - fvec[i]) / h; } return df; } } class NRfmin { private double[] fvec; public double[] Fvec { get { return fvec; } } private Function.ToLastType jv; public NRfmin(Function.ToLastType fv) { jv = fv; } public double Evaluate(double[] x) { int n = x.Length; double sum = 0.0; fvec = jv(x); for (int i = 0; i < n; i++) sum += fvec[i] * fvec[i]; return 0.5 * sum; } } }