计算方法课程设计..doc
《计算方法课程设计..doc》由会员分享，可在线阅读，更多相关《计算方法课程设计..doc（22页珍藏版）》请在冰豆网上搜索。
数理学院
2014级信息与计算科学
课程设计
姓名：
刘金玉
学号：
3141301240
班级：
1402
成绩：
实验要求
1．应用自己熟悉的算法语言编写程序，使之尽可能具有通用性。
2．上机前充分准备，复习有关算法，写出计算步骤，反复检查，调试程序。
（注：
在练习本上写，不上交）
3．完成计算后写出实验报告，内容包括:
算法步骤叙述，变量说明，程序清单，输出计算结果，结构分析和小结等。
（注：
具体题目具体分析，并不是所有的题目的实验报告都包含上述内容！
）
4．独立完成，如有雷同，一律判为零分！
5．上机期间不允许做其他任何与课程设计无关的事情，否则被发现一次扣10分，被发现三次判为不及格！
非特殊情况，不能请假。
旷课3个半天及以上者，直接判为不及格。
目录
一、基本技能训练 4
1、误差分析 4
2、求解非线性方程 6
3、插值 12
4、数值积分 12
二、提高技能训练 16
1、 16
2、 18
三、本课程设计的心得体会（500字左右） 21
一、基本技能训练
1、误差分析
实验1.3求一元二次方程的根
实验目的：
研究误差传播的原因与解决对策。
问题提出：
求解一元二次方程
实验内容：
一元二次方程的求根公式为
用求根公式求解下面两个方程：
实验要求：
（1）考察单精度计算结果（与真解对比）；
（2）若计算结果与真解相差很大，分析其原因，提出新的算法（如先求再根据根与系数关系求）以改进计算结果。
实验步骤：
方程
（1）：
根据求根公式，写出程序：
formatlong
a=1;b=3;c=-2;
x1=（（-1）*b+sqrt（b^2-4*a*c））/2*a
x2=（（-1）*b-sqrt（b^2-4*a*c））/2*a
运行结果：
x1=0.561552812808830
x2=-3.561552812808830
然后
由符号解求的该方程的真值程序为：
symsx;
y=x^2+3*x-2;
s=solve（y,x）;
vpa（s）
运行结果为：
X1=0.56155281280883027491070492798704
X2=-3.561552812808830274910704927987
方程
（2）：
formatlong
a=1;b=-10^10;c=1;
x1=（（-1）*b+sqrt（b^2-4*a*c））/2*a
x2=（（-1）*b-sqrt（b^2-4*a*c））/2*a
运行结果为：
x1=1.000000000000000e+010
x2=0
然后
由符号解求的该方程的真值程序为：
symsx;
y=x^2-10^10*x+1;
s=solve（y,x）;
vpa（s）
运行结果：
X1=10000000000.0
X2=0.000000000116415321827
由此可知，对于方程
（1），使用求根公式求得的结果等于精确值，故求根公式法可靠。
而对于方程
（2），计算值与真值相差很大，故算法不可靠。
改进算法，对于方程
（2）:
先用迭代法求x1，再用，求x2
程序为：
symsk
a=[];
fori=1:
100
a
（1）=4;
a（i+1）=（10^10*a（i）-1）^（1/2）;
end
x1=a（100）
x2=1/（x1）
运行结果为：
x1=1.000000000000000e+010
x2=1.000000000000000e-010
实验结论：
对于方程
（1），两种方法在精确到小数点后15位时相同，说明两种算法的结果都是精确的。
对于方程
（2），两种算法结果有相当大的偏差，求根公式所求的一个根直接为零，求根公式的算法是不精确的。
原因：
方程
（2）用求根公式计算时，公式中，b是大数，出现了大数吃掉小数的误差，也出现了两个相近的数相减的误差，所以出现x2=0这样大的误差。
改进的结果会比较准确。
2、求解非线性方程
实验2.1Gauss消去法的数值稳定性实验
实验目的：
观察和理解高斯消元过程中出现小主元即很小时，引起方程组解的数值不稳定性．
实验内容：
求解线性方程组其中
（1），
（2）
实验要求：
（1）计算矩阵的条件数，判断系数矩阵是良态的还是病态的
（2）用高斯列主元消去法求得L和U及解向量
（3）用不选主元的高斯消去法求得L和U及解向量
（4）观察小主元并分析对计算结果的影响
实验步骤：
（1）计算矩阵的条件数程序：
矩阵A1：
A1=[0.3*10^（-15）59.1431;5.291-6.130-12;11.2952;1211];
cond（A1,1）
cond（A1,2）
cond（A1,inf）
运行结果：
ans=136.294470872045e+000
ans=68.4295577125341e+000
ans=84.3114602051800e+000
由矩阵条件数判断出矩阵A1是病态矩阵。
矩阵A2：
A2=[10-701;-32.0999999999962;5-15-1;0102];
cond（A1,1）
cond（A1,2）
cond（A1,inf）
运行结果：
ans=19.2831683168042e+000
ans=8.99393809015365e+000
ans=18.3564356435280e+000
由矩阵条件数判断出矩阵A2是病态矩阵。
（2）高斯列主元消去法程序：
方程组
（1）：
A1=[0.3*10^（-15）59.1431;5.291-6.130-12;11.2952;1211];
b1=[59.17;46.78;1;2];
n=4;
fork=1:
n-1
a=max（abs（A1（k:
n,k）））;
[p,k]=find（A1==a）;
B=A1（k,:
）;c=b1（k）;
A1（k,:
）=A1（p,:
）;b1（k）=b1（p）;
A1（p,:
）=B;b1（p）=c;
ifA1（k,k）~=0
A1（k+1:
n,k）=A1（k+1:
n,k）/A1（k,k）;
A1（k+1:
n,k+1:
n）=A1（k+1:
n,k+1:
n）-A1（k+1:
n,k）*A1（k,k+1:
n）;
else
break
end
end
L1=tril（A1,0）;
fori=1:
n
L1（i,i）=1;
end
L=L1
U=triu（A1,0）
forj=1:
n-1
b1（j）=b1（j）/L（j,j）;
b1（j+1:
n）=b1（j+1:
n）-b1（j）*L（j+1:
n,j）;
end
b1（n）=b1（n）/L（n,n）;
forj=n:
-1:
2
b1（j）=b1（j）/U（j,j）;
b1（1:
j-1）=b1（1:
j-1）-b1（j）*U（1:
j-1,j）;
end
b1
（1）=b1
（1）/U（1,1）;
x1=b1
运行结果：
方程组
（2）：
A2=[10-701;-32.0999999999962;5-15-1;0102];
b2=[8;5.9000000000001;5;1];
n=4;
fork=1:
n-1
a=max（abs（A2（k:
n,k）））;
[p,k]=find（A2==a）;
B=A2（k,:
）;c=b2（k）;
A2（k,:
）=A2（p,:
）;b2（k）=b2（p）;
A2（p,:
）=B;b2（p）=c;
ifA2（k,k）~=0
A2（k+1:
n,k）=A2（k+1:
n,k）/A2（k,k）;
A2（k+1:
n,k+1:
n）=A2（k+1:
n,k+1:
n）-A2（k+1:
n,k）*A2（k,k+1:
n）;
else
break
end
end
L1=tril（A2,0）;
fori=1:
n
L1（i,i）=1;
end
L=L1
U=triu（A2,0）
forj=1:
n-1
b2（j）=b2（j）/L（j,j）;
b2（j+1:
n）=b2（j+1:
n）-b2（j）*L（j+1:
n,j）;
end
b2（n）=b2（n）/L（n,n）;
forj=n:
-1:
2
b2（j）=b2（j）/U（j,j）;
b2（1:
j-1）=b2（1:
j-1）-b2（j）*U（1:
j-1,j）;
end
b2
（1）=b2
（1）/U（1,1）;
x2=b2
运行结果：
（3）不选主元的高斯消去法程序：
方程组
（1）：
clear
formatlong
A1=[0.3e-1559.1431;5.291-6.130-12;11.2952;1211];
b1=[59.17;46.78;1;2];
n=4;
fork=1:
n-1
A1（k+1:
n,k）=A1（k+1:
n,k）/A1（k,k）;
A1（k+1:
n,k+1:
n）=A1（k+1:
n,k+1:
n）-A1（k+1:
n,k）*A1（k,k+1:
n）;
end
L1=tril（A1,0）;
fori=1:
n
L1（i,i）=1;
end
L=L1
U=triu（A1,0）
forj=1:
n-1
b1（j）=b1（j）/L（j,j）;
b1（j+1:
n）=b1（j+1:
n）-b1（j）*L（j+1:
n,j）;
end
b1（n）=b1（n）/L（n,n）;
forj=n:
-1:
2
b1（j）=b1（j）/U（j,j）;
b1（1:
j-1）=b1（1:
j-1）-b1（j）*U（1:
j-1,j）;
end
b1
（1）=b1
（1）/U（1,1）;
x1=b1
运行结果:
方程组
（2）：
clear
formatlong
A2=[10-701;-32.0999999999962;5-15-1;0102];
b2=[8;5.9000000000001;5;1];
n=4;
fork=1:
n-1
A2（k+1:
n,k）=A2（k+1:
n,k）/A2（k,k）;
A2（k+1:
n,k+1:
n）=A2（k+1:
n,k+1:
n）-A2（k+1:
n,k）*A2（k,k+1:
n）;
end
L1=tril（A2,0）;
fori=1:
n
L1（i,i）=1;
end
L=L1
U=triu（A2,0）
forj=1:
n-1
b2（j）=b2（j）/L（j,j）;
b2（j+1:
n）=b2（j+1:
n）-b2（j）*L（j+1:
n,j）;
end
b2（n）=b2（n）/L（n,n）;
forj
升级会员 