计算方法上机作业集合Word格式文档下载.docx
《计算方法上机作业集合Word格式文档下载.docx》由会员分享，可在线阅读，更多相关《计算方法上机作业集合Word格式文档下载.docx（32页珍藏版）》请在冰豆网上搜索。
上机作业181页第四题
线性方程组为
（1）顺序消元法
A=[1.1348,3.8326,1.1651,3.4017;
0.5301,1.7875,2.5330,1.5435;
3.4129,4.9317,8.7643,1.3142;
1.2371,4.9998,10.6721,0.0147];
b=[9.5342;
6.3941;
18.4231;
16.9237];
上机代码〔函数部分〕如下
function[b]=gaus（A,b）%用b返回方程组的解
B=[A,b];
n=length（b）;
RA=rank（A）;
RB=rank（B）;
dif=RB-RA;
ifdif>
disp（'
此方程组无解'
）;
return
end
ifRA==RB
ifRA==n
formatlong;
此方程组有唯一解'
forp=1:
n-1
fork=p+1:
n
m=B（k,p）/B（p,p）;
B（k,p:
n+1）=B（k,p:
n+1）-m*B（p,p:
n+1）;
end
end%顺序消元形成上三角矩阵
b=B（1:
n,n+1）;
A=B（1:
n,1:
n）;
b（n）=b（n）/A（n,n）;
forq=n-1:
-1:
1
b（q）=（b（q）-sum（A（q,q+1:
n）*b（q+1:
n）））/A（q,q）;
end%回代求解
else
此方程组有无数组解'
上机运行结果为
A=[1.1348,3.8326,1.1651,3.4017;
X=gaus（A,b）
此方程组有唯一解
X=
1.000000000000000
（2）列主元消元法
函数部分代码如下
function[b]=lieZhu（A,b）%用b返回方程组的解
formatlong;
该方程组无解'
ifdif==0
该方程组有唯一解'
c=zeros（1,n）;
fori=1:
max=abs（A（i,i））;
m=i;
forj=i+1:
ifmax<
abs（A（j,i））
max=abs（A（j,i））;
m=j;
end%求出每一次消元时绝对值最大的一行的行号
ifm~=i
fork=i:
c（k）=A（i,k）;
A（i,k）=A（m,k）;
A（m,k）=c（k）;
d1=b（i）;
b（i）=b（m）;
b（m）=d1;
%函数值跟随方程一起换位置
end
fork=i+1:
A（k,j）=A（k,j）-A（i,j）*A（k,i）/A（i,i）;
b（k）=b（k）-b（i）*A（k,i）/A（i,i）;
A（k,i）=0;
end%完成消元操作，形成上三角矩阵
fori=n-1:
sum=0;
sum=sum+A（i,j）*b（j）;
b（i）=（b（i）-sum）/A（i,i）;
%回代求解其他未知数
else
end上机运行，结果为
X=lieZhu（A,b）
该方程组有唯一解
1.000000000000002
0.999999999999999
根据两种方法运算结果，顺序消元法得到结果比列主元消元法得到的结果精度更高。
〔注：
matlab使用的是2015b版本，不知道是保留小数位数太少，还是程序原因，顺序消元输出结果总是等于准确解，请老师指正〕
第四次上机作业
7.分析用以下迭代法解线性方程组
的收敛性，并求出使
0.0001的近似解及相应的迭代次数。
（1）雅可比迭代法
上机编写的函数如下
function[]=Jacobi（X,b）
%雅可比迭代法
D=diag（diag（X））;
%得到对角线元素组成的矩阵
B=inv（D）*（D-X）;
f=inv（D）*b;
b（:
:
）=0;
x1=B*b+f;
num=1;
while（norm（x1-b）>
0.0001）%判断当前的解是否到达精度要求
b=x1;
x1=B*b+f;
num=num+1;
end;
fprintf（'
求得的解为：
\n'
disp（x1）;
迭代次数：
%d次\n'
num）;
上机运行，结果如下
A=[4,-1,0,-1,0,0;
-1,4,-1,0,-1,0;
0,-1,4,-1,0,-1;
-1,0,-1,4,-1,0;
0,-1,0,-1,4,-1;
0,0,-1,0,-1,4];
b=[0;
5;
-2;
6];
Jacobi（A,b）
0.999981765074381
1.99995018125674
0.999975090628368
1.99996353014876
28次
满足要求的解如输出结果所示，总共迭代了28次
（2）高斯-赛德尔迭代法
上机程序如下所示
function[]=Gauss_Seidel（A,b）
%高斯赛德尔迭代法
D=diag（diag（A））;
L=D-tril（A）;
U=D-triu（A）;
B=inv（D-L）*U;
f=inv（D-L）*b;
x0=B*b+f;
while（norm（x0-b）>
0.0001）
b=x0;
x0=B*b+f;
结果为\n'
disp（x0）;
迭代次数为：
Gauss_Seidel（A,b）
结果为
0.999960143810658
1.99995676152139
0.999963508299833
1.99996662162874
0.999972527179715
1.99998400886989
15次
满足精度要求的解如上述程序打印输出所示，迭代了15次
（3）SOR迭代法〔w依次取1.334,1.95,0.95〕
上机代码如下
function[]=SOR（A,b,w）
%SOR迭代法¨
B=inv（D-w*L）*（（1-w）*D+w*U）;
f=w*inv（D-w*L）*b;
迭代次数为%d\n'
上机运行
SOR（A,b,1.334）
1.00001878481009
1.99998698322858
1.00001815013068
2.00000041318053
0.999991858543476
2.0000068413569
迭代次数为13
SOR（A,b,1.95）
0.999984952088107
2.00000960832604
0.999959115182729
2.0000168426006
1.00000443526697
1.99997885113446
迭代次数为241
SOR（A,b,0.95）
0.999961518309351
1.99995706825231
0.999963054838453
1.99996580572033
0.999971141727589
1.9999827300678
迭代次数为17
由以上输出得到w取值不同的情况下，得到的满足精度要求的结果，迭代次数分别如输出所示
第五次上机作业
8.从函数表
x
0.0
0.1
0.2
0.3
0.401
0.5
f（x）
0.39894
0.39695
0.39142
0.38138
0.36812
0.35206
出发，用以下方法计算f（0.15）,f（0.31）及f（0.47）的近似值
（1）分段线性插值
（2）分段二次插值
（3）全区间上拉格朗日插值
〔1〕线性插值
编写函数如下
function[R]=div_line（x0,y0,x）
%线性插值
p=length（x0）;
n=length（y0）;
m=length（x）;
if（p~=n）%x的个数与y的个数不等
error（'
数据输入有误，请重新输入'
return;
fprintf（'
线性插值\n'
fort=1:
m
z=x（t）;
if（z<
x0
（1）||z>
x0（p））
x[%d]不在所给自变量范围内，无法进行插值'
t）;
continue;
p-1
x0（i+1））
break;
end;
R（t）=y0（i）*（x（t）-x0（i+1））/（x0（i）-x0（i+1））+y0（i+1）*（x（t）-x0（i））/（x0（i+1）-x0（i））;
上机运行如下
x0=[0.00.10.1950.30.4010.5];
y0=[0.398940.396950.391420.381380.368120.35206];
x=[0.150.310.47];
div_line（x0,y0,x）
线性插值
ans=
0.394040.380070.35693
即结果为f（0.15）
0.39404,f（0.31）
0.38007,f（0.47）
0.35693
〔2〕分段二次插值
编写的函数如下
function[R]=div2line（x0,y0,x）
%分段二次插值
m=length（y0）;
n=length（x）;
if（p~=m）
输入错误，请重新输入数据'
fort=1:
if（x（t）<
x0
（1）||x（t）>
x[%d]不在所给区间上'
tag=2;
m=abs（x（t）-x0
（1））+abs（x（t）-x0
（2））+abs（x（t）-x0（3））;
fori=3:
sum=abs（x（t）-x0（i-1））+abs（x（t）-x0（i））+abs（x（t）-x0（i+1））;
if（sum<
m）
m=sum;
tag=i;
tag=%d\n'
tag）;
R（t）=y0（tag-1）*（x（t）-x0（tag））*（x（t）-x0（tag+1））/（（x0（tag-1）-x0（tag））*（x0（tag-1）-x0（tag+1）））+y0（tag）*（x（t）-x0（tag-1））*（x（t）-x0（tag+1））/（（x0（tag）-x0（tag-1））*（x0（tag）-x0（tag+1）））+y0（tag+1）*（x（t）-x0（tag-1））*（x（t）-x0（tag））/（（x0（tag+1）-x0（tag-1））*（x0（tag+1）-x0（tag）））;
End
上机运行，执行结果为
div2line（x0,y0,x）
0.394480.380220.35725
即分段二次插值方法下，
f（0.15）
0.39448,f（0.31）
0.38022,f（0.47）
0.35725
〔3〕
上机编写的程序如下
function[R]=lagrange（x0,y0,x）
%全区间上拉格朗日插值
p=length（y0）;
n=length（x0）;
%计算函数表和x的长度
ifp~=nerror（'
%假设函数表的x与y长度不一致则输入有误
elsefprintf（'
拉格朗日插值\n'
m
%利用循环计算每个x的插值
s=0.0;
fork=1:
n
p=1;
ifi~=k
p=p*（z-x0（i））/（x0（k）-x0（i））;
s=s+y0（k）*p;
%根据拉格朗日插值公式求解y
R（t）=s;
%输出插值结果
lagrange（x0,y0,x）
拉格朗日插值
0.394470.380220.35722
0.39447,f（0.31）
0.35722
9.解：
上机程序如下，为方便起见，将所有操作分在四个函数中进行
入口函数
function[]=spline（X,Y,xx,y1_0,y1_18）
%输出自变量所对应的函数值
M=getM（X,Y,y1_0,y1_18）;
%先得到M
s=xx;
k=length（xx）;
fora=1:
k
s（xx（a））=getExp（X,Y,M,xx（a））;
%计算自变量所在小区间对应曲线的表达式，并根据表达式计算对应的函数值
s（%d）=%f\n'
xx（a）,s（xx（a）））;
%输出打印函数值
获取M
function[M]=getM（X,Y,y1_0,y1_1）
%得到M
n=length（X）;
a=0*X;
b=a;
c=a;
h=a;
f=a;
b=b+2;
h（2:
n）=X（2:
n）-X（1:
n-1）;
%h
（1）不用
a（2:
n-1）=h（2:
n-1）./（h（2:
n-1）+h（3:
n））;
c（2:
n-1）=1-a（2:
a（n）=1;
c
（1）=1;
Yf（2:
n）=Y（2:
n）-Y（1:
f（2:
n-1）=6*（Yf（3:
n）./h（3:
n）-Yf（2:
n-1）./h（2:
n-1））./（h（2:
f
（1）=6*（Yf
（2）/h
（2）-y1_0）/h
（2）;
f（n）=6*（y1_1-Yf（n）/h（n））/h（n）;
M=CalM（n,a,b,c,f）;
%计算M
计算M
function[f]=CalM（n,a,b,c,f）
%追赶法求解M
eps=0.1e-15;
%防止参数过小，是的计算误差过大
ifabs（b
（1））<
eps
除数为0，停止计算'
c
（1）=c
（1）/b
（1）;
%追赶法：
根据递推算式计算β
fori=2:
b（i）=b（i）-a（i）*c（i-1）;
ifabs（b（i））<
c（i）=c（i）/b（i）;
b（n）=b（n）-a（n）*c（n-1）;
%追赶法:
根据递推算式计算
f
（1）=f
（1）/b
（1）;
f（i）=（f（i）-a（i）*f（i-1））/b（i）;
%以下求解Ux=y,x的值存入f
fori=n-1:
f（i）=f（i）-c（i）*f（i+1）;
return
得到自变量所在区间的表达式，并求自变量对应的函数值
function[y]=getExp（X,Y,M,x）
%%根据X、Y、M计算表达式，并根据表达式计算对应的函数值
%%判断x落在哪个小区间
n1=1;
n2=n;
whilen2~=n1+1
n5=fix（（n1+n2）/2）;
ifx>
X（n5）
n1=n5;
n2=n5;
%%
%%计算y
y=M（n1）*（X（n2）-x）^3/（6*h（n2））+M（n2）*（x-X（n1））^3/（6*h（n2））;
y=y+（Y（n1）-M（n1）*h（n2）*h（n2）/6）*（X（n2）-x）/h（n2）;
y=y+（Y（n2）-M（n2）*h（n2）*h（n2）/6）*（x-X（n1））/h（n2）;
%结束
上机试运行，过程如下
X=[0.523.18.017.9528.6539.6250.6578104.6156.6208.6260.7312.5364.4416.3468494507520];
Y=[5.287949.413.8420.224.928.4431.13536.536.634.631.026.3420.914.87.83.71.50.2];
xx=[24612163060110180280400515];
y1_0=1.86548;
y1_18=-0.046115;
spline（X,Y,xx,y1_0,y1_18）
s
（2）=7.825123
s（4）=10.481311
s（6）=12.363477
s（12）=16.575574
s（16）=19.091594
s（30）=25.386402
s（60）=32.804283
s（110）=36.647878
s（180）=35.917148
s（280）=29.368462
s（400）=16.799784
s（515）=0.542713
根据上述程序运行结果，可得所有自变量对应的函数值，如上输出结果所示
第六次上机作业
10.已知一组实验数据
i
2
3
4
5
6
7
8
9
升级会员 