计算科学 - 误差计算 - 吾爱随笔录

误差计算

计算科学 matlab 有限差分数值分析误差估计椭圆pde

2021-12-18 07:13:19

我编写了一个代码，在其中我找到了这个椭圆问题的解的近似值。

我使用以下代码部分计算了错误：

但对于出现以下错误： $a=-1, b=1, k_1=k_2=1$

For N=10:  er = 9.9920e-016
For N=20: er = 2.4425e-015
For N=30: er = 7.1054e-015
For N=40: er = 7.7716e-015
For N=50: er = 1.5765e-014
For N=60: er = 1.6764e-014
For N=70: er = 1.0436e-014

但是随着的增加，误差应该会减小。我在上面的代码中做错了什么？ $N$

编辑：使用此代码： http: //pastebin.com/crS4vb1t

我得到以下结果：

[率，错误]=订单

率 =

0.2221   -2.0880   -0.2637   -2.7620

错误 =

1.0e-014 *

0.0281    0.0241    0.1024    0.1230    0.8342

编辑 2：

对于，，矩阵 A 以及我编写的代码如下： $N=4$ $k_1=k_2=1$

一个=

第 1 至 14 列

-4     1     0     0     1     0     0     0     0     0     0     0     0     0
 1    -4     1     0     0     1     0     0     0     0     0     0     0     0
 0     1    -4     1     0     0     1     0     0     0     0     0     0     0
 0     0     1    -4     0     0     0     1     0     0     0     0     0     0
 1     0     0     0    -4     1     0     0     1     0     0     0     0     0
 0     1     0     0     1    -4     1     0     0     1     0     0     0     0
 0     0     1     0     0     1    -4     1     0     0     1     0     0     0
 0     0     0     1     0     0     1    -4     0     0     0     1     0     0
 0     0     0     0     1     0     0     0    -4     1     0     0     1     0
 0     0     0     0     0     1     0     0     1    -4     1     0     0     1
 0     0     0     0     0     0     1     0     0     1    -4     1     0     0
 0     0     0     0     0     0     0     1     0     0     1    -4     0     0
 0     0     0     0     0     0     0     0     1     0     0     0    -4     1
 0     0     0     0     0     0     0     0     0     1     0     0     1    -4
 0     0     0     0     0     0     0     0     0     0     1     0     0     1
 0     0     0     0     0     0     0     0     0     0     0     1     0     0

第 15 至 16 列

编辑：这是与相关的错误的日志日志图： $N$

所以错误不可能是正确的。还是我错了？

1个回答

经过几天的讨论，您的问题是

1）您选择了过于简单的问题，导致机器为零，因此无法观察到正常属性 2）您可能没有正确使用边界值。使用您工作的方案，应将值添加到 RHS，但由于它们为零，因此无需添加任何内容。但是，要记住一件事，边界点不应该是计算的一部分。

请参阅以下代码，了解两个问题（两个不同的来源），一个是您使用的简单问题，另一个是我建议在评论中解决的问题。

for Kind=1:2 % define what problem to solve
    errpre=0;%just to initialize this with something

    for N = 1:7

        %create grid
        n=2^(N+1);
        x=linspace(-1,1,n+2);
        y=linspace(-1,1,n+2);
        h= abs(x(1)-x(2));
        [X,Y]=meshgrid(x(2:end-1),y(2:end-1));

        %create descrete operator
        I = speye(n,n);
        E = sparse(2:n,1:n-1,1,n,n);
        D = E+E'-2*I;
        A = (kron(D,I)+kron(I,D))./(h^2);

        %create RHS and the exact solution
        if Kind==1
            f=2* ((X.^2 - 1) + (Y.^2 - 1));
            ex=(X.^2 - 1) .*(Y.^2 - 1);
        elseif Kind==2
            ex = sin(pi*X).*sin(pi*Y);
            f=-ex*2*pi^2;
        end
        %solve the problem
        ap = (A\f(:));

        %compute the error
        err  = norm(ex(:)-ap(:),inf)/norm(ap(:),inf);

        %print the error and the rate
        fprintf('Kind=%d\t n=%d\t err=%-10.8d\t rate=%-4.2f\n',Kind,n,err,log2(errpre/err));

        %keep the old value for the rate calculation at next step
        errpre=err;

    end
end

运行这段代码得到

Kind=1   n=4     err=2.40933816e-16  rate=-Inf
Kind=1   n=8     err=5.69076036e-16  rate=-1.24
Kind=1   n=16    err=8.94357034e-16  rate=-0.65
Kind=1   n=32    err=3.11434148e-15  rate=-1.80
Kind=1   n=64    err=1.11074876e-15  rate=1.49
Kind=1   n=128   err=3.77521199e-15  rate=-1.77
Kind=1   n=256   err=1.42112850e-14  rate=-1.91

Kind=2   n=4     err=1.24859800e-01  rate=-Inf
Kind=2   n=8     err=3.99615153e-02  rate=1.64
Kind=2   n=16    err=1.13319182e-02  rate=1.82
Kind=2   n=32    err=3.01735099e-03  rate=1.91
Kind=2   n=64    err=7.78424525e-04  rate=1.95
Kind=2   n=128   err=1.97680908e-04  rate=1.98
Kind=2   n=256   err=4.98085147e-05  rate=1.99

其它你可能感兴趣的问题

上一篇我怎样才能显示流量是不可压缩的？下一篇求解器线性系统方程