9

如果有人可以帮助解决以下问题，我将不胜感激。我有以下 ODE：

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

我以两种不同的方式解决了（1）。通过Runge-Kutta 方法（四阶）和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 方法的绝对误差时，为什么它看起来更好？

4

2 回答 2

22

您的 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可能积分。你会在错误中看到很多有趣的行为（在这里绘制相对错误很有用）。你能解释一下吗？你能使用ode45andodeset产生表现良好的结果吗？使用自适应步长方法在大间隔内集成指数函数具有挑战性，并且ode45不一定是这项工作的最佳工具。然而，还有其他选择，但它们可能需要一些编程。

于 2014-02-18T18:04:41.187 回答
0

ode45 耦合 rk4-rk5。我个人认为 ODE45 错误更好。请注意，它保持有界。ode4 在误差幅度过大时得到纠正，每个周期的最小误差约为 1e-10。rk4 正在“逃跑”，没有什么能阻止它。

于 2014-02-18T16:51:36.210 回答