Абсолютная погрешность методов ODE45 и Рунге-Кутты по сравнению с аналитическим решением

Question

Буду признателен, если кто-то может помочь со следующей проблемой. У меня есть следующие ODE:

dr/dt = 4*exp(0.8*t) - 0.5*r   ,r(0)=2, t[0,1]       (1)

Я решил (1) двумя разными способами. С помощью метода Рунге-Кутты (4-й порядок) и с помощью ode45 в Matlab. Я сравнил оба результата с аналитическим решением, которое дает:

r(t) = 4/1.3 (exp(0.8*t) - exp(-0.5*t)) + 2*exp(-0.5*t)

Когда я строю абсолютную ошибку каждого метода относительно точного решения, я получаю следующее:

Для RK-метода мой код:

h=1/50;                                            
x = 0:h:1;                                        
y = zeros(1,length(x)); 
y(1) = 2;    
F_xy = @(t,r) 4.*exp(0.8*t) - 0.5*r;                   
for i=1:(length(x)-1)                              
    k_1 = F_xy(x(i),y(i));
    k_2 = F_xy(x(i)+0.5*h,y(i)+0.5*h*k_1);
    k_3 = F_xy((x(i)+0.5*h),(y(i)+0.5*h*k_2));
    k_4 = F_xy((x(i)+h),(y(i)+k_3*h));
    y(i+1) = y(i) + (1/6)*(k_1+2*k_2+2*k_3+k_4)*h;  % main equation
end

И для ode45:

tspan = 0:1/50:1;
x0 = 2;
f = @(t,r) 4.*exp(0.8*t) - 0.5*r;
[tid, y_ode45] = ode45(f,tspan,x0);

У меня вопрос, почему у меня есть колебания, когда я использую ode45? (Я имею в виду абсолютную ошибку). Оба решения точны (1e-9), но что происходит с ode45 в этом случае?

Когда я вычисляю абсолютную ошибку для RK-метода, почему она выглядит лучше?

Answer 1

Ваша функция RK4 принимает фиксированные шаги, которые намного меньше, чем те, которые ode45 принимает. В действительности вы видите ошибку из-за полиномиальной интерполяции, которая используется для получения точек между истинными шагами, которые ode45 принимает. Это часто называют "плотным выходом" (см. Hairer & Ostermann 1990).

Когда вы указываете TSPAN вектор с более чем двумя элементами, решатели пакета ODE Matlab производят вывод с фиксированным размером шага. Это не означает, что они на самом деле используют фиксированный размер шага или что они используют размеры шагов, указанные в вашем TSPAN тем не мение. Вы можете увидеть фактические используемые размеры шагов и по-прежнему получать желаемый фиксированный размер шага, имея ode45 вывести структуру и использовать deval:

sol = ode45(f,tspan,x0);
diff(sol.x) % Actual step sizes used
y_ode45 = deval(sol,tspan);

Вы увидите это после начального шага 0.02 потому что ваш ODE прост, он сходится к 0.1 для последующих шагов. Допуски по умолчанию в сочетании с пределом максимального размера шага по умолчанию (одна десятая интервала интеграции) определяют это. Давайте нарисуем ошибку на истинных шагах:

exactsol = @(t)(4/1.3)*(exp(0.8*t)-exp(-0.5*t))+2*exp(-0.5*t);
abs_err_ode45 = abs(exactsol(tspan)-y_ode45);
abs_err_ode45_true = abs(exactsol(sol.x)-sol.y);
abs_err_rk4 = abs(exactsol(tspan)-y);
figure;
plot(tspan,abs_err_ode45,'b',sol.x,abs_err_ode45_true,'k.',tspan,abs_err_rk4,'r--')
legend('ODE45','ODE45 (True Steps)','RK4',2)

Как видите, ошибка на истинных шагах растет медленнее, чем ошибка для RK4 (ode45 по сути, метод более высокого порядка, чем RK4, так что вы ожидаете этого). Ошибка растет между точками интегрирования из-за интерполяции. Если вы хотите ограничить это, то вам следует настроить допуски или другие параметры с помощью odeset,

Если ты хотел заставить ode45 использовать шаг 1/50 Вы можете сделать это (работает, потому что ваш ODE прост):

opts = odeset('MaxStep',1/50,'InitialStep',1/50);
sol = ode45(f,tspan,x0,opts);
diff(sol.x)
y_ode45 = deval(sol,tspan);

Для другого эксперимента попробуйте увеличить интервал интегрирования, чтобы интегрировать в t = 10 может быть. Вы увидите много интересного в этой ошибке (здесь может быть показана относительная ошибка). Вы можете это объяснить? Вы можете использовать ode45 а также odeset получить результаты, которые ведут себя хорошо? Интеграция экспоненциальных функций через большие интервалы с адаптивными методами шага является сложной и ode45 не обязательно лучший инструмент для работы. Однако есть альтернативы, но они могут потребовать некоторого программирования.

Answer 2

ode45 связан rk4-rk5. Лично я думаю, что ошибка ODE45 приятнее. Обратите внимание, что он остается ограниченным. Ode4 исправляется, когда величина ошибки становится слишком большой, а минимальная ошибка за цикл составляет около 1e-10. Rk4 "убегает" и ничто не мешает ему.