ODE45和Runge-Kutta方法的绝对误差与解析解的比较 [英] Absolute error of ODE45 and Runge-Kutta methods compared with analytical solution
问题描述
如果有人可以帮助解决以下问题,我将不胜感激. 我有以下ODE:
dr/dt = 4*exp(0.8*t) - 0.5*r ,r(0)=2, t[0,1] (1)
我用两种不同的方法解决了(1).
通过 Runge-Kutta方法(四阶)并通过Matlab中的ode45
.我将这两种结果与解析解进行了比较,解析解由下式给出:
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方法的绝对误差时,为什么它看起来更好?
您的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更高阶的方法,因此您可以期望得到此结果).由于插值,误差在积分点之间增大.如果要限制此范围,则应通过
如果要强制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
不一定是完成此工作的最佳工具.但是,有替代品,但它们可能需要一些编程.
I would appreciate if someone can help with the following issue. I have the following ODE:
dr/dt = 4*exp(0.8*t) - 0.5*r ,r(0)=2, t[0,1] (1)
I have solved (1) in two different ways.
By means of the Runge-Kutta method (4th order) and by means of ode45
in Matlab. I have compared the both results with the analytic solution, which is given by:
r(t) = 4/1.3 (exp(0.8*t) - exp(-0.5*t)) + 2*exp(-0.5*t)
When I plot the absolute error of each method with respect to the exact solution, I get the following:
For RK-method, my code is:
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
And for 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);
My question is, why do I have oscillations when I use ode45
? (I am referring to the absolute error). Both solutions are accurate (1e-9
), but what happens with ode45
in this case?
When I compute the absolute error for the RK-method, why does it looks nicer?
Your RK4 function is taking fixed steps that are much smaller than those that ode45
is taking. What you're really seeing is the error due to polynomial interpolation that is used to produce the points in between the true steps that ode45
takes. This is often referred to as "dense output" (see Hairer & Ostermann 1990).
When you specify a TSPAN
vector with more than two elements, Matlab's ODE suite solvers produce fixed step size output. This does not mean that they that they actually use a fixed step size or that they use the step sizes specified in your TSPAN
however. You can see the actual step sizes used and still get your desired fixed step size output by having ode45
output a structure and using deval
:
sol = ode45(f,tspan,x0);
diff(sol.x) % Actual step sizes used
y_ode45 = deval(sol,tspan);
You'll see that after an initial step of 0.02
, because your ODE is simple it converges to 0.1
for the subsequent steps. The default tolerances combined with the default maximum step size limit (one tenth the integration interval) determine this. Let's plot the error at the true steps:
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)
As you can see, the error at the true steps grows more slowly than the error for RK4 (ode45
is effectively a higher order method than RK4 so you'd expect this). The error grows in between the integration points due to the interpolation. If you want to limit this, then you should adjust the tolerances or other options via odeset
.
If you wanted to force ode45
to use a step of 1/50
you can do this (works because your ODE is simple):
opts = odeset('MaxStep',1/50,'InitialStep',1/50);
sol = ode45(f,tspan,x0,opts);
diff(sol.x)
y_ode45 = deval(sol,tspan);
For another experiment, try enlarging the integration interval to integrate out to t = 10
maybe. You'll see lots of interesting behavior in the error (plotting relative error is useful here). Can you explain this? Can you use ode45
and odeset
to produce results that behave well? Integrating exponential functions over large intervals with adaptive step methods is challenging and ode45
is not necessarily the best tool for the job. There are alternatives however, but they may require some programming.
这篇关于ODE45和Runge-Kutta方法的绝对误差与解析解的比较的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持IT屋!