---

% ============================================================
% РЕШЕНИЕ ДИФФЕРЕНЦИАЛЬНОГО УРАВНЕНИЯ
% dy/dx = -2*x*y, y(0)=1, x ∈ [0, 2]
% 
% Лекция 12: "Решение дифференциальных уравнений"
% ============================================================

clear; clc; close all;

%% ЧАСТЬ 1: АНАЛИТИЧЕСКОЕ РЕШЕНИЕ (функция dsolve)
% В лекции: syms y(x); dy = diff(y,x) == -2*x*y; y0 = y(0)==1; dsolve(dy, y0)

disp('=' * 60);
disp('ЧАСТЬ 1: АНАЛИТИЧЕСКОЕ РЕШЕНИЕ (dsolve)');
disp('=' * 60);

syms y(x)                    % объявляем y как символьную функцию от x
dy = diff(y, x) == -2*x*y;   % определяем уравнение: dy/dx = -2*x*y
y0 = y(0) == 1;              % начальное условие: при x=0, y=1

y_analytical = dsolve(dy, y0);  % точное аналитическое решение
disp(' ');
disp(['Точное решение: y(x) = ', char(y_analytical)]);
disp(' ');

% Преобразуем символьное решение в анонимную функцию для вычислений
y_func = matlabFunction(y_analytical);

%% ЧАСТЬ 2: ЧИСЛЕННОЕ РЕШЕНИЕ (функция ode45)
% В лекции: dydx = @(x,y) -2*x*y; [x, y] = ode45(dydx, [0 2], 1)

disp('=' * 60);
disp('ЧАСТЬ 2: ЧИСЛЕННОЕ РЕШЕНИЕ (ode45)');
disp('=' * 60);

% Определяем правую часть уравнения (как в лекции, стр. 2-3)
dydx = @(x, y) -2 * x * y;   % формируем уравнение: dy/dx = -2*x*y
x_interval = [0, 2];          % интервал интегрирования
y0 = 1;                       % начальное условие y(0)=1

% Решаем ДУ с помощью ode45 (как в лекции, стр. 3)
[x_numerical, y_numerical] = ode45(dydx, x_interval, y0);

disp(['Интервал интегрирования: [0, 2]']);
disp(['Начальное условие: y(0) = ', num2str(y0)]);
disp(['Количество точек решения: ', num2str(length(x_numerical))]);
disp(' ');
disp(['Значение в конце интервала: y(2) = ', num2str(y_numerical(end))]);

%% ЧАСТЬ 3: СРАВНЕНИЕ РЕЗУЛЬТАТОВ И ВИЗУАЛИЗАЦИЯ
% В лекции: fplot(@(t) exp(-t.^2), [0 2]); hold on; plot(x, y, 'ro')

disp('=' * 60);
disp('ЧАСТЬ 3: СРАВНЕНИЕ РЕЗУЛЬТАТОВ');
disp('=' * 60);

% Создаем более плотную сетку для аналитического графика
x_fine = linspace(0, 2, 200);
y_fine = y_func(x_fine);

% Вычисляем погрешность численного метода
% Интерполируем численное решение на ту же сетку, что и аналитическое
y_numerical_interp = interp1(x_numerical, y_numerical, x_fine);
error = max(abs(y_fine - y_numerical_interp));
disp(['Максимальная погрешность ode45: ', num2str(error), ' (должна быть малой)']);

%% ПОСТРОЕНИЕ ГРАФИКОВ (как в лекции, стр. 3)

figure('Position', [100, 100, 800, 600]);

% График 1: Сравнение решений
subplot(2, 1, 1);
fplot(y_analytical, [0, 2], 'b-', 'LineWidth', 2);  % аналитическое решение
hold on;
plot(x_numerical, y_numerical, 'ro', 'MarkerSize', 6, 'MarkerFaceColor', 'r');  % численное
hold off;
grid on;
xlabel('x', 'FontSize', 12);
ylabel('y', 'FontSize', 12);
title('Решение дифференциального уравнения dy/dx = -2*x*y', 'FontSize', 14);
legend('Точное решение y = e^{-x^2}', 'Численное решение (ode45)', 'Location', 'best');

% График 2: Погрешность
subplot(2, 1, 2);
error_plot = abs(y_fine - y_numerical_interp);
plot(x_fine, error_plot, 'k-', 'LineWidth', 1.5);
grid on;
xlabel('x', 'FontSize', 12);
ylabel('|y_точн - y_числ|', 'FontSize', 12);
title('Погрешность численного решения', 'FontSize', 14);

% Сохраняем график
saveas(gcf, 'ODE_solution.png');

%% ДОПОЛНИТЕЛЬНО: ДРУГИЕ РЕШАТЕЛИ (из лекции, стр. 5-6)
% В лекции описаны разные решатели:
% ode45 - методы Рунге-Кутта 4-5 порядка (по умолчанию)
% ode23 - методы Рунге-Кутта 2-3 порядка (для низкой точности)
% ode113 - метод Адамса-Башворта-Мултона (для высокой точности)
% ode15s - для жестких задач

disp(' ');
disp('=' * 60);
disp('ДОПОЛНИТЕЛЬНО: ДРУГИЕ РЕШАТЕЛИ');
disp('=' * 60);

% Решаем тем же уравнением с разными методами
[x23, y23] = ode23(dydx, x_interval, y0);
[x113, y113] = ode113(dydx, x_interval, y0);

% Сравниваем количество точек (эффективность методов)
disp(['ode45: ', num2str(length(x_numerical)), ' точек']);
disp(['ode23: ', num2str(length(x23)), ' точек (меньше, но точность ниже)']);
disp(['ode113: ', num2str(length(x113)), ' точек (адаптивный, может быть эффективнее)']);

%% ВЫВОД ИТОГОВ В КОНСОЛЬ
disp(' ');
disp('=' * 60);
disp('ИТОГОВЫЕ РЕЗУЛЬТАТЫ');
disp('=' * 60);
disp(['Аналитическое решение: y(x) = exp(-x^2)']);
disp(['Значение y(2): ', num2str(y_func(2))]);
disp(['Численное решение y(2): ', num2str(y_numerical(end))]);
disp(['Погрешность в точке x=2: ', num2str(abs(y_func(2) - y_numerical(end)))]);