与解析解相比,ODE45和Runge-Kutta方法的绝对误差

Question

如果有人可以帮助解决以下问题,我将不胜感激.我有以下ODE:

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

我用两种不同的方式解决了(1).通过Runge-Kutta方法(第4顺序)和ode45Matlab中的方法.我将这两个结果与分析解决方案进行了比较,分析解决方案由下式给出:

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两个以上元素的向量时,Matlab的ODE套件求解器会生成固定步长输出.这并不意味着他们实际上使用固定步长或者他们使用您指定的步长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不一定是工作的最佳工具.但是有其他选择,但它们可能需要一些编程.