1

我看到了这个计算错误的欧拉方法：

function [t,le,ge] = euler_errors(h)

f=@(u) u*(2-u); % this is the function for the IVP
t0=0;
tn=1;
t=t0:h:tn;%we want to find the errors along this solutions
%here is the exact solution of the IVP
u_exact=(0.2*exp(2*t))./(2+0.1*(exp(2*t)+1)); %the with initial value u(0)=0.1

n=length(t);
u_e=zeros(1,n);
u_g=zeros(1,n);

i=1;
u_e(i)=0.1;
u_g(i)=0.1;

%u_e and u_g are both values given by Euler method
%u_e is for the local error
%u_g is for the global error
while (i<n)
    u_e(i+1)=u_e(i)+h*f(u_e(i));
    u_g(i+1)=u_g(i)+h*f(u_exact(i));
    i=i+1;
end;
%le1 is the local error
%ge1 is the global error
le=abs(u_e-u_exact);
ge=abs(u_g-u_exact);
end

我试图将该方法转换为 Heun 的方法——这是我的尝试：

function [t,le,ge] = heun_errors(h)

f=@(u) u*(2-u); % this is the function for the IVP
t0=0;
tn=1;
t=t0:h:tn;%we want to find the errors along this solutions
%here is the exact solution of the IVP
u_exact=(0.2*exp(2*t))./(2+0.1*(exp(2*t)+1)); %the with initial value u(0)=0.1

n=length(t);
u_e=zeros(1,n);
u_g=zeros(1,n);

i=1;
u_e(i)=0.1;
u_g(i)=0.1;

%u_e and u_g are both values given by Euler method
%u_e is for the local error
%u_g is for the global error
while (i<n)
    u_e1(i+1)=u_e(i)+h*f(u_e(i));
u_e(i+1)=u_e(i)+(h/2)*(f(u_e(i))+f(u_e1(i+1)));
u_g1(i+1)=u_g(i)+h*f(u_exact(i));
u_g(i+1)=u_g(i)+(h/2)*(f(u_exact(i))+f(u_g1(i+1)));

    i=i+1;
end;
%le1 is the local error
%ge1 is the global error
le=abs(u_e-u_exact);
ge=abs(u_g-u_exact);
end

但现在误差实际上增加了。有人可以告诉我我做错了什么以及如何解决吗？

4

1 回答 1

0

精确解的修正

首先，精确解有一个符号错误，它使所有错误计算无效。作为伯努利方程，可以通过以下方式获得解

(e^(2t)/u)'=e^(2t)*(-u'/u^2 + 2/u) = e^(2t)

e^(2t)/u(t)-1/u0 = (e^(2t)-1)/2  

u(t) = 2*u0*exp(2*t) ./ ( 2 + u0*(exp(2*t)-1) )

关于局部和全局误差的计算

用于计算局部和全局误差的欧拉代码的组织在概念上已经混淆了一些东西。完整的正常积分给出了全局结果，因此也给出了精确解的全局误差。因此所谓的u_e应该是u_g。

局部误差将开始的精确解与从同一点u_i(t)开始的方法步骤进行比较。(t_i,u_i)由于只有一个精确的解决方案，它应该执行步骤 for(t_i,u_i) = (t(i), u_exact(i))然后将存储结果的方法步骤u_e(i+1)与精确值进行比较u_exact(i+1)。

两者都通过改变循环来解决

while (i<n)
  u_g(i+1)=u_g(i)+h*f(u_g(i));
  u_e(i+1)=u_exact(i)+h*f(u_exact(i));
  i=i+1;
end;

类似地，这些问题也必须在 Heun 代码中解决，因此循环更改为（在一般 Runge-Kutta 方法中使用阶段的中间变量）

while (i<n)
  k1=f(u_g(i));
  k2=f(u_g(i)+h*k1);
  u_g(i+1)=u_g(i)+(h/2)*(k1+k2);
  k1=f(u_exact(i));
  k2=f(u_exact(i)+h*k1);
  u_e(i+1)=u_exact(i)+(h/2)*(k1+k2);
  i=i+1;
end;

错误曲线图

然后绘制le/h^(p+1)局部和lg/h^p全局误差随时间的误差曲线，使用p=1Euler 和p=2Heun 的误差阶数，并覆盖不同步长的点图，给出

误差曲线图

它们在视觉上相对于步长是不变的。

作为一致性检查，局部 Euler 误差在 左右0.3*h^2，而 ODE 的 Lipschitz 常数在相关区域中在 2 左右，因此全局误差遵循 Grönwall 引理近似于(exp(2*t)-1)*0.15*h小的定律t。类似地，Heun 方法在给small0.15*h^3的全局误差附近有局部误差，这两者都与全局误差图兼容。(exp(2*t)-1)*0.075*h^2t

对于较大t的 ，局部误差的符号可以改变，从而补偿部分早期的全局误差。然而，这并不能保证会发生，也不容易预测。

于 2019-01-01T17:31:13.400 回答