数值分析上机实验报告课程名称：数值分析上机实验学院：机械工程学院专业：机械制造姓名：张法光学号：2012021691 年级：12级任课教师：代新敏老师2012年12月30日一．已知A 与b12.38412 2.115237 -1.061074 1.112336 -0.1135840.718719 1.742382 3.067813 -2.031743 2.11523719.141823-3.125432 -1.012345 2.1897361.563849-0.784165 1.112348 3.123124 -1.061074 -3.125A =43215.567914 3.123848 2.031454 1.836742-1.056781 0.336993 -1.010103 1.112336 -1.012345 3.12384827.108437 4.101011-3.741856 2.101023 -0.71828 -0.037585 -0.1135842.189736 2.031454 4.10101119.8979180.431637-3.111223 2.121314 1.784317 0.718719 1.563849 1.836742 -3.741856 0.4316379.789365-0.103458 -1.103456 0.238417 1.742382 -0.784165 -1.056781 2.101023-3.111223-0.10345814.7138465 3.123789 -2.213474 3.067813 1.112348 0.336993-0.71828 2.121314-1.103456 3.12378930.719334 4.446782 -2.031743 3.123124 -1.010103-0.037585 1.7843170.238417-2.213474 4.44678240.00001[ 2.1874369 33.992318 -25.173417 0.84671695 1.784317 -86.612343 1.1101230 4.719345 -5.6784392]TB ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦=（2）用超松弛法求解Bx=b （取松弛因子ω=1.4，x （0）=0，迭代9次）。

（3）用列主元素消去法求解Bx=b 。

解：（3）、用列主元素消去法求解Bx=b（一）、理论依据：其基本思想是选取绝对值尽量大的元素作为主元素，进行行与列的交换，再进行回代，求出方程的解。

将方阵A 和向量b 写成C=（A b ）。

将C 的第1列中第1行的元素与其下面的此列的元素逐一进行比较，找到最大的元素1j c ，将第j 行的元素与第1行的元素进行交换，然后通过行变换，将第1列中第2到第n 个元素都消成0。

将变换后的矩阵(1)C 的第二列中第二行的元素与其下面的此列的元素逐一进行比较，找到最大的元素(1)2k c ，将第k 行的元素与第2行的元素进行交换，然后通过行变换，将第2列中第3到第n 个元素都消成0。

以此方法将矩阵的左下部分全都消成0。

（二）、计算程序： #include "math.h" #include "stdio.h" void main() { double u[9],x1[9],y[9],q[9],b1[9][10],x[9],a[9][9]={{12.38412,2.115237,-1.061074,1.112336,-0.113584,0.718719,1.742382,3.067813,-2.031743},{2.115237,19.141823,-3.125432,-1.012345,2.189736,1.563849,-0.784165,1.112348,3.123124},{-1.061074,-3.125432,15.567914,3.123848,2.031454,1.836742,-1.056781,0.336993,-1.010103}, {1.112336,-1.012345,3.123848,27.108437,4.101011,-3.741856,2.101023,-0.71828,-0.037585}, {-0.113584,2.189736,2.031454,4.101011,19.897918,0.431637,-3.111223,2.121314,1.784317}, {0.718719,1.563849,1.836742,-3.741856,0.431637,9.789365,-0.103458,-1.103456,0.238417}, {1.742382,-0.784165,-1.056781,2.101023,-3.111223,-0.103458,14.713846,3.123789,-2.213474}, {3.067813,1.112348,0.336993,-0.71828,2.121314,-1.103456,3.123789,30.719334,4.446782},{-2.031743,3.123124,-1.010103,-0.037585,1.784317,0.238417,-2.213474,4.446782,40.00001}};int sign(double x);double k,t,s,w,e,c,z;int i,j,n,r;doubleb[9]={2.1874369,33.992318,-25.173417,0.84671695,1.784317,-86.612343,1.1101230,4.719345,-5.6784392};for(r=0;r<=6;r++) /*Household 变换*/{e=0.0;for(i=r+1;i<=8;i++)e=e+a[i][r]*a[i][r];s=sqrt(e);t=s*s+fabs(a[r+1][r])*s;for(i=0;i<=r;i++) u[i]=0;c=a[r+1][r]; /*求u[i]的值*/u[r+1]=a[r+1][r]+s*sign(c);for(i=r+2;i<=8;i++)u[i]=a[i][r];for(i=0;i<=8;i++){y[i]=0;for(j=0;j<=8;j++)y[i]+=a[i][j]*u[j]/t;} /*求出y向量*/ k=0;for(i=0;i<9;++i)k+=0.5*(u[i]*y[i])/t;for(i=0;i<=8;i++)q[i]=y[i]-k*u[i];for(i=0;i<=8;i++)for(j=0;j<=8;j++)a[i][j]=a[i][j]-(q[j]*u[i]+u[j]*q[i]);} /*求结果*/printf("Household变换:\n");for(i=0;i<9;++i)for(j=0;j<9;++j) /*打印转化后的矩阵*/ {if (j%9==0)printf("\n");printf("%-9.5f",a[i][j]);}printf("\n");printf("超松弛变量法得解:\n");w=1.4; /*超松弛法*/for(i=0;i<9;i++)x1[i]=0;for(i=0;i<9;i++)for(j=0;j<9;j++){if(i==j)b1[i][j]=0;else b1[i][j]=-a[i][j]/a[i][i];} /*求出矩阵b1[9][9]和b1[i][9]的值*/ for(i=0;i<9;i++)b1[i][9]=b[i]/a[i][i];for(n=0;n<9;n++)for(i=0;i<9;i++){z=0;for(j=0;j<9;j++)z=z+b1[i][j]*x1[j]; /*执行本算法*/z=z+b1[i][9];x1[i]=x1[i]*(1-w)+w*z;}for(i=0;i<9;i++){if (i==5)printf("\n");printf("x%d=%-10.6f",i,x1[i]);}printf("\n");printf("列主元消去法得解:\n");u[0]=a[0][0]; /*以下是消去法*/y[0]=b[0];for(i=1;i<9;i++){q[i]=a[i][i-1]/u[i-1];u[i]=a[i][i]-q[i]*a[i-1][i];y[i]=b[i]-q[i]*y[i-1];}x[8]=y[8]/u[8]; /*执行本算法*/for(i=7;i>=0;i--)x[i]=(y[i]-a[i][i+1]*x[i+1])/u[i]; /*求出x的值*/ for(i=0;i<9;i++){if (i==5)printf("\n");printf(" x%d=%-10.6f",i,x[i]);}printf("\n");}int sign(double x){int z;z=(x>=(1e-6)?1:-1);return(z);}（三）、计算结果打印：（四）、问题讨论：a ．由于选主元，使|lij|最小，这样去乘方程的每一系数时，系数的舍入误差不至扩大，并防止溢出与停机。

b ． 由于选主元，回代时作除数主元aii(i)的绝对值也是最大,这样扩大的误差也是最小.c ． 若detA ≠0,则选主元后的主元aii(i) ≠0,这是因为若aii(i)=0则必有aji(i)=0 (j>i)这样按行列式的Laplace 展开式就有detA=0而矛盾,因此用主元消去法中进行下去,不止中断停机.同时也由于选主元, aii(i)接近零的概率减少.运行结果基本与准确值无异，因为这种算法的无误差的算法。

三．试用三次样条插值求()4.563f 及(4.563)f '的近似值。

（一）方法由s(x i )=y i ,I=1,2,…N s’(x 0)=y 0’,s ’(x N )=y N ’可得S(x i )= ∑+-=11N j c j Ω3(x i -x j /h)=y iS ’(x 0)=1/h ∑+-=11N j c j Ω3’(x 0-x j /h)=y ’0S ’(x N )=1/h ∑+-=11N j c j Ω3’(x N -x j /h)=y ’N其中j c 必须的方程 ：⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫⎝⎛=⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎭⎫⎝⎛---2.081551.1318338.1347667.12675460.11750557.106566268.93177664.85916738.6158839.4021011411411411411411411411411411411011098765432101c c c c c c c c c c c cq[0]=0 , u[0]=0 ,1,,2,1]),[][][(][][-⋅⋅⋅=⋅+-=n i i q i a i b i c i qn i i q i a i b i u i a i d i u ,,2,1]),1[][][(])[][][(}[⋅⋅⋅=-⋅+⋅-= x[9]=u[9]1,,2,1],[]1[][][⋅⋅⋅--=++⋅=n n i i u i x i q i x []333333)2()1(46)1(4)2(61)(+++++-+--++-+=Ωx x x x x x s ' (x)= ∑∑-=-=--Ω--+Ω10121012)21()21(j j j j j x c j x c（二）原程序及运行结果：#include<stdio.h> #include<math.h> #define N 9 void main() {double x[N+1]={1,2,3,4,5,6,7,8,9,10},y[N+1]={0,0.69314718,1.0986123,1.3862944,1.6094378,1.7917595,1.9459101,2.079445,2.1 972246,2.3025851},h[N+1],d[N+1],a[N+1],c[N+1],b[N+1]={2,2,2,2,2,2,2,2,2,2},s[N+1],t[N+1],l[N+1],M[N+1], f,f1;int i;for(i=1;i<=N;i++) /*使用对任意分化的三弯矩插值法*/h[i-1]=x[i]-x[i-1];d[0]=6/h[0]*((y[1]-y[0])/h[0]-1);d[N]=6/h[N-1]*(0.1-(y[N]-y[N-1])/h[N-1]);for(i=1;i<=N-1;i++){d[i]=6/(h[i-1]+h[i])*((y[i+1]-y[i])/h[i]-(y[i]-y[i-1])/h[i-1]);a[i]=h[i-1]/(h[i-1]+h[i]);c[i]=1-a[i];}c[0]=1;a[N]=1;s[0]=b[0];t[0]=d[0]; /*用追赶法求解三对角方程组*/for(i=0;i<=N-1;i++){l[i+1]=a[i+1]/s[i];s[i+1]=b[i+1]-l[i+1]*c[i];t[i+1]=d[i+1]-l[i+1]*t[i];}M[N]=t[N]/s[N];for(i=N-1;i>=0;i--)M[i]=(t[i]-c[i]*M[i+1])/s[i];f=M[3]*pow((x[4]-4.563),3)/6/h[3] /*求算4.563这点的函数值*/+M[4]*pow((4.563-x[3]),3)/6/h[3]+(y[3]-M[3]*h[3]*h[3]/6)*(x[4]-4.563)/h[3]+(y[4]-M[4]*h[3]*h[3]/6)*(4.563-x[3])/h[3];f1=-3*M[3]*pow((x[4]-4.563),2)/6/h[3] /*求算4.563这点的一阶导数值*/+3*M[4]*pow((4.563-x[3]),2)/6/h[3] -(y[3]-M[3]*h[3]*h[3]/6)/h[3] +(y[4]-M[4]*h[3]*h[3]/6)/h[3];printf("f(4.563)=%lf f'(4.563)=%lf\n",f,f1);}（三） 输出结果：（四） 问题讨论：其基本思想是对均匀分划的插值函数的构造,三次样条函数空间中不取1,x,,x 2,x 3,(x-x j )+3为基函数,而取B 样条函数Ω3(x-x j /h)为基函数.由于三次样条函数空间是N+3维的,故我们把分点扩大到X -1,X N+1,则任意三次样条函数可用Ω3(x-x j /h)线性组合来表示 S(x)= ∑+-=11N j c j Ω3(x-x j /h) 这样对不同插值问题,若能确定c j 由解的唯一性就能求得解S(x).同时样条插值效果比Lagrange 插值好,近似误差较小.没有Runge 现象。