Решение нелинейного уравнения. Метод хорд и касательных - C++ - Ответ 1578627

typedef double(*foo)(double,double);
 
typedef struct _item
{
    double a,b,p,x;         //[a,b] - интервал, на котором есть корень для уравнения с параметром p. Решение - x
}item;
 
double f(double x, double p)
{
    return log(x+1)-p/(x*x);
}
 
class HK
{
    private:
        foo f;              //функция по которой считаем
        double c,d;         //границы
        int n;              //на сколько разбиваем
        double eps;
        vector<item> s;
        double df(double x, double p);
        double ddf(double x, double p);
        item m(vector<item> x, int f=-1);   //f=1 - максимум, f=-1 - минимум
    public:
        HK(foo function, double c=-1, double d=1, int n=20, double e=0.01);
        vector<item> solve();
        double clarify(double a, double b, double p);
};
 
    HK::HK(foo function, double c, double d, int n, double e)
    {
        f=function;
        this->c=c;
        this->d=d;
        this->n=n;
        eps=e;
        s.clear();
    }
 
    double HK::ddf(double x, double p)
    {
        double dx=pow(10.0,-6);
        return (((f)(x,p)-2*(f)(x-dx,p)+f(x-2*dx,p))/(pow(dx,2)));
    }
 
    double HK::df(double x, double p)
    {
        double dx=pow(10.0,-6);
        return (((f)(x+dx,p)-(f)(x,p))/dx);
    }
 
    item HK::m(vector<item> x, int f)
    {
        item r;
        r.x=(-f)*32000;
        for(int i=0;i<x.size();i++)
            if(f*x[i].x>f*r.x) r=x[i];
        return r;
    }
 
    double HK::clarify(double a, double b, double p)
    {
        double x=(a+b)/2;
        while(fabs(a-b)>eps)
        {
            double ai=a, bi=b;      //старые значения - нужны для расчета новых
            if(df(x,p)*ddf(x,p)>0)
            {
                a=ai-((f)(ai,p)*(bi-ai))/((f)(bi,p)-f(ai,p));
                b=bi-((f)(bi,p))/(df(bi,p));
            }
            else
            {
                a=ai-((f)(ai,p))/(df(ai,p));
                b=bi-((f)(bi,p)*(bi-ai))/((f)(bi,p)-f(ai,p));
            }
            x=(a+b)/2;
        }
        return x;
    }
 
    vector<item> HK::solve()
    {
        double h=(d-c)/n;           //шаг p
        vector<item> x_all;     //вектор всех корней
        for(double p=c;p<=d;p+=h)   //для каждого p от с до d
        {
            x_all.clear();
            for(double x=-2;x<=2;x+=0.2)//бесконечный цикл нахождения корня, начиная от -2. Ищем интервалы по 0.2 и корень в нём
            {
                double a=(f)(x,p), b=(f)(x+0.2,p), a_d=df(x,p), b_d=df(x+0.2,p);
                if((fabs(a-b)>fabs(a+b))&&(fabs(a_d-b_d)<fabs(a_d+b_d)))
                //если знаки разные, то корень есть, а если производная не меняет свой знак - корень один
                //если на интервале [x;x+0.2] есть корень, то его нужно уточнить методом хорд и касательных
                //и занести в значения x_all
                {
                    item temp;
                    temp.a=x;
                    temp.b=x+0.2;
                    temp.p=p;
                    temp.x = clarify(a,b,p);
                    x_all.push_back(temp);
                }
            }
            //если корень не один, то отдаем предпочтение положительным. Если все корни одного знака - наибольшего по модулю
            if(x_all.size()>1)
            {
                vector<item> less, large;
                less.clear(); large.clear();
                for(int i=0;i<x_all.size();i++)
                {
                    if(x_all[i].x<0) less.push_back(x_all[i]);
                    else large.push_back(x_all[i]);
                }
                x_all.clear();
                if(!large.size())
                    x_all.push_back(m(less,-1));
                else
                    x_all.push_back(m(large,1));
            }
            s.push_back(x_all[0]);
        }
        return s;
    }