数值分析算法编程..doc
文本预览下载声明
2013年工程数值分析作业
Matlab编程
一编程求解
试分别选用Jacobi 迭代法,Gauss-Seidol迭代法计算其解。
Jacobi 迭代法
function[x,n]=jacobi(A,b,x0,eps,varargin)
%采用Jacobi迭代法求线性方程组Ax=b的解
%线性方程组的系数矩阵: A
%线性方程组的常数向量: b
%迭代初始向量: x0
%解的精度控制: eps
%迭代步数控制: varargin
%线性方程组的解: x
%求出所需精度的解实际的迭代步数: n
if nargin==3
eps=1.0e-6;
M=200;
elseif nargin3
return
elseif nargin==5
M=varargin{1};
end
D=diag(diag(A)); %求A的对角矩阵
L=-tril(A,-1); %求A的下三角阵
U=-triu(A,1); %求A的上三角阵
B=D\(L+U);
f=D\b;
x=B*x0+f;
n=1; %迭代次数
%迭代过程
while norm(x-x0)=eps
x0=x;
x=B*x0+f;
n=n+1;
if(n=M)
disp(‘warning:迭代次数太多,可能不收敛!’);
return;
end
end
在MATLAB命令窗口中输入求解程序:
A=[7 2 1 -2;9 15 3 -2;-2 -2 11 5;1 3 2 13];
b=[4 7 -1 0];
x0=ones(4,1);
[x,n]=jacobi(A,b,x0)
得
x =
0.4981
0.1443
0.0629
-0.0812
n =
8
Gauss-Seidol迭代法
clear;
clc;
A=[7 2 1 -2;9 15 3 -2;-2 -2 11 5;1 3 2 13];
b=[4 7 -1 0];
N=length(b); %解向量的维数
fprintf(库函数计算结果:);
x=inv(A)*b; %库函数计算结果
x=zeros(N,1);%迭代初始值
%-----(A=D-E-F)------
D=diag(diag(A));
E=-tril(A,-1);%下三角
F=-triu(A,1);%上三角
B=inv(D-E)*F;
g=inv(D-E)*b;
eps=0.0001;%相邻解的距离小于该数时,结束迭代
%--------开始迭代-------
for k=1:100 %最大迭代次数为100
fprintf(第%d次迭代:,k);
y=B*x+g;
fprintf(\n与上次计算结果的距离(2范数):%f \n,norm(x-y)^2);
if norm(x-y)eps
break;
end
x=y
end
x
结果
第1次迭代:
与上次计算结果的距离(2范数):0.349202
x =
0.571428571428571
0.123809523809524
0.035497835497836
-0.077988677988678
第2次迭代:
与上次计算结果的距离(2范数):0.005125
x =
0.508700822986537
0.143948115376687
0.063203751515440
-0.082073282472883
第3次迭代:
与上次计算结果的距离(2范数):0.000119
x =
0.497821921826489
0.144389658937968
0.063162688535757
-0.081332021208608
第4次迭代:
与上次计算结果的距离(2范数):0.000000
x =
0.497913421595870
0.144441806507512
0.062851869295437
-0.081303275362252
第5次迭代:
与上次计算结果的距离(2范数):0.000000
x =
0.497913421595870
0.144441806507512
0.062851869295437
-0.081303275362252
对于精度,从程序计算结果来看,收敛最快,最慢。
二 欧拉法 MATLAB
取步长,用Eu
显示全部