Gauss-Seidal solver for 1D heat equation
1 次查看(过去 30 天)
显示 更早的评论
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)
回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Matrix Indexing 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!