Gauss-Seidal solver for 1D heat equation

2 次查看（过去 30 天）

Dereje 2018-3-22

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/389974-gauss-seidal-solver-for-1d-heat-equation

评论： Dereje 2018-3-23

I tried using Gauss-Seidal solver by computing A and the right hand side. But still got an error:) asking for help.

 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)