Решение задачи Коши усовершенствованным методом Эйлера с автоматическим выбором шага - C (СИ)

Формулировка задачи:

Есть задача: Решить задачу Коши методом Эйлера второго порядка (усовершенствованный метод Эйлера), в программе должен быть автоматический подбор шага. Вот само задание: xy' + y = y*y*ln(x) начальное условие: 0.5; а = 1; b = 2; точное решение: 1 / (1 + х + ln(x)) формулы метода: y_k+1 = y_k + h*f(x_k + h/2, y_k+1 + k1); k1 = h/2 * f(x_k, yp) Алгоритм решения с автоматическим выбором шага: 1. По значению y_k в точке x_k по формуле метода решение y1 в точке x_k + h с шагом h; 2. По значению y_k в точке x_k по формуле метода решение y_пром в точке x_k + h/2 с шагом h/2; 3. по полученному решению y_пром вычисляем по тем же формулам решение у2 в точке x_k + h с шагом h/2; 4. Проверяем условие неравенства g < eps, где g = |y1-y2| при |y1| <= 1 и g = |y1-y2|/|y1| при |y1| > 1; 5. Если неравенство не выполняется, делим шаг пополам и проводим пересчет с с пункта 1, начинаю с точки x_k 6. если получили нужную точность, тоесть неравенство выполняется, то за решение в точке x_k+h берем величину y_k+1 = y2 - (y1 - y2)/(2^s - 1), где s - порядок метода (то есть у меня s = 2); 7. выбираем величину нового шага по формулам: alpha = 0.8 * pow(eps/g, 1/(s+1)), h_новы = alpha * h_стар 8. преходим к следующему шагу начиная с пункта 1 вот мой код, пробовал в него подставлять обычный метод Эйлера, метод Хойне и метод Рунге-Кутта 4-го порядка, постоянно погрешность до 3%. что делать уже не знаю( если можете помогите.

#include "stdafx.h"
#include <stdio.h>
#include <conio.h>
#include <math.h>
#include <stdlib.h>
 
double ft(double x, double y)
{
    return 1 / (1 + x + log(x));
}
 
double f (double x, double y)
{
    return (y*y * log(x) - y) / x;
}
 
double formula(double x, double y, double h)
{
    double k = h/2 * f(x, y);
    return y + h * f(x + h/2, y + k);
}

void main()
{
    double a = 1.;
    double b = 2.;
    double yk = 0.5;
    double x = a;
    double eps = 0.001;
    double MinStep = (b - a) / 1000000;
    double MaxStep = (b - a) / 4;
    double Step = (b - a) / 1000;
    double g = 0;
    double alpha = 0;
    double y1, y2, yp;
    int s = 2;
    printf("   x  |   y   | f(x) |   D   |  D/f  |   h   |   g   |   a   \n");
    do
    {
        y1 = formula(x, yk, Step);
        yp = formula(x, yk, Step/2);
        y2 = formula(x, yp, Step/2);
        if (fabs(y1) <= 1) g = fabs(y1-y2);
        else{ g = fabs(y1 - y2) / fabs(y1); }
        if (g < eps)
        {
            yk = y2 - (y1 - y2)/(pow(2., s) - 1);
            double D = fabs(ft(x, yk) - yk);
            printf("%5.3lf | %5.3lf | %5.3lf | %5.3lf | %5.3lf | %7.5lf | %5.3lf | %5.3lf\n", 
                x, yk, ft(x, yk), D, D * 100 / ft(x, yk), Step, g, alpha );
 
            alpha = 0.8 * pow(eps / g, 1. / (s+1));
            Step *= alpha;
 
            if (Step < MinStep) { Step = MinStep;   }
            else
            {
                if(Step > MaxStep) Step = MaxStep;
            }
 
            x += Step;
        }
        else
        { 
            Step *= 0.5; 
            if (Step < MinStep) { Step = MinStep;   }
        }
        
    } while (x <= b+Step);
    system("pause");
}

Решение задачи: «Решение задачи Коши усовершенствованным методом Эйлера с автоматическим выбором шага»

Листинг программы

#include "stdafx.h"
#include <stdio.h>
#include <conio.h>
#include <math.h>
#include <stdlib.h>
 
double ft(double x, double y)
{
    return 1 / (1 + x + log(x));
}
 
double f (double x, double y)
{
    return (y*y * log(x) - y) / x;
}
 
double formula(double x, double y, double h)
{
    double k = h/2 * f(x, y);
    return y + h * f(x + h/2, y + k);
}
 
 
void main()
{
    double a = 1.;
    double b = 2.;
    double yk = ft(a);
    double x = a;
    double eps = 0.001;
    double MinStep = (b - a) / 1000000;
    double MaxStep = (b - a) / 4;
    double Step = (b - a) / 1000;
    double g = 0;
    double alpha = 0;
    double y1, y2, yp;
    int s = 2;
    printf("   x  |   y   | f(x) |   D   |  D/f  |   h   |   g   |   a   \n");
    do
    {
        y1 = formula(x, yk, Step);
        yp = formula(x, yk, Step/2);
        y2 = formula(x, yp, Step/2);
        if (fabs(y1) <= 1) g = fabs(y1-y2);
        else{ g = fabs(y1 - y2) / fabs(y1); }
        if (g < eps)
        {
            yk = y2 - (y1 - y2)/(pow(2., s) - 1);
            double D = fabs(ft(x, yk) - yk);
            printf("%5.3lf | %5.3lf | %5.3lf | %5.3lf | %5.3lf | %7.5lf | %5.3lf | %5.3lf\n", 
                x, yk, ft(x, yk), D, D * 100 / ft(x, yk), Step, g, alpha );
 
            alpha = 0.8 * pow(eps / g, 1. / (s+1));
 
            x += Step;
            Step *= alpha;
 
            if (Step < MinStep) { Step = MinStep;   }
            else
            {
                if(Step > MaxStep) Step = MaxStep;
            }
 
        }
        else
        { 
            Step *= 0.5; 
            if (Step < MinStep) { Step = MinStep;   }
        }
        
    } while (x <= b+Step);
    system("pause");
}

Объяснение кода листинга программы

Включаемые заголовочные файлы:
- stdafx.h
- stdio.h
- conio.h
- math.h
- stdlib.h
Функция ft(double x, double y) возвращает значение функции 1 / (1 + x + log(x))
Функция f(double x, double y) возвращает значение функции (y*y * log(x) - y) / x
Функция formula(double x, double y, double h) использует метод Эйлера для решения задачи Коши. Возвращаемое значение функции является приближенным значением корня.
В функции main() определены следующие переменные:
- a (начальное приближение)
- b (конечное приближение)
- yk (начальное значение функции)
- x (текущее значение переменной)
- eps (требуемая точность)
- MinStep (минимальный шаг)
- MaxStep (максимальный шаг)
- Step (текущий шаг)
- g (погрешность на предыдущем шаге)
- alpha (коэффициент уменьшения шага)
- y1, y2, yp (временные переменные для хранения промежуточных значений)
- s (порядок метода Эйлера)
В цикле do-while происходит итеративное решение задачи Коши методом Эйлера. На каждой итерации выполняются следующие действия:
- Вычисляются значения y1, yp и y2 с помощью функции formula()
- Вычисляется погрешность g на предыдущем шаге
- Если погрешность g меньше требуемой точности eps, выполняются следующие действия:
  - Обновляется значение yk
  - Вычисляется значение функции ft(x, yk)
  - Вычисляется значение функции f(x, yk)
  - Вычисляется значение D (погрешность на предыдущем шаге)
  - Вычисляется значение h (новый шаг)
  - Обновляется значение g
  - Обновляется значение alpha
  - Увеличивается значение x на текущий шаг Step
  - Уменьшается значение Step на alpha
  - Если Step меньше MinStep, устанавливается Step равным MinStep
  - Если Step больше MaxStep, устанавливается Step равным MaxStep
- Если погрешность g больше или равна требуемой точности eps, выполняются следующие действия:
  - Уменьшается значение Step в два раза
В конце программы вызывается функция system(pause), чтобы приостановить выполнение программы до нажатия клавиши.