数值试验程序.docx
- 文档编号:24707391
- 上传时间:2023-05-31
- 格式:DOCX
- 页数:75
- 大小:1.21MB
数值试验程序.docx
《数值试验程序.docx》由会员分享,可在线阅读,更多相关《数值试验程序.docx(75页珍藏版)》请在冰豆网上搜索。
数值试验程序
数值计算方法
数值试验
目录
1.数值试验程序—第二章第1题
2.数值试验程序—第二章第2题
3.数值试验程序—第二章第4题
4.数值试验程序—第三章第1题
5.数值试验程序—第四章第4题
6.数值试验程序—第五章第1题
7.数值试验程序—第五章第2题
8.数值试验程序—第五章第3题
9.数值试验程序—第五章第4题
10.数值试验程序—第六章第1题
11.数值试验程序—第六章第2题
12.数值试验程序—第六章第3题
13.数值试验程序—第七章第1题
14.数值试验程序—第七章第2题
15.数值试验程序—第七章第3题
16.数值试验程序—第八章第1题
17.数值试验程序—第八章第3题
18.数值试验程序—第八章第4题
1.数值试验程序—第二章第1题
(1)程序源代码
functionXA=liezhuyuan2_1(A)
%第二章数值试验第1题
%=====================================================================
%列主元分解法
%=====================================================================
%A-------输入矩阵
%XA-------消元后的矩阵
%=====================================================================
N=size(A);
ifN
(1)~=N
(2)
disp('error:
输入矩阵不是方阵!
');
XA=-1;
return;
end
n=N
(1);
index=0;%标记
fori=1:
(n-1)
me=max(abs(A(1:
n,i)));%第i列的最大值
fork=i:
n
if(abs(A(k,i))==me)
index=k;%找出第i列最大值所在的行k,标记
break;
end
end
temp=A(i,1:
n);%交换第i行和第k行
A(i,1:
n)=A(index,1:
n);
A(index,1:
n)=temp;
forj=(i+1):
n%判断对角元素是否为0
if(A(i,i)==0)
disp('对角元素为0!
')
return;
end
l=A(j,i);%消元
m=A(i,i);
A(j,1:
n)=A(j,1:
n)-l*A(i,1:
n)/m;
end
end
XA=A;
%=====================================================================
(2)matlab运算结果
在matlab中输入以下代码:
>>A=[11111;12345;1361015;14102035;15153570];
X=liezhuyuan2_1(A)
得到A的列主元消去法后的矩阵:
2.数值试验程序—第二章第2题
(1)程序源代码
追赶法分解对三角矩阵
functionx=sjzhuigan2_2(a,b,c,d)
%第二章数值试验第2题
%=====================================================================
%追赶法解三对角方程组
%a---------下对角线元素,行向量
%b---------对角线元素,行向量
%c---------上对角线元素,行向量
%d---------方程组等号右边常数列向量
%x---------方程组的解,列向量
%=====================================================================
if(length(a)+1~=length(b))||(length(c)+1~=length(b))
disp('errorin==>输入向量长度有误!
');
x=-1;
return;
end
n=length(b);
u
(1)=b
(1);y
(1)=d
(1);
fori=1:
n-1
l(i)=a(i)/b(i);%求得l(i)
u(i+1)=b(i+1)-c(i)*l(i);%求得u(i)
y(i+1)=d(i+1)-l(i)*y(i);%由Ly=d求得y
end
x(n)=y(n)/u(n);
fori=n-1:
-1:
1
x(i)=(y(i)-c(i)*x(i+1))/u(i);%由Ux=d求得x
end
x=x’;
%=====================================================================
(2)matlab运行结果
输入以下语句:
>>a=[-2-2-2-2-2-2-2];
b=[25555555];
c=a;
d=[220/27;0;0;0;0;0;0;0];
x=sjzhuigan2_2(a,b,c,d)
得到结果为:
3.数值试验程序—第二章第4题
(1)程序源代码
%第二章数值试验第4题
%=====================================================================
%第1问,求det(A),cond(A),A的所有特征值
%=====================================================================
clc;
A=[10787;
7565;
86109;
75910];
b=[32;23;33;31];
det_A=det(A);%矩阵A的行列式值
cond_A2=cond(A,2);%矩阵A的2条件数
cond_Ainf=cond(A,inf);%矩阵A的无穷条件数
lambda=eig(A);%矩阵A的特征值
fprintf('\n计算结果\n\n第1问:
\n\n');
fprintf('矩阵A的行列式值:
|A|=%f\n\n',det(A));
fprintf('矩阵A的2条件数:
cond(A,2)=%f\n\n',cond_A2);
fprintf('矩阵A的无穷条件数:
cond(A,inf)=%f\n\n',cond_Ainf);
fprintf('矩阵A的特征值:
λ1=%f,λ2=%f,λ3=%f,λ4=%f\n\n',lambda);
%=====================================================================
%第2问,求δx和δx的2范数,分析方程组解的误差与A的误差的关系
%=====================================================================
x=A\b;
delta_A=[00.20.1-0.1;
0.080.070.020;
0.2-0.11-0.040.01;
-0.020.04-0.03-0.02];%向量δA
delta_x=(A+delta_A)\b-x;%向量δx的值
norm_delta_x2=norm(delta_x,2);%向量δx的2范数
norm_x2=norm(x,2);%向量x的2范数
norm_A2=norm(A,2);%矩阵A的2范数
norm_delta_A2=norm(delta_A,2);%矩阵δA的2范数
wucha_A=norm_delta_A2/norm_A2;%矩阵A的相对误差
wucha_solution=norm_delta_x2/norm_x2;%方程组解的相对误差
fprintf('第2问:
\n\n');
fprintf('向量δx:
δx=[%f;%f;%f;%f]\n\n',delta_x);
fprintf('向量δx的2范数:
‖δx‖2=%f\n\n',norm_delta_x2);
fprintf('矩阵A的相对误差:
‖δA‖2/‖A‖2=%f\n\n',wucha_A);
fprintf('方程组解的相对误差:
‖δx‖2/‖x‖2=%f\n\n',wucha_solution);
%=====================================================================
%第3问,在新的随机扰动矩阵δAA下重复第2问
%=====================================================================
delta_AA=0.5*0.0001*rand(4);%新的扰动矩阵δAA
delta_xx=(A+delta_AA)\b-x;%新的扰动下向量δxx的值
norm_delta_xx2=norm(delta_xx,2);%向量δxx的2范数
norm_delta_AA2=norm(delta_AA,2);%矩阵δAA的2范数
wucha_AA=norm_delta_AA2/norm_A2;%新扰动下矩阵A的相对误差
wucha_solution1=norm_delta_xx2/norm_x2;%新扰动下方程组解的相对误差
fprintf('第3问:
\n\n');
fprintf('新的随机扰动矩阵δAA:
');
delta_AA
fprintf('新的扰动下的向量δxx:
δxx=[%e;%e;%e;%e]\n\n',delta_xx);
fprintf('新的扰动下向量δxx的2范数:
‖δxx‖2=%e\n\n',norm_delta_xx2);
fprintf('新的扰动下矩阵A的相对误差:
‖δAA‖2/‖A‖2=%e\n\n',wucha_AA);
fprintf('新的扰动下方程组解的相对误差:
‖δxx‖2/‖x‖2=%e\n\n',wucha_solution1);
fprintf('结论\n比较第2,3问的结果发现,A的相对误差越小,方程组解的相对误差越小.\n');
%=====================================================================
(2)matlab运算结果
(3)分析
从理论结果式来分析
由书本46页内容可知道,当δA充分小,使得‖A-1‖*‖δA‖<1时,可得到公式(2-39),即:
(2-39)
已经计算出:
矩阵A的2条件数:
cond(A,2)=2984.092702(很大,说明A很病态)
矩阵A的相对误差:
‖δA‖2/‖A‖2=0.009683,
代入2-39式可得:
;
得出一个负的结果,这明显不对。
造成这个结果的原因是‖A-1‖*‖δA‖>1,不满足前提条件,造成了分母为负值。
实际上,‖A-1‖2*‖δA‖2=28.8964,远远大于1.
从实际计算方面分析
(a).原扰动下矩阵A的相对误差:
‖δA‖2/‖A‖2=0.009683;
方程组解的相对误差:
‖δx‖2/‖x‖2=0.822463;
(b).新的扰动下矩阵A的相对误差:
‖δAA‖2/‖A‖2=4.149810e-006;
新的扰动下方程组解的相对误差:
‖δxx‖2/‖x‖2=2.091492e-005;
比较结果可以得出结论:
不同的扰动δA下,矩阵A的相对误差越小,方程组解的相对误就差越小.
实际上,第三问采用了能够产生均匀分布的随机矩阵函数rand(n)来产生新的扰动矩阵,并且使其元素的绝对值不超过0.00005。
每一次运行程序时,由rand函数产生的扰动矩阵都是不一样的,经过多次运行发现,计算结果都是满足上述结论的。
4.数值试验程序—第三章第1题
(1)程序源代码
a.Jacobi迭代法自定义函数
function[x,n]=jacobi3_1(A,b,x0,ep,N)
%第三章数值试验第1题
%=====================================================================
%jacobi法解线性方程组
%=====================================================================
%A------------系数矩阵
%b------------常数向量
%x0------------迭代初始向量
%ep------------精度控制(可选择)
%N------------迭代步数控制(可选择)
%x------------线性方程组的解
%n------------满足精度要求的实际迭代步数
%=====================================================================
ifnargin==3%判断输入参数的个数
ep=1e-6;
N=10000;
elseifnargin==4
N=10000;
end
D=diag(diag(A));%求A的对角矩阵
L=-tril(A,-1);%求A的下三角阵
U=-triu(A,1);%求A的上三角阵
M=D\(L+U);
g=D\b;
x=x0;
n=0;%迭代次数
tol=1;%前后两次迭代的差值
whiletol>=ep
x=M*x0+g;%迭代公式
n=n+1;
tol=norm(x-x0);
x0=x;
if(n>=N)
disp('Warning:
迭代次数太多,可能不收敛!
')
return;
end
end
%=====================================================================
b.Guass-seidel迭代法自定义函数
function[x,n]=guass_seidel3_1(A,b,x0,ep,N)
%第三章数值试验第1题
%=====================================================================
%Guass-seidel法解线性方程组
%=====================================================================
%A------------系数矩阵
%b------------常数向量
%x0------------迭代初始向量
%ep------------精度控制(可选择)
%N------------迭代步数控制(可选择)
%x------------线性方程组的解
%n------------满足精度要求的实际迭代步数
%=====================================================================
ifnargin==3%判断输入参数的个数
ep=1e-6;
N=10000;
elseifnargin==4
N=10000;
end
D=diag(diag(A));%求A的对角矩阵
L=-tril(A,-1);%求A的下三角阵
U=-triu(A,1);%求A的上三角阵
M=(D-L)\U;%迭代矩阵M
g=(D-L)\b;
x=x0;
n=0;%迭代次数
tol=1;%前后两次迭代的差值
whiletol>=ep
x=M*x0+g;%迭代公式
n=n+1;
tol=norm(x-x0);
x0=x;
if(n>=N)
disp('Warning:
迭代次数太多,可能不收敛!
')
return;
end
end
%=====================================================================
c.共轭梯度迭代法自定义函数
function[x,n]=gongetidu3_1(A,b,x0,ep)
%第三章数值试验第1题
%=====================================================================
%共轭梯度法解线性方程组
%=====================================================================
%A------------系数矩阵
%b------------常数向量
%x0------------迭代初始向量
%ep------------精度要求
%x------------线性方程组的解
%n------------满足精度要求的实际迭代步数
%=====================================================================
d0=b-A*x0;%求初值d0,r0,lambda0,x1
r0=d0;
lambda0=dot(r0,d0)/dot(d0,A*d0);
x1=x0+lambda0*d0;
n=1;%迭代次数
r=r0;
whilenorm(r,2)>ep%停机条件
r=b-A*x1;
belta=-dot(r,A*d0)/dot(d0,A*d0);
d=r+belta*d0;
lambda=dot(r,d)/dot(d,A*d);
x=x1+lambda*d;
d0=d;
x1=x;
n=n+1;
end
%=====================================================================
(2).matlab运行结果
输入以下代码:
>>A=[101234;19-12-3;2-173-5;32312-1;4-3-5-115];
b=[12;-27;14;-17;12]
x0=[0;0;0;0;0];
ep=1e-6;
[x1,n1]=jacobi3_1(A,b,x0,ep);
[x2,n2]=guass_seidel3_1(A,b,x0,ep);
[x3,n3]=gongetidu3_1(A,b,x0,ep);
X=[x1,x2,x3]
N=[n1,n2,n3]
运行结果为:
其中,矩阵X第一列,第二列和第三列分别为Jacobi法,Guass_seidel法和共轭梯度法得到的线性方程组的解;向量N中的元素分别为三种方法的迭代次数。
可见,共轭梯度法的迭代次数很少,收敛非常快。
5.数值试验程序—第四章第4题
(1)程序源代码
%第四章第四题
%==================================================================================
%第一问
%==================================================================================
clc;
A=[19066-8430;
6630342-36;
336-168147-112;
30-3628291];
N=100;%迭代次数
ep=1e-6;%精度要求
x=[0;0;0;1];%初始向量x
k=0;%从k=0开始迭代
mu=0;%μ置初值0
whilek<=N%当k小于迭代次数N时
alpha=max(abs(x));%α置向量x(k)的无穷范数
y=x/alpha;%y=x(k)/α
x=A*y;%x(k+1)=A*y
lambda=alpha;%最大特征值λ=α
ifabs(lambda-mu) break; end k=k+1;mu=lambda;%反之,k加1,λ的值置给μ,进行下一次循环 end x=x/norm(x); fprintf('\n------------------第一问结果------------------------------\n') fprintf('\n矩阵A的按模最大的特征值λ=%f\n\n',lambda); fprintf('λ对应的特征向量=[%f;%f;%f;%f]\n',x); %================================================================================== %第二问 %================================================================================== v=x;%按模最小的特征向量v m=0; l=0;%按模最小的特征值 for(k=1: N) yy=A\v; m=max(abs(yy)); v=yy/m; if(abs(m-l) l=l/m; return; end end fprintf('\n-----------------------第二问结果-------------------------\n') fprintf('\n矩阵A的按模最小的特征值λ=%f\n\n',l); fprintf('λ对应的特征向量=[%f;%f;%f;%f]\n',v); %=======================================================================
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值 试验 程序
![提示](https://static.bdocx.com/images/bang_tan.gif)