Как построить график сходимости итерационного метода?

Question

Я начал учиться программировать в MatLab и пытаюсь понять, как быстро сходится метод Ньютона. В настоящее время я хочу построить абсолютную ошибку | xn − x | из итераций xn, где x = sqrt (10) - истинное решение. Ошибка должна стремиться к нулю, поскольку метод сходится, но при его построении точки не отображаются. Может кто-нибудь помочь мне разобраться в этом вопросе?

f = @(x) (x).^2 - 10;
f_deriv = @(x) 2*x;
x0 = 1;
x_true = sqrt(10);
x_save = (x0);

for jj = 1:20
    plot(jj,abs(x_save(jj)-x_true),'k.');
    x_new = x0 - f(x0)/f_deriv(x0);
    x_save(1+jj,1) = x_new;
    x0 = x_new;
end

xlim([0 20]);
ylim([0 2.5]);
xticks(0:20);
yticks(0:0.1:2.5);

Answer 1

Как вы говорите, часть метода Ньютона работает так, как задумано. Проблема с вашей графикой.

В вашей реализации вы генерируете 1 график за итерацию в цикле for. Это приводит к 20 участкам с 1 очком за участок, что не является желаемым результатом.

Ниже приведен код для генерации 1 графика с 21 точкой.

f = @(x) (x).^2 - 10;
f_deriv = @(x) 2*x;
x0 = 1;
x_true = sqrt(10);
x_save = (x0);

for jj = 1:20
    % Don't plot here. Wait until all data is collected.
    x_new = x0 - f(x0)/f_deriv(x0);
    x_save(1+jj,1) = x_new;
    x0 = x_new;
end

% Plot here. All data has been collected.
plot(0:20',abs(x_save-x_true),'kx');
xlim([0 20]);
ylim([0 2.5]);
xticks(0:20);
yticks(0:0.1:2.5);

---

РЕДАКТИРОВАТЬ: Для решения проблемы пропущенных точек с semilogy вместо plot.

1. Точки отсутствуют, поскольку semilogy вычисляет журнал значений. log(0) = -Inf
2. Точность станка приводит к тому, что очень маленькие числа рассматриваются как ноль.
3. Точность вычитания с плавающей запятой пропорциональна величине входных данных. Вычитание двух больших чисел имеет низкую точность.

Ниже приведен пример изменения координат, чтобы избежать вычитания двух больших чисел.

% Original problem parameters
x0 = 1;
x_true = sqrt(10);

% Change of co-ordinates
% u = x - sqrt(10);
% f(x) = g(u(x))
g = @(u) u^2 + 2*x_true*u;
g_deriv = @(u) 2*u + 2*x_true;
u0 = x0 - x_true;
u_save = u0;

for jj = 1:20
    u_new = u0 - g(u0)/g_deriv(u0);
    u_save = [u_save u_new];
    u0 = u_new;
end

semilogy(0:20, abs(u_save - 0), 'kx')
xlim([0 20]);

Вы заметите, что точки за 6 итераций по-прежнему исчезают, потому что ошибка все еще меньше точности машины (~ 1e-15). Однако ошибка снова не увеличивается.

Если вы хотите повысить точность своих расчетов, вам может понадобиться функция vpa.

https://www.mathworks.com/help/symbolic/increase-precision-of-numeric-calculations.html

Answer 2

Вам нужно только hold on, чтобы каждый график добавлялся к предыдущим. По умолчанию каждый новый сюжет заменяет предыдущий.

Кроме того, я изменил маркер с '.' до 'o' для лучшей видимости.

Вы также можете включить drawnow (и, возможно, подходящий pause) после каждого plot, чтобы цифра была немедленно обновлена. Это полезно, если вы хотите увидеть линию, обозначенную маркерами «расти».

f = @(x) (x).^2 - 10;
f_deriv = @(x) 2*x;
x0 = 1;
x_true = sqrt(10);
x_save = (x0);

hold on %%% New line
for jj = 1:20
    plot(jj,abs(x_save(jj)-x_true),'ko'); %%% Modified line
    x_new = x0 - f(x0)/f_deriv(x0);
    x_save(1+jj,1) = x_new;
    x0 = x_new;
end

xlim([0 20]);
ylim([0 2.5]);
xticks(0:20);
yticks(0:0.1:2.5);