数理学院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 求一元二次方程的根实验目的：研究误差传播的原因与解决对策。

问题提出：求解一元二次方程20ax bx c ++=实验内容：一元二次方程的求根公式为1,22b x a-+= 用求根公式求解下面两个方程：2210(1)320(2)1010x x x x +-=-+=实验要求：（1） 考察单精度计算结果（与真解对比）；（2） 若计算结果与真解相差很大，分析其原因，提出新的算法（如先求1x 再根据根与系数关系求2x ）以改进计算结果。

实验步骤：方程（1）：根据求根公式，写出程序：format longa=1;b=3;c=-2;x1=((-1)*b+sqrt(b^2-4*a*c))/2*ax2=((-1)*b-sqrt(b^2-4*a*c))/2*a运行结果：x1 = 0.5630x2 = -3.5630然后由符号解求的该方程的真值程序为：syms x;y=x^2+3*x-2;s=solve(y,x);vpa(s)运行结果为：X1= 0.5632798704X2= -3.56327987方程（2）：format longa=1;b=-10^10;c=1;x1=((-1)*b+sqrt(b^2-4*a*c))/2*ax2=((-1)*b-sqrt(b^2-4*a*c))/2*a运行结果为：x1 = 1.0000e+010x2 = 0然后由符号解求的该方程的真值程序为：syms x;y=x^2-10^10*x+1;s=solve(y,x);vpa(s)运行结果：X1= 10000000000.0X2= 0.6415321827由此可知，对于方程（1），使用求根公式求得的结果等于精确值，故求根公式法可靠。

而对于方程（2），计算值与真值相差很大，故算法不可靠。

改进算法，对于方程（2）:先用迭代法求x1，再用11ax x，求x2c程序为：syms ka=[ ];for i=1:100a(1)=4;a(i+1)=(10^10*a(i)-1)^(1/2);endx1=a(100)x2=1/(x1)运行结果为：x1 = 1.0000e+010x2 = 1.0000e-010实验结论：对于方程（1），两种方法在精确到小数点后15位时相同，说明两种算法的结果都是精确的。

对于方程（2），两种算法结果有相当大的偏差，求根公式所求的一个根直接为零，求根公式的算法是不精确的。

原因：方程（2）用求根公式计算时，b -b 是大数，出现了大数吃掉小数的误差，也出现了两个相近的数相减的误差，所以出现x2=0这样大的误差。

改进的结果会比较准确。

2、求解非线性方程实验2.1 Gauss 消去法的数值稳定性实验实验目的：观察和理解高斯消元过程中出现小主元即()k kk a 很小时，引起方程组解的数值不稳定性．实验内容：求解线性方程组(1,2)i i A x b i == 其中（1）1510.31059.14315.291 6.1301211.29521211A -⎛⎫⨯ ⎪-- ⎪= ⎪ ⎪ ⎪⎝⎭，159.1746.7812b ⎛⎫ ⎪ ⎪= ⎪ ⎪⎝⎭（2）210701 3 2.09999999999962 5151 0102A-⎛⎫⎪- ⎪=⎪--⎪⎝⎭185.900000000000151b⎛⎫⎪⎪=⎪⎪⎝⎭实验要求：（1）计算矩阵的条件数，判断系数矩阵是良态的还是病态的（2）用高斯列主元消去法求得L和U及解向量（3）用不选主元的高斯消去法求得L和U及解向量（4）观察小主元并分析对计算结果的影响实验步骤：（1）计算矩阵的条件数程序：矩阵A1：A1=[0.3*10^(-15) 59.14 3 1;5.291 -6.130 -1 2;11.2 9 5 2;1 2 1 1];cond(A1,1)cond(A1,2)cond(A1,inf)运行结果：ans = 136.294470872045e+000ans = 68.4295577125341e+000ans = 84.3114602051800e+000由矩阵条件数判断出矩阵A1是病态矩阵。

矩阵A2：A2=[10 -7 0 1;-3 2. 6 2;5 -1 5 -1;0 1 0 2];cond(A1,1)cond(A1,2)cond(A1,inf)运行结果：ans = 19.2831683168042e+000ans = 8.99393809015365e+000ans = 18.3564356435280e+000由矩阵条件数判断出矩阵A2是病态矩阵。

（2）高斯列主元消去法程序：方程组（1）：A1=[0.3*10^(-15) 59.14 3 1;5.291 -6.130 -1 2;11.2 9 5 2;1 2 1 1];b1=[59.17;46.78;1;2];n=4;for k=1:n-1a=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;if A1(k,k)~=0A1(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);elsebreakendendL1=tril(A1,0);for i=1:nL1(i,i)=1;endL=L1U=triu(A1,0)for j=1:n-1b1(j)=b1(j)/L(j,j);b1(j+1:n)=b1(j+1:n)-b1(j)*L(j+1:n,j);endb1(n)=b1(n)/L(n,n);for j=n:-1:2b1(j)=b1(j)/U(j,j);b1(1:j-1)=b1(1:j-1)-b1(j)*U(1:j-1,j);endb1(1)=b1(1)/U(1,1);x1=b1运行结果：18100026.79101000.47240.1755100.08930.02020.49291L ⎡⎤⎢⎥⨯⎢⎥=⎢⎥-⎢⎥--⎣⎦ 11.2952059.143100 2.835 1.2310000.801U ⎡⎤⎢⎥⎢⎥=⎢⎥-⎢⎥⎣⎦ 1[18.9882;3.3378;34.747;33.9865]x =--方程组（2）：A2=[10 -7 0 1;-3 2. 6 2;5 -1 5 -1;0 1 0 2];b2=[8;5.91;5;1];n=4;for k=1:n-1a=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;if A2(k,k)~=0A2(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); elsebreakendendL1=tril(A2,0);for i=1:nL1(i,i)=1;endL=L1U=triu(A2,0)for j=1:n-1b2(j)=b2(j)/L(j,j);b2(j+1:n)=b2(j+1:n)-b2(j)*L(j+1:n,j);endb2(n)=b2(n)/L(n,n);for j=n:-1:2b2(j)=b2(j)/U(j,j);b2(1:j-1)=b2(1:j-1)-b2(j)*U(1:j-1,j);endb2(1)=b2(1)/U(1,1);x2=b2运行结果：1210000.51000.30.4101000.40.3331L ⎡⎤⎢⎥⎢⎥=⎢⎥--⨯⎢⎥-⎣⎦ 107010 2.55 1.5006 2.30000 3.3667U -⎡⎤⎢⎥-⎢⎥=⎢⎥⎢⎥⎣⎦ 112[0.44410; 1.0;1.0;1.0]x -=⨯-（3） 不选主元的高斯消去法程序：方程组（1）：clearformat longA1=[0.3e-15 59.14 3 1;5.291 -6.130 -1 2;11.2 9 5 2;1 2 1 1]; b1=[59.17;46.78;1;2];n=4;for k=1:n-1A1(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); endL1=tril(A1,0);for i=1:nL1(i,i)=1;endL=L1U=triu(A1,0)for j=1:n-1b1(j)=b1(j)/L(j,j);b1(j+1:n)=b1(j+1:n)-b1(j)*L(j+1:n,j);endb1(n)=b1(n)/L(n,n);for j=n:-1:2b1(j)=b1(j)/U(j,j);b1(1:j-1)=b1(1:j-1)-b1(j)*U(1:j-1,j); endb1(1)=b1(1)/U(1,1); x1=b1 运行结果:151515100017.63710100'37.3310 2.1168103.333100.18900L ⎡⎤⎢⎥⨯⎢⎥=⎢⎥⨯⎢⎥⨯⎣⎦151815150.31059.14310 1.0431052.911017.63710'001680000.5U ⎡⎤⨯⎢⎥-⨯-⨯-⨯⎢⎥=⎢⎥--⎢⎥⎣⎦'1[23.6848;1.0005;0;0]x =方程组（2）： clear format longA2=[10 -7 0 1;-3 2. 6 2;5 -1 5 -1;0 1 0 2]; b2=[8;5.91;5;1]; n=4; for k=1:n-1A2(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); endL1=tril(A2,0); for i=1:n L1(i,i)=1; end L=L1U=triu(A2,0)for j=1:n-1b2(j)=b2(j)/L(j,j);b2(j+1:n)=b2(j+1:n)-b2(j)*L(j+1:n,j); endb2(n)=b2(n)/L(n,n); for j=n:-1:2b2(j)=b2(j)/U(j,j);b2(1:j-1)=b2(1:j-1)-b2(j)*U(1:j-1,j); endb2(1)=b2(1)/U(1,1); x2=b2 运行结果：121210000.3100'0.5 2.4998101000.9999100.4001L ⎡⎤⎢⎥-⎢⎥=⎢⎥-⨯⎢⎥-⨯⎣⎦121212107010 1.01062.3'0014.998710 5.749510000.400 3.3667U -⎡⎤⎢⎥-⨯⎢⎥=⎢⎥⨯⨯⎢⎥⎣⎦'52[0.36210; 1.0;1.0;1.0]x -=⨯-（4） 分析小元对计算结果的影响通过观察计算结果，分析可知，小元对计算结果的影响很大，小元的存在会使得到的计算结果有很大的误差。