1-D Heat equation

I have no idea what went wrong. Two days since looking for the errors, please help

 if
close all
x_min= 0;
x_max = 1;
N = 10; 
L = x_max-x_min; 
h = L/N; 
x = linspace(x_min,x_max,N+1)'
uexact= @(x) 4*(x-x.^2);
f = @(x) 8*15; % Source term 
u_exact = uexact(x); 
ua = uexact(x_min); 
ub = uexact(x_max);
u = zeros(N-1,1); 
A = zeros(N-1,N-1); 
b = f(x(2:N)); % 
b(1) = b(1) + ua/h^2; 
b(N-1) = b(N-1) + ub/h^2;
dA = diag( 2*ones(1,N-1) ); 
dAp1 = diag( -1*ones(1,N-2), 1 ); 
dAm1 = diag( -1*ones(1,N-2), -1 ); 
A = (dA + dAp1 + dAm1);
A = 15*A/h^2;
u = A\b; % solving the linear system
g = [ua; u; ub]; 
plot(x,u_exact,'b',x,g,'ro-');
legend('exact','numerical');

Thanks for the help.

Answer 1

[Birdman]

x_min= 0;
x_max = 1;
N = 10; 
L = x_max-x_min; 
h = L/N; 
x = linspace(x_min,x_max,N+1)'
uexact= @(x) 4*(x-x.^2);
f = @(x) 8*15*ones(1,N-1); % Source term 
u_exact = uexact(x); 
ua = uexact(x_min); 
ub = uexact(x_max);
u = zeros(N-1,1); 
A = zeros(N-1,N-1); 
b = f(x(2:N)); %
b(1) = b(1) + ua/h^2; 
for i=2:N-1
    b(i) = b(i) + ub/h^2;
end
dA = diag( 2*ones(1,N-1) ); 
dAp1 = diag( -1*ones(1,N-2), 1 ); 
dAm1 = diag( -1*ones(1,N-2), -1 ); 
A = (dA + dAp1 + dAm1);
A = 15*A/h^2;
u = A\b.'; % solving the linear system
g = [ua; u; ub]; 
plot(x,u_exact,'b',x,g,'ro-');
legend('exact','numerical');

[Dereje]

You did it with 5 min which took me days :) Thanks a lot!!

[Birdman]

You are welcome :)

[Dereje]

When I entered new

   uexact= @(x) 8*(x-3*x.^2+4*x.^3-2*x.^4)
   f = @(x) 48*15*(2*x-1).^2

The previous error appeared.

[Birdman]

I assume you got the error on the line

 u=A\b.';

Change that line to

 u=A\b;

Be careful to do that operation according to whether b is a row or column vector.

[Dereje]

Yes I see that. Thanks again!!

[Dereje]

Now I wanted to do the Gauuss-Sidel solver by taking the above A and the right side as a n input and tried to run the code below. But still got an error:) and asking for your help again.

if 
 x_min= 0;
x_max = 1;
N = 10; 
lamda=15;
L = x_max-x_min; 
h = L/N; 
x = linspace(x_min,x_max,N+1)';
uexact= @(x) 8*(x-3*x.^2+4*x.^3-2*x.^4);
% uexact= @(x) 4*(x-x.^2);
% f = @(x) 8*lamda*ones(1,N-1); % Source term 
f = @(x) 48*lamda*(2*x-1).^2;
u_exact = uexact(x); 
ua = uexact(x_min); 
ub = uexact(x_max);
u = zeros(N-1,1); 
A = zeros(N-1,N-1); 
b = f(x(2:N)); %
b(1) = b(1) + ua/h^2; 
for i=2:N-1
  b(i) = b(i) + ub/h^2;
end
A1 = diag( 2*ones(1,N-1) );
A2 = diag( -1*ones(1,N-2), 1 ); 
A3 = diag( -1*ones(1,N-2), -1 ); 
A = (A1 + A2 + A3);
A = lamda*A/h^2;
x0=zeros(N-1,1);
tol=0.0001;
M=1000000;
k=1;
x2=x0;
while(k<M)
  n = length(b);
  x2(1)=(b(1)-A(1,2)*x0(2)-A(1,3)*x0(3))/A(1,1);
for j = 2 : n
   x2(j)=(b(j)-A(j,1:j-1).*x2(1:j-1)-A(j,j+1:n).*x0(j+1:n))./A(j,j);
end
  x1=x2';
  error=abs(x1(j))-x0;
  if norm(error)< tol%if norm(x1-x0)< tol
       disp('succefful');

      break  
    end
end    
plot(x,x1)

[Birdman]

This looks like a new question. Ask it to the forum as a new question so more people can contribute.

[Dereje]

Ok, I will do that.