数值分析报告文档格式.docx
《数值分析报告文档格式.docx》由会员分享，可在线阅读，更多相关《数值分析报告文档格式.docx（13页珍藏版）》请在冰豆网上搜索。
m=length（y）;
%输入的插值点与它的函数值应有相同的个数
ifn~=m
error（'
ThelengthsofXandYmustbeequal!
'
）;
return;
end
p=zeros（1,n）;
fork=1:
n
t=ones（1,n）;
forj=1:
ifj~=k
%输入的插值点必须互异
ifabs（x（k）-x（j））<
eps
error（'
theDATAiserror!
return;
end
t（j）=（xx-x（j））/（x（k）-x（j））;
p（k）=prod（t）;
yi=sum（y.*p）;
测试数据：
>
x=[01491625364964]
y=[112345678]
xx=8
yi=Lagrange（x,y,xx）
实验结果：
x=01491625364964
y=112345678
xx=8
xx=8
yi=2.9036
截图如下：
实验二：
用不同的方法计算积分
取不同的步长h，分别用复合梯形公式及复合辛普森公式计算积分，给出误差中关于h的函数，并与积分精确值比较两个公式的精度，是否存在一个最小的h，使得精度不能再被改善。
复合梯形公式：
functionI=T_quad（x,y）
n=length（x）;
m=length（y）;
ifn~=m
ThelengthsofXandYmustbeequal'
h=（x（n）-x
（1））/（n-1）;
a=[12*ones（1,n-2）1];
I=h/2*sum（a.*y）;
复合辛普森公式：
functionI=S_quad（x,y）
ifrem（n-1,2）~=0
I=T_quad（x,y）;
N=（n-1）/2;
h=（x（n）-x
（1））/N;
a=zeros（1,n）;
fork=1:
N
a（2*k-1）=a（2*k-1）+1;
a（2*k）=a（2*k）+4;
a（2*k+1）=a（2*k+1）+1;
I=h/6*sum（a.*y）;
复合梯形测试数据：
x=0.00001:
0.0001:
0.99999
y=sqrt（x）.*log（x）
I=T_quad（x,y）
复合辛普森测试数据：
I=S_quad（x,y）
复合梯形实验结果：
x=Columns1through8
0.00000.00010.00020.00030.00040.00050.00060.0007
Columns9through16
…………
Columns9993through10000
0.99920.99930.99940.99950.99960.99970.99980.9999
y=Columns1through8
-0.0364-0.0956-0.1227-0.1422-0.1579-0.1712-0.1828-0.1932
-0.0008-0.0007-0.0006-0.0005-0.0004-0.0003-0.0002-0.0001
I=T_quad（x,y）
I=-0.4444
复合辛普森实验结果：
0.00000.00010.00020.00030.00040.00050.00060.0007
Columns9993through10000
0.99920.99930.99940.99950.99960.99970.99980.9999
I=S_quad（x,y）
实验三：
用LU分解和列主元消去法解线性方程组
输出Ax=b中系数A=LU分解的矩阵L和U，解向量x及detA；
列主元法的行交换次序，解向量x及detA；
比较两种方法所得的结果。
LU分解：
function[L,U,x]=lux（A,b）
[n,n]=size（A）;
p=eye（n）;
fork=1:
n-1
[r,m]=max（abs（A（k:
n,k）））;
m=m+k-1;
if（A（m,k）~=0）
if（m~=k）
A（[km],:
）=A（[mk],:
p（[km]）=p（[mk]）;
fori=k+1:
n
A（i,k）=A（i,k）/A（k,k）;
j=k+1:
n;
A（i,j）=A（i,j）-A（i,k）*A（k,j）;
L=tril（A,-1）+eye（n,n）;
U=triu（A）;
newb=p*b;
y=zeros（n,1）;
n
j=1:
k-1;
y（k）=（newb（k）-L（k,j）*y（j））/L（k,k）;
end
x=zeros（n,1）;
fork=n:
-1:
1
x（k）=（y（k）-U（k,j）*x（j））/U（k,k）;
高斯列主元消去法：
functionx=gauss_lie（A,b）
%UNTITLED5Summaryofthisfunctiongoeshere
%Detailedexplanationgoeshere
%采用高斯列主元法求解方程组Ax=b
n=length（b）;
p=1:
lu=A;
y=[];
[c,i]=max（abs（lu（k:
ik=i+k-1;
ifik~=k
m=p（k）;
p（k）=p（ik）;
p（ik）=m;
ck=lu（k,:
lu（k,:
）=lu（ik,:
lu（ik,:
）=ck;
ifk==n
break;
lu（k+1:
n,k）=lu（k+1:
n,k）/lu（k,k）;
n,k+1:
n）=lu（k+1:
n）-lu（k+1:
n,k）*lu（k,k+1:
n）;
l=diag（ones（n,1））+tril（lu,-1）;
u=triu（lu）;
y
（1）=b（p
（1））;
fori=2:
y（i）=b（p（i））-l（i,1:
i-1）*y（1:
i-1）'
;
x（n）=y（n）/u（n,n）;
fori=n-1:
x（i）=（y（i）-u（i,i+1:
n）*x（i+1:
n）'
）/u（i,i）;
x=x'
LU分解测试数据：
A=[10,-7,0,1;
-3,2.099999,6,2;
5,-1,5,-1;
2,1,0,2]
b=[8;
5.900001;
5;
1]
[l,u,x]=lux（A,b）
det（A）
高斯列主元消去法测试数据：
A=[10,-7,0,1;
b=[8;
x=gauss_lie（A,b）
det（A）
LU分解实验结果：
A=10.0000-7.000001.0000
-3.00002.10006.00002.0000
5.0000-1.00005.0000-1.0000
2.00001.000002.0000
b=
8.0000
5.9000
5.0000
1.0000
l=
1.0000000
0.50001.000000
-0.3000-0.00001.00000
0.20000.9600-0.80001.0000
u=
10.0000-7.000001.0000
02.50005.0000-1.5000
006.00002.3000
0005.0800
x=
0.1949
-0.7661
0.9695
0.6882
ans=
-762.0001
高斯列主元分解实验结果：
A=
8.0000
x=gauss_lie（A,b）
0.0000
-1.0000
升级会员 