24小时热门版块排行榜

返回列表

duinichixin

金虫 (著名写手)

应助: 1 (幼儿园)
金币: 892
散金: 596
红花: 4
帖子: 1009
在线: 338.5小时
虫号: 2475163
注册: 2013-05-18
性别: GG
专业: 计算数学与科学工程计算

[求助] 小弟求教几个程序-关于迭代法的（多谢）已有1人参与

小弟编写了Table4.1SOR中的r=0.01及n=8的程序，现有如下几个疑问：1、我运行出it为8，cpu=0.042，w(exp)=1.02，为什么与作者不一样？2、作者终止条件是F范数，我的是2范数，应该怎么修改？3.当n=64时出现‘Out of memory. Type HELP MEMORY for your options.’请问是我程序太复杂了吗？
注：test为确定w(exp)程序，SOR为调用程序，test2为根据w(exp)计算it程序，各位大神不吝赐教。
test.m:
n=8;
r=0.01;
a=ones(n,1);
b=ones(n-1,1);
I=eye(n);
MM=2*diag(a)-diag(b,-1)-diag(b,1);
N=0.5*diag(b,-1)-0.5*diag(b,1);
A=MM+2*r*N+100*I/(n+1);
B=A;
X=rand(n,n);
F=A*X+X*B;
E=kron(I,A)+kron(B',I);
b=F(

;
%x=E\b;
x0=b;
ww=0.01:0.01:1.99;
D=diag(diag(E)); %求E的对角矩阵
L=-tril(E,-1);    %求E的下三角矩阵
U=-triu(E,1);    %求E的上三角矩阵
for i=1:length(ww)
w=ww(i);
B=inv(D-L*w)*((1-w)*D+w*U);
p(i)=max(max(abs(eig(B))));
%plot(ww,p);
[x,n1]=SOR(E,b,x0,w);
mm(i)=n1;
end
plot(ww,mm);
Xe=reshape(x,n,n);
errror=norm((Xe-X),2)

SOR.m:
function [x,n1]=SOR(E,b,x0,w,eps,M)

if nargin==4
eps=1.0e-6;
M=2000;
elseif nargin<4
error
return
elseif nargin==5
M=200;
end
if(w<=0 || w>=2) %收敛条件要求
error;
return;
end
D=diag(diag(E)); %求E的对角矩阵
L=-tril(E,-1);    %求E的下三角矩阵
U=-triu(E,1);    %求E的上三角矩阵
B=inv(D-L*w)*((1-w)*D+w*U);
f=w*inv((D-L*w))*b;
x=B*x0+f;
n1=1;             %迭代次数
%迭代过程
while norm(x-x0)>=eps
x0=x;
x=B*x0+f;
n1=n1+1;
if(n1>=M)
      disp('Warining:迭代次数太多，可能不收敛！');
      return;
end
end

test2.m
w=1.02
n=8;
r=0.01;
a=ones(n,1);
b=ones(n-1,1);
I=eye(n);
M=2*diag(a)-diag(b,-1)-diag(b,1);
N=0.5*diag(b,-1)-0.5*diag(b,1);
A=M+2*r*N+100*I/(n+1);
B=A;
X=rand(n,n);
F=A*X+X*B;
E=kron(I,A)+kron(B',I);
b=F(

;
%x=E\b;
x0=b;
[x,n1]=SOR(E,b,x0,w)
Xe=reshape(x,n,n);
errror=norm((Xe-X),2)