贵州大学数值分析上机实验答案
数值分析上机实验报告
课程名称：数值分析上机实验
学院：机械工程学院专业：机械制造
姓名：******* 学号：********** 年级：12级任课教师：***老师
2012年12月30日
一．已知A 与b
12.38412 2.115237 -1.061074 1.112336 -0.1135840.718719 1.742382 3.067813 -2.031743 2.11523719.141823 -3.125432 -1.012345 2.189736 1.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.113584
2.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.1034581
4.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 3
3.992318 -25.173417 0.84671695 1.784317 -86.612343 1.1101230
4.719345 -
5.6784392]T B ?????
???
?
?
????
?
?
????
?
?
?????
?
=（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;
double
b[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 )= ∑+-=1
1N j c j Ω3(x i -x j /h)=y i
S ’
(x 0)=1/h ∑+-=1
1
N j c j Ω3’(x 0-x j /h)=y ’0
S ’
(x N )=1/h ∑+-=1
1
N j c j Ω3’(x N -x j /h)=y ’N
其中j c 必须的方程 ：
???
??
?
??
?????
?
??
?
?
?
??=??????????????????? ?????????????????????
??---2.081551.1318338.1347667.12675460.11750557.106566268.93177664.85916738.6158839.4021011411411411411411
411411411411411011098765432101c c c c c c c c c c c c
q[0]=0 , u[0]=0 ,
1,,2,1]),[][][(][][-???=?+-=n i i q i a i b i c i q
n 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 []
3
33333)2()1(46)1(4)2(6
1)(+++++-+--++-+=
Ωx x x x x x s ' (x)= ∑∑-=-=--Ω--+Ω10
1
210
12)21
()21(j j j j j x c j x c
（二）原程序及运行结果：
#include
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)= ∑+-=1
1N j c j Ω3(x-x j /h) 这样对不同插值问题,
若能确定c j 由解的唯一性就能求得解S(x).同时样条插值效果比Lagrange 插值好,近似误差较小.没有Runge 现象。
四．用Newton 法求方程：x 7
-28x 4
+14=0 在(0.1,1.9)中的近似根(初始近似值取为区间端点,迭代6次或误差小于0.00001.)
解：（一） 理论依据：
（1）Newton 迭代法的几何意义为1k x +是()f x 在点k x 的切线与x 轴的交点，故也称为切线法，它是平方收敛的，但它是局部收敛的，即要求初始值与根α充分接近，对于任意初始近似值0x 的非局部收敛，有下面定理；
（2）Newton 迭代法定理：设函数在有限区间[],a b 上二阶导数存在，且满足条件 1、()()0f a f b <；
2、"()f x 在区间[],a b 上不变号；
3、'()0f x ≠；
4、
'()
(),f c f c b a
≤-其中c 是,a b 中使()
''min (),()f a f b 达到的一个。 则对任意初始近似值[]0,x a b ∈，由Newton 迭代过程1'()
(),()
k k k k k f x x x x f x ?+==-
0,1,2
k =
所生的迭代序列{}k x 平方收敛于方程()0f x =在区间[],a b 上的唯一解α。
（二）、计算程序（C 语言）：
#include
}
while(y>=0.00001); printf("x=%lf\n",x1);
}
}（三）、计算结果打印：
（四）、问题讨论：
当x=1.9时运行结果在[0.1，1.9]之间，是要求的解；当x=0.1时运行结果不在[0.1，1.9]之间，不符合要求。这并不是说当x=0.1时所得解不是原方程的解，而是所得解不在所要求的范围内。 五 。用Romberg 算法求
3
1.421
3(57)sin x x x x dx +?
（允许误差ε=0.000 01）
。 解：（一）、理论依据：
Romberg 算法的计算步骤如下：
（1）先求出按梯形公式所得的积分值：(0)
1
[()()]2
b a
T f a f b -=
+ （2）把区间2等分，求出两个小梯形面积之和，记为(
1)
1T ，即
(1)
1[()()2()]
42
b a b a
T f a f b f -+=
++ 这样由外推法可得(0)2T ,(1)2(0)
(1)(0)11(0)
112
21
()421411()2
T T T T T --==
--。 （3）把区间再等分（即22等分），得复化梯形公式(2)1T ，由(1)1T 与(2)1T 外推可得
(2)(1)(1)
112
441T T T
-=-，2(1)(0)(0)2232
441
T T T -=-，如此，若已算出2k 等分的复化梯形公式()1k T ，则由Richardson 外推法，构造新序列()(1)(1)
1
441
m k k
k m m m m T T T
--+-=
-， m=1,2,…,l, k=1,2,…,l-m+1,最后求得(0)
1l T +。 （4）(0)
(0)1l
l T T +≈或|(0)(0)
1
l l T T +-| <ε就停止计算，否则回到（3），计算(1)1l T +，一般可
用如下算法：
1(0)12()(1)1111()(1)(1)
1[()()]21{[(21)]}2224,1,2,,,1,2,,1
41l l l l l
i m k k k m m m m
b a
T f a f b b a
b a T T f a i T T T m l k l m ---=--+?-=+??
?--?=++-??
?-?===-+-??
∑
其具体流程如下：
(0)1(1)T
(0)2(3)T (0)3(6)T (1)1(2)T (1)2(5)T (2)1(4)T
（5）当ε<-+)
0(1
)
0(l l
T T 就停机
（二）、计算程序：
#include
f=pow(3.0,x)*pow(x,1.4)*(5*x+7)*sin(x*x); return (f);
} main() { int i=1,j,k,n=12;
float T[12],a=1.0,b=3.0,s=0.0; T[0]=0.5*(b-a)*(f(a)+f(b)); for(j=1;j { for(k=1;k<=pow(2,j-1);k++) s+=f(a+(2*k-1)*(b-a)/pow(2,j)); T[j]=0.5*(T[j-1]+(b-a)*s/pow(2,j-1)); s=0.0; } T[11]=(4*T[1]-T[0])/(float)3; for(;fabs(T[11]-T[0])>0.00001;i++) { T[0]=T[11]; for(j=1;j } printf("%f\n",T[11]); } （三）、计算结果打印： （四）、问题讨论： 程序中由count 计算循环次数，经过5次循环得到计算结果。 Romberg 算法的优点是： 1）把积分化为代数运算,而实际上只需求T1(i),以后用递推可得。 2）算法简单且收敛速度快,一般4或5次即能达到要求。 3）节省存储量，算出的(0) m T 可存入(0) 1 T 。 Romberg 算法的缺点是： 1）对函数的光滑性要求较高。 2）计算新分点的值时,这些数值的个数成倍增加 六 ．用定步长四阶Runge-Kutta 求解 ?????????? ?===--===0 )0(0)0(0)0(10010001000//1/321 3 2332 1y y y y y dt dy y dt dy dt dy h=0.0005,打印yi(0.025),yi(0.045),yi(0.085),yi(0.1),(i=1,2,3) 解：（一）、理论依据： 四阶Runge-Kutta 方法： 112341213243(22)611(,)2211 (,)22(,) n n n n n n n n n h y y k k k k k f k f t h y hk k f t h y hk k f t h y hk +?=++++?? =? ?? =++?? ?=++??=++?? （二）、计算程序（C 语言）： #include double F(double x,double y[4],double f[4]) { f[1]=0*x+0*y[1]+0*y[2]+0*y[3]+1; f[2]=0*x+0*y[1]+0*y[2]+1*y[3]+0; f[3]=0*x+0*y[1]-1000*y[2]-100*y[3]+1000; return(1); } void main() { double F(double x,double y[4],double f[4]); double h=0.0005,x=0,Y[4],k[5][4],s[4],f[4],det,m[4]={0.025,0.045,0.085,0.1}; int i,j,t; for(t=0;t<=3;t++) /*龙格—库塔算法*/ { for(j=0;j<=3;j++) Y[j]=0; //每求一组值后将初值条件还原为0 for(i=1;i<=int(m[t]/h);i++) { for(j=1;j<=3;j++) s[j]=Y[j]; det=F(x,s,f); for(j=1;j<=3;j++) k[1][j]=h*f[j]; /*四阶古典公式中的k值和求和的计算*/ for(j=1;j<=3;j++) s[j]=Y[j]+0.5*k[1][j]; det=F(x+0.5*h,s,f); for(j=1;j<=3;j++) k[2][j]=h*f[j]; for(j=1;j<=3;j++) s[j]=Y[j]+0.5*k[2][j]; det=F(x+0.5*h,s,f); for(j=1;j<=3;j++) k[3][j]=h*f[j]; for(j=1;j<=3;j++) s[j]=Y[j]+k[3][j]; det=F(x+h,s,f); for(j=1;j<=3;j++) k[4][j]=h*f[j]; for(j=1;j<=3;j++) Y[j]=Y[j]+(k[1][j]+2*k[2][j]+2*k[3][j]+k[4][j])/6; x+=h; } for(j=1;j<=3;j++) printf("y[%d](%f)=%f ",j,m[t],Y[j]); printf("\n"); } } （三）、结果打印： （四）、问题讨论： Runge_Kutta方法的优点： （1）精度高，不必用别的方法求开始几点的函数值。 （2）可根据f’(t,y)变化的情况与需要的精度自动修改步长。（3）程序简单，存储量少。 （4）方法稳定