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

Буду признателен, если кто-то поможет со следующей проблемой. У меня есть следующий 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)

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

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

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?

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


matlab ode differential-equations numerical-integration runge-kutta
person Sergio Haram    schedule 18.02.2014    source источник

Ответы (2)


arrow_upward
22
arrow_downward

Ваша функция 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, поскольку ваше ОДУ простое, оно сходится к 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 не обязательно является лучшим инструментом для этой работы. Однако существуют альтернативы, но они могут потребовать некоторого программирования.

person horchler    schedule 18.02.2014
comment
Привет @horchler. Жаль, что я могу поставить только один балл вашему ответу. На данный момент я не уверен, чему я больше завидую: вашим навыкам программирования или вашему математическому пониманию. - person Sergio Haram; 19.02.2014
comment
Я не знал, что на самом деле шаги от ode45 где точки снизу. Если это так, то я с вами согласен на 100%, ода точнее. Не могли бы вы принять «тик-так» к обоим методам? Я обнаружил, что метод RK использует 0,018 секунды, а ode45 использует 0,5 секунды. Могу ли я сделать вывод, что ode45 точнее, но медленнее? И я полагал, что адаптивный контроллер размера шага от ode45 очень хорош(?) - person Sergio Haram; 19.02.2014
comment
Я еще не пытался увеличить интервал интегрирования, но мне очень любопытен результат. В то же время меня заинтриговали ваши комментарии о том, каких хороших результатов я могу ожидать от ode45 и odeset. Я попробую это как сын, как я могу! Спасибо, что поделился! - person Sergio Haram; 19.02.2014
comment
@SergioHaram: переоценка времени. ode45 требует некоторой инициализации. Например, необходимо выяснить, какие размеры шагов использовать. Это требует времени. Вы можете дать ему подсказку через opts = odeset('InitialStep',0.1);, а затем передать opts as the last argument to ode45. Also, because the ODE is so simple (at least it looks simple to the integrator, but as I explain at the end exponential growth can be challenging) you could try using ode23`. Более простые методы с фиксированным размером шага могут быть быстрее во многих случаях, но обычно не в тех случаях, когда ОДУ более сложное, например, много генераторов. - person horchler; 19.02.2014
comment
@ZheyuanLi: я бы посоветовал отправиться в чат-комнату Matlab для дальнейшего обсуждения/вопросов по это. Там много тех, кто может вам помочь. - person horchler; 09.12.2016

arrow_upward
0
arrow_downward

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

person EngrStudent    schedule 18.02.2014
comment
Это не то, что на самом деле показывает график ошибок. ode45 в этом случае сходится к постоянному размеру шага вторым шагом. И график также не показывает, что ошибка ограничена, она просто растет медленнее, как и следовало ожидать. - person horchler; 18.02.2014
comment
Вы правы, и я дал положительную оценку вашему ответу за это. Более точный способ, и тот, который я должен был использовать в описании, убегает быстрее, а частота ошибок остается более ограниченной. - person EngrStudent; 18.02.2014