# 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)

Show older comments

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

##### 0 Comments

### Answers (1)

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

h = (b-a)/(N-2); % is changing this going to mess up usage for matrix calcs?

x = a:h:b

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
on 9 Oct 2021

### See Also

### Categories

### Community Treasure Hunt

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

Start Hunting!