$x_1 = x_0 - \displaystyle \frac {2 f(x_0) f'(x_0)} {2 {[f'(x_0)]}^2 - f(x_0) f''(x_0)}$
Uzasadnienie jest dla mnie dość pokrętne. Weźmy funkcję:
$g(x) = \displaystyle\frac {f(x)} {\sqrt{|f'(x)|}}$
Miejsca zerowe g(x) to miejsca zerowe f(x). Miejsca zerowe f(x) o ile nie są miejscami zerowymi $f'(x)$ są miejscami zerowymi g(x). Funkcja g(x) ma wartości nieoznaczone w ekstremach. Funkcja ta została zbudowana w taki sposób nieprzypadkowo. Jeśli podstawimy ją do wzoru na metodę Newtona to otrzymamy:
$x_1=x_0-\displaystyle\frac{g(x_0)}{g'(x_0)}=\frac{f(x_0)}{\sqrt{\left|f'(x_0)\right|}}\left(\frac{\sqrt{\left|f'(x_0)\right|}}{f(x_0)}\right)'$
$x_1=\displaystyle\frac{f(x_0)}{\sqrt{\left|f'(x_0)\right|}} \frac{2 f'(x) \sqrt{|f'(x)|}} {2 {[f'(x)]}^2 - f(x) f''(x)}$
Podobnie jak w metodzie Newtona staramy się poszukiwania tą metodą przeprowadzić dwa razy, za każdym razem za punkt startowy obierając inny koniec zakresu (a,b). I podobnie jak w metodzie Newtona nie ma gwarancji na znalezienie miejsca zerowego.
public class CentralDifferenceMethod
{
private double m_h;
private Func<double, double> m_func;
public CentralDifferenceMethod()
: base()
{
}
public CentralDifferenceMethod(Func<double, double> a_func,
double a_h = Constants.DOUBLE_PRECISION * 2)
{
m_h = a_h;
m_func = a_func;
}
public double Derive1(double a_x)
{
return (m_func(a_x + m_h) - m_func(a_x - m_h)) / (2 * m_h);
}
public double Derive2(double a_x)
{
return (m_func(a_x + m_h) - 2 * m_func(a_x) +
m_func(a_x - m_h)) / (m_h * m_h);
}
public static Func<double, double> Derive1(
Func<double, double> a_func, double a_h)
{
CentralDifferenceMethod cdm =
new CentralDifferenceMethod(a_func, a_h);
return (x) =>
{
return cdm.Derive1(x);
};
}
public static Func<double, double> Derive2(
Func<double, double> a_func, double a_h)
{
CentralDifferenceMethod cdm =
new CentralDifferenceMethod(a_func, a_h);
return (x) =>
{
return cdm.Derive2(x);
};
}
public static Func<double, double> Derive2(
Func<double, double> a_func,
Func<double, double> a_d1, double a_h)
{
CentralDifferenceMethod cdm =
new CentralDifferenceMethod(a_d1, a_h);
return (x) =>
{
return cdm.Derive1(x);
};
}
}
public abstract class RootFinder
{
public int MaxIterations;
public double AbsoluteError;
public RootFinder()
{
MaxIterations = Constants.ROOT_FINDER_MAXIMUM_ITERATIONS;
AbsoluteError = Constants.ROOT_FINDER_ABSOLUTE_ERROR;
}
public static RootFinder Create(RootFindingMethod a_method)
{
switch (a_method)
{
case RootFindingMethod.Bisection:
return new BisectionRootFinder();
case RootFindingMethod.ModifedRegulaFalsi:
return new ModifiedRegulaFalsiRootFinder();
case RootFindingMethod.RegulaFalsi:
return new RegulaFalsiRootFinder();
case RootFindingMethod.Secant:
return new SecantRootFinder();
case RootFindingMethod.Newton:
return new NewtonRootFinder();
case RootFindingMethod.Halley:
return new HalleyRootFinder();
default: throw new NotImplementedException();
}
}
public abstract double FindRoot(double a_a, double a_b,
Polynomial a_poly);
public abstract double FindRoot(double a_a, double a_b,
Func<double, double> a_func);
public abstract double FindRoot(double a_a, double a_b,
Func<double, double> a_func, Func<double, double> a_d1);
public abstract double FindRoot(double a_a, double a_b,
Func<double, double> a_func, Func<double, double> a_d1,
Func<double, double> a_d2);
}
public abstract class DerivativeRootFinder : RootFinder
{
protected struct RootState
{
public double Root;
public bool Finded;
public RootState(double a_x, bool a_finded = true)
{
Debug.Assert(a_x.IsNumber());
Root = a_x;
Finded = a_finded;
}
}
protected double FindRootFromBothSide(double a_a, double a_b,
Func<double, double> a_func, Func<double,
RootState> a_root_finder)
{
Debug.Assert(a_a < a_b);
Debug.Assert(a_a.IsNumber());
Debug.Assert(a_b.IsNumber());
double ya = a_func(a_a);
double yb = a_func(a_b);
if (Math.Abs(ya) < AbsoluteError)
return a_a;
if (Math.Abs(yb) < AbsoluteError)
return a_b;
if (Math.Sign(ya) * Math.Sign(yb) > 0)
return Double.PositiveInfinity;
RootState x1 = a_root_finder(a_a);
if (x1.Finded && (x1.Root > a_a) && (x1.Root < a_b))
return x1.Root;
RootState x2 = a_root_finder(a_b);
if (x2.Finded && (x2.Root > a_a) && (x2.Root < a_b))
return x2.Root;
return Double.PositiveInfinity;
}
}
public class HalleyRootFinder : DerivativeRootFinder
{
public override double FindRoot(double a_a, double a_b,
Func<double, double> a_func)
{
return FindRoot(a_a, a_b, a_func,
CentralDifferenceMethod.Derive1(a_func, AbsoluteError),
CentralDifferenceMethod.Derive2(a_func, AbsoluteError));
}
public override double FindRoot(double a_a, double a_b,
Func<double, double> a_func, Func<double, double> a_d1)
{
return FindRoot(a_a, a_b, a_func, a_d1,
CentralDifferenceMethod.Derive2(a_func, a_d1, AbsoluteError));
}
public override double FindRoot(double a_a, double a_b,
Func<double, double> a_func, Func<double, double> a_d1,
Func<double, double> a_d2)
{
return FindRootFromBothSide(a_a, a_b, a_func,
(s) => FindRoot(s, a_func, a_d1, a_d2));
}
private RootState FindRoot(double a_x, Func<double, double> a_func,
Func<double, double> a_d1, Func<double, double> a_d2)
{
double prev_y = a_func(a_x);
for (int i = 0; i < MaxIterations; i++)
{
double dy = a_d1(a_x);
double ddy = a_d2(a_x);
double den = 2 * dy * dy - prev_y * ddy;
if (den.IsAlmostEquals(0))
return new RootState(a_x);
a_x = a_x - 2 * prev_y * dy / den;
double y = a_func(a_x);
prev_y = y;
if (Math.Abs(a_x) > AbsoluteError)
{
if (Math.Abs(y / a_x) < AbsoluteError)
return new RootState(a_x);
}
else if (Math.Abs(y) < AbsoluteError)
return new RootState(a_x);
}
return new RootState(a_x, false);
}
public override double FindRoot(double a_a, double a_b,
Polynomial a_poly)
{
Polynomial d1 = a_poly.Differentiate();
Polynomial d2 = d1.Differentiate();
return FindRoot(a_a, a_b, (x) => a_poly.Evaluate(x),
(x) => d1.Evaluate(x), (x) => d2.Evaluate(x));
}
}
Brak komentarzy:
Prześlij komentarz