How to subtract two column vectors properly?
2 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
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 is large but should be small considering unew and u are the same. Am I subtracting the columns here incorrectly? Thanks
0 Kommentare
Antworten (1)
Matt J
am 16 Jan. 2024
Bearbeitet: Matt J
am 16 Jan. 2024
The relative error is the only thing you can expect to be small, which it is:
DEsol = norm(unew-u)/norm(u) %relative error
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.
0 Kommentare
Siehe auch
Kategorien
Mehr zu Dialog Boxes finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!