该楼层疑似违规已被系统折叠 隐藏此楼查看此楼
%SOR迭代法
function [y,k,err,w]=sor(eps,a,n)
h=1/n;
A=zeros(n-1,n-1); %定义系数矩阵A
for i=1:n-1
A(i,i)=-(2*eps+h);
end
for i=1:n-1
for j=1:n-1
if i==j+1
A(i,j)=eps;
end
if i==j-1
A(i,j)=eps+h;
end
end
end
b=zeros(n-1,1); %定义常数项b
for i=1:n-2
b(i,1)=a*h^2;
end
b(n-1,1)=a*h^2-eps-h; %考虑边值问题
D=diag(diag(A));
L=triu(A)-A;
U=tril(A)-A;
B=D\(L+U);
x=eig(B); %B的特征值
p=abs(max(x)); %B谱半径
w=2/(1+sqrt(1-p^2)); %最佳松弛因子
if w<0||w>2
disp('迭代不收敛!');
return;
end;
L=(D-w*L)\((1-w)*D+w*U);
f=w*(D-w*L)\b;
q=1-h*h;
delta=(1.0e-4)*(1-q)/q;
y=zeros(n-1,1);
z=zeros(n-1,1);
k=0;
while 1
z=L*y+f;
if norm(z-y,inf)
break;
end
y=z;
k=k+1;
end
x=[h:h:(n-1)*h];
true=(1-a)/(1-exp(-1/eps))*(1-exp(-x/eps))+x*a;
t=true';
err=norm(t-y,inf);