数值分析报告编程及运行结果高斯顺序消元法.docx
《数值分析报告编程及运行结果高斯顺序消元法.docx》由会员分享，可在线阅读，更多相关《数值分析报告编程及运行结果高斯顺序消元法.docx（25页珍藏版）》请在冰点文库上搜索。
数值分析报告编程及运行结果高斯顺序消元法
高斯消元法
1.程序:
clear
formatrat
A=input（'输入增广矩阵A='）
[m,n]=size（A）;
fori=1:
（m-1）
numb=int2str（i）;
disp（['第',numb,'次消元后的增广矩阵']）
forj=（i+1）:
m
A（j,:
）=A（j,:
）-A（i,:
）*A（j,i）/A（i,i）;
end
A
end
%回代过程
disp（'回代求解'）
x（m）=A（m,n）/A（m,m）;
fori=（m-1）:
-1:
1
x（i）=（A（i,n）-A（i,i+1:
m）*x（i+1:
m）'）/A（i,i）;
end
x
2.运行结果：
高斯选列主元消元法
1.程序:
clear
formatrat
A=input（'输入增广矩阵A='）
[m,n]=size（A）;
fori=1:
（m-1）
numb=int2str（i）;
disp（['第',numb,'次选列主元后的增广矩阵']）
temp=max（abs（A（i:
m,i）））;
[a,b]=find（abs（A（i:
m,i））==temp）;
tempo=A（a
（1）+i-1,:
）;
A（a
（1）+i-1,:
）=A（i,:
）;
A（i,:
）=tempo
disp（['第',numb,'次消元后的增广矩阵']）
forj=（i+1）:
m
A（j,:
）=A（j,:
）-A（i,:
）*A（j,i）/A（i,i）;
end
A
end
%回代过程
disp（'回代求解'）
x（m）=A（m,n）/A（m,m）;
fori=（m-1）:
-1:
1
x（i）=（A（i,n）-A（i,i+1:
m）*x（i+1:
m）'）/A（i,i）;
end
x
2.运行结果：
追赶法
1.程序：
function[x,L,U]=zhuiganfa（a,b,c,f）
a=input（'输入矩阵-1对角元素a='）;
b=input（'输入矩阵对角元素b='）;
c=input（'输入矩阵+1对角元素c='）;
f=input（'输入增广矩阵最后一列元素f='）;
n=length（b）;
%对A进行分解
u
（1）=b
（1）;
fori=2:
n
if（u（i-1）~=0）
l（i-1）=a（i-1）/u（i-1）;
u（i）=b（i）-l（i-1）*c（i-1）;
else
break;
end
end
L=eye（n）+diag（l,-1）;
U=diag（u）+diag（c,1）;
x=zeros（n,1）;
y=x;
%求解Ly=b
y
（1）=f
（1）;
fori=2:
n
y（i）=f（i）-l（i-1）*y（i-1）;
end
%求解Ux=y
if（u（n）~=0）
x（n）=y（n）/u（n）;
end
fori=n-1:
-1:
1
x（i）=（y（i）-c（i）*x（i+1））/u（i）;
end
2.运行结果：
高斯-塞德尔迭代格式
1.程序：
functionx=Gauss_Seidel（a,b）
a=input（'输入系数矩阵a='）
b=input（'输入增广矩阵最后一列b='）;
e=0.5e-7;
n=length（b）;
N=50;
x=zeros（n,1）;
t=zeros（n,1）;
fork=1:
N
sum=0;
E=0;
t（1:
n）=x（1:
n）;
fori=1:
n
x（i）=（b（i）-a（i,1:
（i-1））*x（1:
（i-1））-a（i,（i+1）:
n）*t（（i+1）:
n））/a（i,i）;
end
ifnorm（x-t）k
break;
end
end
2.运行结果：
雅戈比迭代格式
1.程序：
functionx=Jocabi（a,b）
a=input（'输入系数矩阵a='）;
b=input（'输入增广矩阵最后一列b='）;
e=0.5e-7;
n=length（b）;
N=100;
x=zeros（n,1）;
y=zeros（n,1）;
fork=1:
N
sum=0;
fori=1:
n
y（i）=（b（i）-a（i,1:
n）*x（1:
n）+a（i,i）*x（i））/a（i,i）;
end
fori=1:
n
sum=sum+（y（i）-x（i））^2;
end
ifsqrt（sum）k
break;
else
fori=1:
n
x（i）=y（i）;
end
end
end
ifk==Nwarning（'未能找到近似解'）;
end
2.运行结果：
逐次超松弛法（SOR）
1.程序：
function[n,x]=sor22（A,b,X,nm,w,ww）
%用超松弛迭代法求解方程组Ax=b
%输入：
A为方程组的系数矩阵，b为方程组右端的列向量，X为迭代初值构成的列向量，nm为最大迭代次数，w为误差精度，ww为松弛因子
%输出：
x为求得的方程组的解构成的列向量，n为迭代次数
A=input（'输入系数矩阵A='）;
b=input（'输入方程组右端的列向量b='）;
X=input（'输入迭代初值构成的列向量X='）;
nm=input（'输入最大迭代次数nm='）;
w=input（'输入误差精度w='）;
ww=input（'输入松弛因子ww='）;
n=1;
m=length（A）;
D=diag（diag（A））;%令A=D-L-U,计算矩阵D
L=tril（-A）+D;%令A=D-L-U,计算矩阵L
U=triu（-A）+D;%令A=D-L-U,计算矩阵U
M=inv（D-ww*L）*（（1-ww）*D+ww*U）;%计算迭代矩阵
g=ww*inv（D-ww*L）*b;%计算迭代格式中的常数项
%下面是迭代过程
whilen<=nm
x=M*X+g;%用迭代格式进行迭代
ifnorm（x-X,'inf'）disp（'迭代次数为'）;n
disp（'方程组的解为'）;x
return;
%上面：
达到精度要求就结束程序，输出迭代次数和方程组的解
end
X=x;n=n+1;
end
%下面：
如果达到最大迭代次数仍不收敛，输出警告语句及迭代的最终结果（并不是方程组的解）
disp（'在最大迭代次数内不收敛!
'）;
disp（'最大迭代次数后的结果为'）;
2.运行结果：
二分法求解方程的根
1.程序：
%其中a，b表示查找根存在的范围，M表示要求解函数的值
%f（x）表示要求解根的方程
%eps表示所允许的误差大小
functiony=er_fen_fa（a,b,M）
k=0;
eps=0.05
whileb-a>eps
x=（a+b）/2;
%检查是否大于值
if（（x^3）-3*x-1）>M
b=x
else
a=x
end
k=k+1
end
2.运行结果：
Newton迭代法（切线法）
1.程序：
functionx=nanewton（fname,dfname,x0,e,N）
%newton迭代法解方程组
%fname和dfname分别表示F（x）及其导函数的M函数句柄或内嵌函数，x0为迭代初值，e为精度要求
x=x0;x0=x+2*e;k=0;
ifnargin<5,N=500;end
ifnargin<4e=1e-4;end
whileabs（x0-x）>e&kk=k+1;
x0=x;x=x0-feval（fname,x0）/feval（dfname,x0）;
disp（x）
end
ifk==N,warning（'已达迭代次数上限'）;
end
2.运行结果：
割线方式迭代法
1.程序：
functionx=ge_xian_fa（fname,dfname,x0,x1,e,N）
%割线方式迭代法解方程组
%fname和dfname分别表示F（x）及其导函数的M函数句柄或内嵌函数，x0,x1分别为迭代初值，e为精度要求
k=0;a=x1;b=x0;
ifnargin<5,N=500;end
ifnargin<4e=1e-4;end
whileabs（a-b）>e&kk=k+1;
x=x1-（（x1-x0）/（feval（fname,x1）-feval（fname,x0）））*feval（fname,x1）;
iffeval（fname,x）*feval（fname,x0）>0,
x0=x;b=x0;
else
x1=x;a=x1;
end
x=x1-（（x1-x0）/（feval（fname,x1）-feval（fname,x0）））*feval（fname,x1）;
numb=int2str（k）;
disp（['第',numb,'次计算后x=']）
fprintf（'%f\n\n',x）;
end
ifk==N,warning（'已达迭代次数上限'）;
end
2.运行结果：
Newton插值
1.程序：
%保存文件名为New_Int.m
%Newton基本插值公式
%x为向量，全部的插值节点
%y为向量，差值节点处的函数值
%xi为标量，是自变量
%yi为xi出的函数估计值
functionyi=newton_chazhi（x,y,xi）
n=length（x）;
m=length（y）;
ifn~=m
error（'ThelengthsofXangYmustbeequal!
'）;
return;
end
%计算均差表Y
Y=zeros（n）;
Y（:
1）=y';
fork=1:
n-1
fori=1:
n-k
ifabs（x（i+k）-x（i））error（'theDATAiserror!
'）;
return;
end
Y（i,k+1）=（Y（i+1,k）-Y（i,k））/（x（i+k）-x（i））;
end
end
%计算牛顿插值公式
yi=0;
fori=1:
n
z=1;
fork=1:
i-1
z=z*（xi-x（k））;
end
yi=yi+Y（1,i）*z;
end
2.运行结果：
Lagrange插值
1.程序：
functiony0=Language（x,y,x0）
symstl;
iflength（x）==length（y）
n=length（x）;
else
disp（'x和y的维数不相等！
'）;
return;%检错
end
h=sym（0）;
fori=1:
n
l=sym（y（i））;
forj=1:
i-1
l=l*（t-x（j））/（x（i）-x（j））;
end;
forj=i+1:
n
l=l*（t-x（j））/（x（i）-x（j））;
end;
h=h+l;
end
simplify（h）;
ifnargin==3
y0=subs（h,'t',x0）;%计算插值点的函数值
else
y0=collect（h）;
y0=vpa（y0,6）;%将插值多项式的系数化成6位精度的小数
end
2.运行结果：
最小二乘法
1.程序：
functionp=nafit（x,y,m）
%多项式拟合
%x,y为已知数据点向量,分别表示横,纵坐标,m为拟合多项式的次数,结果返回m次拟合多项式系数,从高次到低次存放在向量p中.
A=zeros（m+1,m+1）;
fori=0:
m
forj=0:
m
A（i+1,j+1）=sum（x.^（i+j））;
end
b（i+1）=sum（x.^i.*y）;
end
a=A\b';
p=fliplr（a'）;
t=0:
0.1:
1.6;
S3=polyval（p,t）;
plot（x,y,'pk'）;
holdon
plot（t,S3,'r'）;
2.运行结果：
升级会员 