该楼层疑似违规已被系统折叠 隐藏此楼查看此楼

%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);