数值分析报告.docx
《数值分析报告.docx》由会员分享，可在线阅读，更多相关《数值分析报告.docx（26页珍藏版）》请在冰豆网上搜索。
数值分析报告
《数值分析》
课程名称数值分析实验报告
专业班级应用数学111班
姓名陈伟
学号8
教学教师岳雪芝
实验二函数逼近与曲线拟合
一、问题提出
从随机的数据中找出其规律性，给出其近似表达式的问题，在生产实践和科学实验中大量存在，通常利用数据的最小二乘法求得拟合曲线。
在某冶炼过程中，根据统计数据的含碳量与时间关系，试求含碳量与时间t的拟合曲线。
t（分）
0510152025303540455055
01.272.162.863.443.874.154.374.514.584.024.64
二、要求
1、用最小二乘法进行曲线拟合；
2、近似解析表达式为
；
3、打印出拟合函数
，并打印出
与
的误差，
；
4、另外选取一个近似表达式，尝试拟合效果的比较；
5、*绘制出曲线拟合图。
三、目的和意义
1、掌握曲线拟合的最小二乘法；
2、最小二乘法亦可用于解超定线代数方程组；
3、探索拟合函数的选择与拟合精度间的关系。
四、实验步骤:
输入代码：
t=[0510152025303540455055];
y=[01.272.162.863.443.874.154.374.514.584.024.64];
fori=1:
12
A（i,1）=t（i）;
A（i,2）=t（i）.^2;
A（i,3）=t（i）.^3;
end
p_star=A\y';
y_star=p_star
（1）*t+p_star
（2）*t.^2+p_star（3）*t.^3;
s_star=0;
fori=1:
12
s_star=s_star+（y_star（i）-y（i））^2;
end
plot（t,y,'*',t,y_star,'g'）;
比较拟和结果，分别取二次，三次和四次多项式，再做最小二乘法，代码如下：
t=[0510152025303540455055];
y=[01.272.162.863.443.874.154.374.514.584.024.64];
%3次拟合
p3=polyfit（t,y,3）
y_new=polyval（p3,t）;
s_3=0;
fori=1:
12
s_3=s_3+（y_new（i）-y（i））^2;
end
%2次拟合
p2=polyfit（t,y,2）;
y_new2=polyval（p2,t）;
s_2=0;
fori=1:
12
s_2=s_2+（y_new2（i）-y（i））^2;
end
%4次拟合
p4=polyfit（t,y,4）;
y_new4=polyval（p4,t）;
s_4=0;
fori=1:
12
s_4=s_4+（y_new4（i）-y（i））^2;
end
比较4种拟合函数，结论：
常数项为0的三次拟合函数在保证拟合精度的同时，保证0点的函数值仍为0，是四种拟合方法中最好的一种。
原图与曲线拟合图如下所示：
分析:
从上面的拟合中也可以知道多项式拟合误差平方和随着拟合多项式次数的增加而逐渐减少，拟合曲线更靠近实际数据，拟合更准确。
实验三数值积分与数值微分
一、基本题
选用复合梯形公式，复合Simpson公式，Romberg算法，计算
（1）
二、要求
1、编制数值积分算法的程序；
2、分别用两种算法计算同一个积分，并比较其结果；
3、分别取不同步长
，试比较计算结果（如n=10,20等）；
4、给定精度要求ε，试用变步长算法，确定最佳步长。
三、目的和意义
1、深刻认识数值积分法的意义；
2、明确数值积分精度与步长的关系；
3、根据定积分的计算方法，结合专业考虑给出一个二重积分的计算问题。
四、实验步骤：
代码1，复合梯形公式
functionresult=trapezoid_integration（f,count_node,a,b）
%compositetrapezoidintegration,f,函数表达式,count_node插入节点数量,a,开始点,b结束点
h=（b-a）/count_node;%h步长
np1=count_node+1;%包括端点,总共的点的数量
x=a:
h:
b;%x自变量
c=ones（np1,1）;
c（1,1）=0.5;c（np1,1）=0.5;%两个端点都取一半
symssymbol_x;
fx=subs（f,symbol_x,x）;
result=vpa（h*fx*c,9）;
代码2，复合simpson公式
functionresult=simpson_integration（f,count_node,a,b）
%compositesimpsonintegration,f,函数表达式,count_node插入节点数量,a,开始点,b结束点
h=（b-a）/count_node;%h步长
np1=2*count_node+1;%包括端点,总共的2n+1点
x=a:
h/2:
b;%x自变量
c=ones（np1,1）;
fori=1:
np1
if（mod（i,2）==0）%偶数
c（i,1）=2/3;
else%奇数
c（i,1）=1/3;
end
end
c（1,1）=1/6;c（np1,1）=1/6;%两个端点都取1/6
symssymbol_x;
fx=subs（f,symbol_x,x）;
result=vpa（h*fx*c,9）;
代码3，检验复合梯形公式和复合simpson公式
function[result1]=test_integration（）
%测试复合梯形公式，复合simpson公式,用到trapezoid_integration,simpson_integration
symssymbol_x;
f1=sin（symbol_x）/symbol_x;
symsA1;%用于输出f1积分结果
fori=1:
9
A1（i,1）=trapezoid_integration（f1,i,1e-9,1）;
A1（i,2）=simpson_integration（f1,i,1e-9,1）;
end
k=10;
fori=10:
10:
100
A1（k,1）=trapezoid_integration（f1,i,1e-9,1）;
A1（k,2）=simpson_integration（f1,i,1e-9,1）;
k=k+1;
end
IS
（2）=int（f1,symbol_x,0,1）;
vpa（IS
（2））
result1=A1;
结果分析
从上面的计算结果可以看出复合梯形公式和复合simpson公式的稳定性，并且步长越小精度越高，当n趋于正无穷，即步长趋于0时，上述两公式的计算值都收敛到积分值。
从收敛的速度来看，复合simpson公式优于复合梯形公式。
实验四线方程组的直接解法
一、问题提出：
1、三对角形线性方程组
二、要求
1、对上述三个方程组分别利用Gauss顺序消去法与Gauss列主元消去法；平方根法与改进平方根法；追赶法求解（选择其一）；
2、应用结构程序设计编出通用程序；
3、比较计算结果，分析数值解误差的原因；
4、尽可能利用相应模块输出系数矩阵的三角分解式。
三、目的和意义
1、通过该课题的实验，体会模块化结构程序设计方法的优点；
2、运用所学的计算方法，解决各类线性方程组的直接算法；
3、提高分析和解决问题的能力，做到学以致用；
通过三对角形线性方程组的解法，体会稀疏线性方程组解法的特点。
四、实验步骤：
高斯顺序消去法:
functionresult=Gauss（A,d）
%高斯顺序消元法
[n,m]=size（A）
k=ones（1,n）;%用于记录-A（j,i）/A（i,i）
fori=1:
n
if（A（i,i）==0）
'主元为0，不能使用高斯顺序消元法'
return;
end
forj=i+1:
n
k（j）=-A（j,i）/A（i,i）
A（j,:
）=A（i,:
）.*k（j）+A（j,:
）;%加到第j行
d（j）=d（i）.*k（j）+d（j）;
end
end
x=zeros（n,1）
x（n）=d（n）/A（n,n）;
fori=n-1:
-1:
1
x（i）=（d（i）-A（i,:
）*x）/A（i,i）;
end
result=x
Gauss列主元消去法
functionresult=Gauss2（A,d）
%高斯列主元消元法
[n,m]=size（A）;
[d_length,d_width]=size（d）;
if（d_length~=n||d_width~=1）
'方程右端向量输入有误'
return;
end
if（n~=m）
'系数矩阵不是方阵'
return;
end
AandD=[Ad];%增广阵
k=ones（1,n）;%用于记录-A（j,i）/A（i,i）
record=1:
n;%用于记录行的交换情况
fori=1:
n
%选取第i列的最大元进行换位
[l,index]=max（abs（AandD（i:
n,i）））;
if（index~=i）
temp_record=record（index）;
record（index）=record（i+index-1）;
record（i+index-1）=temp_record;
temp=AandD（i,:
）;
AandD（i,:
）=AandD（i+index-1,:
）;
AandD（i+index-1,:
）=temp;
end
%顺序消元
if（AandD（i,i）==0）
'主元为0，不能使用高斯列主元顺序消元法'
return;
end
forj=i+1:
n
k（j）=-AandD（j,i）/AandD（i,i）;
AandD（j,:
）=AandD（i,:
）.*k（j）+AandD（j,:
）;%加到第j行
end
end
A=AandD（1:
n,1:
n）;
d=AandD（1:
n,n+1）;
x=zeros（n,1）;
x（n）=d（n）/A（n,n）;
fori=n-1:
-1:
1
x（i）=（d（i）-A（i,:
）*x）/A（i,i）;
end
%由于不是顺序消元，所以还要再换回来
fori=1:
n
result（i）=x（record（i））;
end
追赶法
functionresult=chase_after（A,d）
[n,m]=size（A）;
[d_length,d_width]=size（d）;
if（d_length~=n||d_width~=1）
'方程右端向量输入有误'
return;
end
if（n~=m）
'系数矩阵不是方阵'
return;
end
fori=2:
n
toCheck1=diag（A,i）;
toCheck2=diag（A,-i）;if（norm（toCheck1）~=0||norm（toCheck2）~=0）
'系数矩阵不是三对角阵'
return;
end
end
b=diag（A）;
a=zeros（n,1）;
a（2:
n）=diag（A,-1）;
c=diag（A,1）;
if（b
（1）==0）
'因顺序主子式为0,不能进行Crout分解'
return;
end
l
（1）=b
（1）;
y
（1）=d
（1）/l
（1）;
u
（1）=c
（1）/l
（1）;
fori=2:
n
l（i）=b（i）-a（i）*u（i-1）;
if（l（i）==0）
'因顺序主子式为0,不能进行Crout分解'
return;
end
y（i）=（d（i）-y（i-1）*a（i））/l（i）;
if（i~=n）
u（i）=c（i）/l（i）;
end
end
x（n）=y（n）;
fori=n-1:
-1:
1
x（i）=y（i）-u（i）*x（i+1）;
end
y=y';result=x';
检验三种直接解线性方程组的方法
functiontest_chase_after
A=[4-100000000;
-14-10000000;
0-14-1000000;
00-14-100000;
000-14-10000;
0000-14-1000;
00000-14-100;
000000-14-10;
0000000-14-1;
00000000-14;];
d=[75-1326-1214-45-5]';%按列存储
[L,U]=lu（A）
%检验高斯顺序消去法
result_Gauss=Gauss（A,d）
%检验高斯列主元消去法
result_Gauss2=Gauss2（A,d）
%检验追赶法
result_chase=chase_after（A,d）
运行结果：
高斯顺序消去法的计算结果：
result_Gauss=
2
1
-3
1.7849e-016
1
-2
3
-1.4874e-016
1
-1
高斯列主元消去法的计算结果：
result_Gauss=
2
1
-3
1.7849e-016
1
-2
3
-1.4874e-016
1
-1
追赶法的计算结果：
result_Gauss=
2
1
-3
1.7849e-016
1
-2
3
-1.4874e-016
1
-1
结果分析：
高斯顺序消去法与精确解
，比较其误差是1.7849*10^（-16）,该误差可以忽略，事实上这是由计算机的除法引起的。
而在这个问题中，因为顺序消元的除法中用到的数k的绝对值都小于1，这就避免使用绝对值较大的数，使得总体舍入误差得到有效控制，所以高斯顺序消元法的效果也很好。
同样高斯列主元消去法与精确解的误差都是
（1）式的结果，说明在这个问题中误差都控制得很好。
实验五解线性方程组的迭代法
一、问题提出
对实验四所列目的和意义的线性方程组，试分别选用Jacobi迭代法，Gauss-Seidel迭代法和SOR方法计算其解。
二、要求
1、体会迭代法求解线性方程组，并能与消去法做以比较；
2、分别对不同精度要求，如
由迭代次数体会该迭代法的收敛快慢；
3、对方程组2，3使用SOR方法时，选取松弛因子ω=0.8，0.9，1，1.1，1.2等，试看对算法收敛性的影响，并能找出你所选用的松弛因子的最佳者；
4、给出各种算法的设计程序和计算结果。
三、目的和意义
1、通过上机计算体会迭代法求解线性方程组的特点，并能和消去法比较；
2、运用所学的迭代法算法，解决各类线性方程组，编出算法程序；
3、体会上机计算时，终止步骤
或k>（给予的迭代次数），对迭代法敛散性的意义；
4、体会初始解
，松弛因子的选取，对计算结果的影响。
四、实验步骤：
Jacobi迭代法
functionresult=Jacobi（A,b,x0,yupsilong）
%Jacobi迭代方法,x0是迭代用的初始向量,yupsilong是解的误差上限
[n,m]=size（A）;
[b_length,b_width]=size（b）;
if（b_length~=n||b_width~=1）
'方程右端向量输入有误'
return;
end
if（n~=m）
'系数矩阵不是方阵'
return;
end
if（det（A）==0）
'系数矩阵是奇异阵'
return
end
%分解A=D-L-U;D是对角元，L是对角元为0的下三角阵，U是对角元为0的上三角阵
L=tril（A,-1）;
U=triu（A,1）;
D=A-L-U;
L=-L;
U=-U;
%计算迭代矩阵
B0=D^（-1）;
B=B0*（L+U）;
g=B0*b;
x=x0;
delta=1;
count_iteration=0;
whiledelta>=yupsilong
if（count_iteration>1000）
'超过迭代次数上限，认为算法此时不收敛'
return;
end
x1=B*x+g;
delta=norm（x1-x）;
x=x1;
count_iteration=count_iteration+1;
end
count_iteration
result=x;
Gauss_Seidol迭代法
functionresult=Gauss_Seidol（A,b,x0,yupsilong）
%Gauss_Seidol方法,x0是迭代用的初始向量,yupsilong是解的误差上限
[n,m]=size（A）;
[b_length,b_width]=size（b）;
if（b_length~=n||b_width~=1）
'方程右端向量输入有误'
return;
end
if（n~=m）
'系数矩阵不是方阵'
return;
end
if（det（A）==0）
'系数矩阵是奇异阵'
return
end
%分解A=D-L-U;D是对角元，L是对角元为0的下三角阵，U是对角元为0的上三角阵
L=tril（A,-1）;
U=triu（A,1）;
D=A-L-U;
L=-L;
U=-U;
%计算迭代矩阵
B0=（D-L）^（-1）;
B=B0*U;
g=B0*b;
x=x0;
delta=1;
count_iteration=0;
whiledelta>=yupsilong
if（count_iteration>1000）
'超过迭代次数上限，认为算法此时不收敛'
return;
end
x1=B*x+g;
delta=norm（x1-x）;
x=x1;
count_iteration=count_iteration+1;
end
count_iteration
result=x;
SOR迭代法
functionresult=SOR（A,b,x0,yupsilong,omiga）
%SOR迭代方法,x0是迭代用的初始向量,yupsilong是解的误差上限,omiga是松弛因子
[n,m]=size（A）;
[b_length,b_width]=size（b）;
if（b_length~=n||b_width~=1）
'方程右端向量输入有误'
return;
end
if（n~=m）
'系数矩阵不是方阵'
return;
end
if（det（A）==0）
'系数矩阵是奇异阵'
return
end
%分解A=D-L-U;D是对角元，L是对角元为0的下三角阵，U是对角元为0的上三角阵
L=tril（A,-1）;
U=triu（A,1）;
D=A-L-U;
L=-L;
U=-U;
%计算迭代矩阵
B0=（D-omiga*L）^（-1）;
B=B0*（（1-omiga）*D+omiga*U）;
g=omiga*B0*b;
x=x0;
delta=1;
count_iteration=0;
whiledelta>=yupsilong
if（count_iteration>1000）
'超过迭代次数上限，认为算法此时不收敛'
return;
end
x1=B*x+g;
delta=norm（x1-x）;
x=x1;
count_iteration=count_iteration+1;
end
count_iteration
result=x;
检验Jocobi,GS,SOR迭代方法
functiontest_Gauss_Seidol
clc;
A=[4-100000000;
-14-10000000;
0-14-1000000;
00-14-100000;
000-14-10000;
0000-14-1000;
00000-14-100;
000000-14-10;
0000000-14-1;
00000000-14;];
b=[75-1326-1214-45-5]';
[n,m]=size（A）;
yupsilong=[1e-31e-41e-51e-61e-71e-81e-91e-10];
x_real=[21-301-2301-1];%精确解
'测试初始点对迭代的影响'
fori=1:
8
x0=ones（n,1）*10^（i-1）;
x0'
result_Jacobi=（Jacobi（A,b,x0,yupsilong（3）））';
if（norm（result_Jacobi-x_real）>yupsilong（3））
'误差过大'
end
result_GS=（Gauss_Seidol（A,b,x0,yupsilong（3）））';
if（norm（result_GS-x_real）>yupsilong（3））
'误差过大'
end
result_SOR=（SOR（A,b,x0,yupsilong（3）,1.1））';
if（norm（result_GS-x_real）>yupsilong（3））
'误差过大'
end
end
'测试yupsilong对迭代的影响'
x0=ones（n,1）;
fori=1:
length（yupsilong）
（yupsilong（i））'
result_Jacobi=（Jacobi（A,b,x0,yupsilong（i）））';
result_GS=（Gauss_Seidol（A,b,x0,yup
升级会员 