How to subtract two column vectors properly?

Question

Janee 2024-1-16

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2070301-how-to-subtract-two-column-vectors-properly

编辑： Matt J 2024-1-16

cheb.m

I am trying to calculate the discretization error in my code which is simply the difference between my numerical solution and exact solution. This number should be real small (close to zero) which I expect here, but for some reason I am getting a large number.

Code:

N=16;
[D,x] = cheb(N); 
ylow = 0; %a
yupp = 16; %b
Ly = yupp-ylow;
eta_ygl  = 2/Ly;
x = (1/2)*((yupp-ylow)*x + (yupp+ylow));
D=D*eta_ygl; 
D2 = D^2;
D2 = D2(2:N,2:N);
u = (yupp*exp(4*ylow)-ylow*exp(4*yupp)+ylow*exp(4*x)-yupp*exp(4*x)-exp(4*ylow)*x+exp(4*yupp)*x)/(16*ylow-16*yupp);
n = x.^2+10;
dndx = D *n; 
dudx = D *u;
prod1 = dndx .* dudx;
du2dx = D * dudx;
prod2 = n .*du2dx;
invN = (1./n) ;
source = prod1(2:N) + prod2(2:N);
uold = ones(size(u(2:N)));
max_iter =500; 
err_max = 1e-8;
 for iterations = 1:max_iter
    phikmax_old = (max(abs(uold)));
    duoldx = D(2:N,2:N) *uold;
    dnudx = dndx(2:N) .* duoldx;
    ffh = source; 
    RHS = ffh - dnudx; 
    Stilde = invN(2:N) .* RHS;
    unew = D2\Stilde;
    phikmax = (max(abs(unew)));
    if phikmax < err_max 
        it_error = err_max /2;
    else
        it_error = abs( phikmax - phikmax_old) / phikmax;
    end
    if it_error < err_max
        break;
    end
    uold = unew;
 end
unew = [0;unew;0];
DEsol = (unew-u) %HERE
DEsol = 17×1
1.0e+19 *

         0
    0.0077
    0.0340
    0.0876
    0.1811
    0.3252
    0.5203
    0.7432
    0.9367
    1.0166

DEsol is large but should be small considering unew and u are the same. Am I subtracting the columns here incorrectly? Thanks

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Matt J 2024-1-16

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2070301-how-to-subtract-two-column-vectors-properly#answer_1390376

编辑：Matt J 2024-1-16

在 MATLAB Online 中打开

cheb.m

The relative error is the only thing you can expect to be small, which it is:

DEsol = norm(unew-u)/norm(u)   %relative error
DEsol = 2.5021e-08

I don't know if 2.5e-8 is as small as you'd hoped, but it's hard to know how from looking at your code what the sources and magnitudes of the errors are.