I have a problem with dimensions of vektor and matrix, i do not know why my f vektor is 11x1, it should be 9x1 to be compatible with my matrix A which is 9x9.

3 views (last 30 days)
Ammar Anayi
Ammar Anayi on 9 Oct 2021
Commented: DGM on 9 Oct 2021
u_exact = @(x) ( 300 * exp(-(x-1).^2) ); % exakt lösning
L = 2;
a = 0; b = L;
ga = 300; gb = 400; % RV
px = @(x) ( 5 + ((1/6) * x.^2) ); %p(x)
p_primx = @(x) ((2/6)*x); %p'(x)
N = 10; % antal delintervall
h = (b-a)/N; % nätsorlek
x = a:h:b; % punkter x_j
p = px(x); % p(x_j)
p_prim = p_primx(x); % p'(x_j)
% Matrisen A
A = zeros( N-1 , N-1 );
for j = 2:(N-2)
A(j,j+1)= -p(j)/h^2-p_prim(j)/(2*h);
A(j,j-1)= -p(j)/h^2+p_prim(j)/(2*h);
A(j,j) = 2*p(j)/h^2 + 1; %q_j = 1.
end
A(1,1) = 1;
A(N-1,N-1)=1;
% A(1,1) = 2*px(x(2))/h^2;
% A(1,2) = -p_primx((x(2)))/h^2;
% A(N-1,N-2) = p_primx((x(N)))/(2 * h) - px(x(2))/h^2;
% A(N-1,N-1) = 2*px(x(N))/h^2;
Q_start = 300*exp(-(((a+h)-(L/2))^2));
Q_slut = 300*exp(-(((b-h)-(L/2))^2));
f = [( 300 * exp(-(x-1).^2) ) ]';
% f(1) = ga;
% f(N-1)=gb;
f(1) = Q_start - (((((p_primx(a+h))/(2*h))-((px(a+h))/(h^2)))*ga));
f(N-1) = Q_slut - ((((-(p_primx(b-h))/(2*h))-((px(b-h))/(h^2)))*gb));
u = A\f;
plot( x , u )
% felet beräknat med trapetsregeln
error = 0;
for j = 1:N
error = error + h/2 * ( ( u_exact(x(j))-u(j) )^2 + ( u_exact(x(j+1))-u(j+1) )^2 );
end
error = sqrt( error )
hold on
fplot( u_exact, [a,b] )
hold off

Answers (1)

DGM
DGM on 9 Oct 2021
Edited: DGM on 9 Oct 2021
The problem is in how you construct x. This looks to me like a classic fencepost error problem, but it's compounded by the fact that you set N to the matrix size +1 and are using h for both the matrix calculations and the vector construction.
a = 0; b = 2;
N = 10;
h = (b-a)/N;
x = a:h:b
x = 1×11
0 0.2000 0.4000 0.6000 0.8000 1.0000 1.2000 1.4000 1.6000 1.8000 2.0000
h = (b-a)/(N-2); % is changing this going to mess up usage for matrix calcs?
x = a:h:b
x = 1×9
0 0.2500 0.5000 0.7500 1.0000 1.2500 1.5000 1.7500 2.0000
I'm confused as to why you set N to 10 and then always use N-1. Why not set N to 9 and just use N?
On that note, you'll have to change the terminal index for the error vector loop to N-2.
  2 Comments
DGM
DGM on 9 Oct 2021
Well if it's dictated that you use N-1, then N would indeed need to be 10, I guess. You'll have to work around that accordingly in order for x to be the right length.

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!