数值分析上机作业修改.docx
《数值分析上机作业修改.docx》由会员分享，可在线阅读，更多相关《数值分析上机作业修改.docx（15页珍藏版）》请在冰豆网上搜索。
数值分析上机作业修改
西南交通大学
数值分析2015上机实习报告
2015年11月
目录
第2题1
1.程序1
2.结果分析2
第3题5
1.程序6
2.结果分析7
第5题8
1.程序8
2.结果分析9
第2题
2.某过程测涉及两变量x和y,拟分别用插值多项式和多项式拟合给出其对应规律的近似多项式，已知xi与yi之间的对应数据如下，xi=1,2,…,10
yi=34.658840.371914.6448-14.2721-13.357024.823475.2795
103.574397.484778.2392
（1）请用次数分别为3，4，5，6的多项式拟合并给出最好近似结果f（x）。
（2）请用插值多项式给出最好近似结果
下列数据为另外的对照记录，它们可以作为近似函数的评价参考数据。
xi=
Columns1through7
1.50001.90002.30002.70003.10003.50003.9000
Columns8through14
4.30004.70005.10005.50005.90006.30006.7000
Columns15through17
7.10007.50007.9000
yi=
Columns1through7
42.149841.462035.118224.385211.2732-1.7813-12.3006
Columns8through14
-18.1566-17.9069-11.02262.028419.854940.362661.0840
Columns15through17
79.568893.7700102.3677
1.程序
（1）多项式拟合程序
n=input（'输入所要拟合的阶数n='）;
holdon;
x=1:
10;
y=[34.658840.371914.6448-14.2721-13.357024.823475.2795103.574397.484778.2392];
P=polyfit（x,y,n）
xi=1:
.2:
10;
yi=polyval（P,xi）;
plot（xi,yi,x,y,'r*'）;
legend（'拟合曲线','原始数据'）
（2）拉格朗日插值多项式拟合程序
clc;
x=1:
10;
y=[34.658840.371914.6448-14.2721-13.357024.823475.2795103.574397.484778.2392];
plot（x,y,'r*'）;
holdon;
symst;
n=length（x）;
f=0.0;
fori=1:
n
l=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;
f=f+l;
simplify（f）;
f=collect（f）;
f=vpa（f,6）;
end
ti=1.0:
0.4:
10;
f=subs（f,'t',ti）;
plot（ti,f）
legend（'拟合曲线','原始数据'）
title'拉格朗日插值多项式拟合'
2.结果分析
（1）请用次数分别为3，4，5，6的多项式拟合并给出最好近似结果f（x）。
图2-13次多项式拟合结果
图2-24次多项式拟合结果
图2-35次多项式拟合结果
图2-46次多项式拟合结果
从绘制的图形来看，当采用6次多项式拟合的时候，拟合的曲线已经与所给出的点非常逼近了。
6次多项式拟合曲线为：
f（x）=0.0194x6-0.5408x5+5.1137x4-16.8973x3-0.8670x2+66.3750x-18.6991
（2）拉格朗日插值多项式给出最好近似结果
图1-4拉格朗日插值多项式拟合
第3题
3．用雅格比法与高斯－赛德尔迭代法解下列方程组Ax=b，研究其收敛性，上机验证理论分析是否正确，比较它们的收敛速度，观察右端项对迭代收敛有无影响。
（1）A行分别为A1=[6,2,-1],A2=[1,4,-2],A3=[-3,1,4]；b1=[-3,2,4]T,b2=[100,-200,345]T,
（2）A行分别为A1=[1,0,8,0.8],A2=[0.8,1,0.8],A3=[0.8,0.8,1]；b1=[3,2,1]T,b2=[5,0,-10]T,
（3）A行分别为A1=[1,3],A2=[-7,1]；b=[4,6]T,
1.程序
（1）雅格比法程序
functionjacobi（）
clc;
clear;
A=[62-1;14-2;-314];
B=[-3;2;4];
Err_user=0.01;
N=500;
[m,n]=size（A）;
X=zeros（n,1）;
k=1;
whilek<=N;
Xk=X;
fori=1:
n
forj=1:
n
ifi~=j
AX（j）=A（i,j）*Xk（j）;
end
end
Sum_AX=sum（AX）;
AX=0;
X（i）=（B（i）-Sum_AX）/A（i,i）;
end
E=max（abs（Xk-X））;
ifEbreak;
end
k=k+1;
end
disp（X）;%显示迭代结果
disp（k）;%显示迭代次数
（2）高斯－赛德尔迭代法程序
functionGS（）
clc;
clear;
A=[62-1;14-2;-314];
B=[-3;2;4];
A=[10.80.8;0.810.8;0.80.81];
B1=[3;2;1];
B2=[5;0;-10];
Err_user=0.0001;
N=50;
[m,n]=size（A）;
X=zeros（n,1）;
k=1;
whilek<=N
Xk=X;
fori=1:
n
forj=1:
n
ifi~=j
AX（j）=A（i,j）*X（j）;%ÓëJacobi·¨Ö÷ÒªÇø±ð
end
end
sum_AX=sum（AX）;
AX=0;X（i）=（B1（i）-sum_AX）/A（i,i）;
end
Er=max（abs（Xk-X））;
ifErbreak;
end
k=k+1;
end
disp（X）;
disp（k）;
end
2.结果分析
（1）
1）A行分别为A1=[6,2,-1],A2=[1,4,-2],A3=[-3,1,4]；b1=[-3,2,4]T
雅克比迭代：
计算结果x=[-0.63630.59600.3737]T,.迭代次数16次。
高斯赛德迭代：
计算结果x=[-0.63630.59590.3738]T,.迭代次数10次。
2）A行分别为A1=[1,0,8,0.8],A2=[0.8,1,0.8],A3=[0.8,0.8,1]；b1=[3,2,1]T
雅克比迭代：
计算不收敛
高斯赛德迭代：
计算结果x=[5.76900.7694-4.2307]T,.迭代次数30次。
（2）
A行分别为A1=[1,3],A2=[-7,1]；b1=[4,6]T。
雅克比迭代：
计算结果不收敛
高斯赛德迭代：
计算结果不收敛
通过对雅克比法和高斯法的上机编程实习，分析对比实验结果可得：
在方程组Ax=b中，右端项对迭代收敛是有影响的，即当b增大时，迭代次数增加，收敛速度降低。
并且通过对比可知，在相同条件下，高斯－赛德尔迭代法比雅克比迭代法收敛速度快；方法的选择也很重要，比如第二问，用雅克比迭代法是发散的，而用高斯迭代法则是收敛的。
通过上机验证，我们也得出结论：
理论分析是正确的，即当迭代矩阵的谱半径小于1时，迭代法是收敛的，而当迭代矩阵谱半径大于1的时，迭代法都是发散的。
第5题
5.用Runge-Kutta4阶算法对初值问题y/=-20*y，y（0）=1按不同步长求解，用于观察稳定区间的作用，推荐两种步长h=0.1,0.2。
注：
此方程的精确解为：
y=e-20x
1.程序
%%%%%%%%%%%%%%%%%%%%%
%Runge-Kutta4阶算法
%f=-20y
%y（0）=1
%%%%%%%%%%%%%%%%%%%%%
clc;
clear;
N=10;%设定节点个数
h=0.05;%设定步长
x
（1）=0;%x0=0
y
（1）=1;%y（0）=1
yy
（1）=exp（-20*x
（1））;%y（0）的精确解
%%%%%%%%%%%%%%%%%%%%%%%开始用runge-kutta法计算
fori=2:
N
K1=-20*（y（i-1））;%（xi,yi）点的导数为f=-20*y
K2=-20*（y（i-1）+K1*h/2）;
K3=-20*（y（i-1）+K2*h/2）;
K4=-20*（y（i-1）+K3*h）;
delta=h*（K1+2*K2+2*K3+K4）/6;
y（i）=y（i-1）+delta;%计算y（i）值
x（i）=x（i-1）+h;%计算下个节点的x（i）值
yy（i）=exp（-20*x（i））;%计算y（i）的精确值
end;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%存储计算结果
dlmwrite（'f:
\计算结果.xls',x,'delimiter','\t','precision',8）;
dlmwrite（'f:
\计算结果.xls',y,'-append','delimiter','\t','precision',8）;
dlmwrite（'f:
\计算结果.xls',yy,'-append','delimiter','\t','precision',8）;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%画图功能
plot（x,y,'o',x,yy,'*'）;
xlabel（'x'）;
ylabel（'y'）;
legend（'四阶Runge-Kutta法','精确解'）;
pause;
close;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2.结果分析
（1）步长h=0.1时结果
当步长取为0.1时，，计算结果如图1.2-2：
图5-1h=0.05时计算结果
计算结果整理如表5-1所示：
表5-1h=0.1时结果
1.2.3步长h=0.2时结果
当步长取为0.2时，，计算结果如图5-2：
图5-2h=0.1时计算结果
计算结果整理如表5-2所示：
图5-2h=0.2时计算结果
通过以上对Runge-Kutta法的应用，计算结果表明步长h的取值会影响算法的收敛性和稳定性：
（1）Runge-Kutta法的步长h越长，计算结果的精确度越低，甚至计算结果不收敛；
（2）当步长h在稳定区间时，误差逐步衰减。
升级会员 