数值分析编程
习题一：
已知A=
??
?
??
?
????
?
??
?
?
?
?-------------------------------00001.40446782.4213474.2238417.0784317.1037585.0010103.1123124.3031743.2446782.4719334.30123789.3103456.1121314.271828.0336993.0112384.1067813.3213474.2123789.37138465.14103458.0111223.3101023.2056781.1784165.0742382.1238417.0103456.1103458.0789365.9431637.0741856.3836742.1563849.1718719.0784317.1121314.2111223.3431637.0897918.19101011.4031454.2189736.2113584.0037585.071828.0101023.2741856.3101011.4108437.27123848.3012345.1112336.1010103.1336993.0056781.1741856.3031454.2123848.3567914.15125432.3061074.1123124.3112348.1784165.0836742.1189736.2012345.1125432.3141823.19115237.2031743.2067813.3742382.1718719.0113584.0112336.1061074.1115237.238412.12b=(2.1874369,33.992318,-25.173417,0.84671695,1.784317,-86.612343,1.1101230,4.719345,-5.6784392);
1． 用Household 变换，把A 化为三对角阵B 。（并打印B ）。 算法：
（1） )(,,,)
1(11)
1(0ij r r ij ij a A A a a A A ===--已知令。 （2） S r =（
21
1
2
)())(∑+=n r i r ir
a
（3） r r r r r r S a S )
(,1++=α
（4） T
r r n r r r r r r r r r r r a a s a sign a u ??
?
??????+???=+++)(,)(,2)(,1)(,1,,,)(,0,0,0,0
（5） r
r r r u
A y α1-=
（6） r r r r u k y q -=
（7）2,,2,1),(1-???=+-=-n r u q u q A A r T
r T
r r r r C 源程序:
#include
#include
#define N
main()
{ float a[N][N]=
{{12.38412,2.115237,-1.061074,1.112336,-0.113584,0.718719,1.7 42382,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.1010 23,-0.71828,-0.037585},
{-0.113584,2.189736,2.031454,4.101011,19.897918,0.431637,-3.1112 23,2.121314,1.784317},
{0.718719,1.563849,1.836742,-3.741856,0.431637,9.789365,-0.10 3458,-1.103456,0.238417},
{1.742382,-0.784165,-1.056781,2.101023,-3.111223,-0.103458,1 4.7138465,3.123789,-2.213474},
{3.067813,1.112348,0.336993,-0.71828,2.121314,-1.103456,3.123 789,30.719334,4.446782},
{-2.031743,3.123124,-1.010103,-0.037585,1.784317,0.238417,-2. 213474,4.446782,40.00001}
};
int i,j,r,n=N; /*分别定义变量*/
float m=0;
float s,p,k=0;
float u[N],y[N],q[N];
for(r=1;r<=n-2;r++) /*求S的值*/
{ for(i=r+1;i<=N;i++)
m=m+a[i-1][r-1]*a[i-1][r-1];
s=sqrt(m);
p=s*s+(fabs(a[r][r-1]))*s; /*求?的值,由P代替*/
for(i=n;i>=r+2;i--) /*U阵赋值*/
u[i-1]=a[i-1][r-1];
if(a[r][r-1]>0) u[r]=a[r][r-1]+s;
else u[r]=a[r][r-1]-s;
for(i=0;i<=N-1;i++) /*y阵赋值*/
{for(j=0;j<=N-1;j++)
y[i]=y[i]+a[i][j]*u[j];
y[i]=y[i]/p;
}
for(i=0;i<=n-1;i++) /*求K值*/
k+=y[i]*u[i];
k=0.5*k/p;
for(i=0;i<=n-1;i++)
u[i]=k*u[i];
for(i=0;i<=n-1;i++) /*求q阵*/
q[i]=y[i]-u[i];
for(i=0;i<=8;i++) /*求出新的变换阵a */
for(j=0;j<=N-1;j++)
a[i][j]=a[i][j]-(q[i]*u[j]+u[i]*q[j]);
}
cout<<"A"< for(i=0;i<=N-1;i++) for(j=0;j<=N-1;j++) { cout< if(j>=N-1) cout< } }运行结果: A= 12.384120 -4.893077 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 -4.893077 25.398416 6.494097 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 6.494097 20.611500 8.243927 -0.000000 -0.000000 -0.000000 0.000000 0.000000 0.000000 0.000000 8.243927 23.422845 -13.880073 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 -0.000000 -13.880073 29.698272 4.534501 0.000000 0.000000 0.000000 0.000000 0.000000 -0.000000 0.000000 4.534501 16.006116 4.881435 0.000000 -0.000000 0.000000 0.000000 -0.000000 0.000000 0.000000 4.881435 26.013316 -4.503634 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 -4.503634 21.254059 4.504498 0.000000 0.000000 0.000000 0.000000 0.000000 -0.000000 0.000000 4.504498 14.534124 问题总结: 在此过程中，由于开始对本题的算法理解不够透彻，出现了很多的错误，到后来仔细看书厚才明白其算法并实现了三对角阵的输出。感觉在实现三对角阵的输出的编程中实现U 阵是一个难点，我在此过程中费了好大力气才实现。 2. 用超松驰法求解BX=b (取松驰因子w=1.4,)0(X =0, 迭代9次)。算法: 1 () () (1) 1 1 i n m m m i ij j ij j i j j i x b x b x g --==+=+ +∑∑ (i=1,2,…,n) ， 当0ii a ≠时 ij ij ii a b a =- （i j ≠） ，i i g b = ，()ij B b = ，0ii b =，i,j=1,2,…n ； 当0ii a =时，ij ij b a =- （i j ≠） ，i i g b = ，()ij B b = ，1ii b = , i,j=1,2,…n 。 如果把已算出的()m i x 作为中间量并记为()m i x ，再用它与(1)m i x -作组合，得 ()()(1)(1)()(1)(1)()m m m m m m i i i i i i x x x x x x ωωω---=+-=+- ， 即 1()()(1)11 ()(1)()(1)()i n m m m i ij j ij j i j j i m m m m i i i i x b x b x g x x x x ω--==+--?=++???=+-?∑∑ 这里ω为松弛因子。当ω>1时称为超松弛法。 3. 用列主元素消去法求解BX=b . 算法: 采用列主元素消去法首先要求系数矩阵Ai 非奇异，并且在消元前通过行变换使得主对角线上的元素的绝对值比该列它以下的任何元素的绝对值都要大，接着再用Gauss 消元法。 第1列中第1行的元素与其下面的此列的元素逐一进行比较，找到最大的元素1j c ，将第j 行的元素与第1行的元素进行交换，然后通过行变换，将第1 列中第2到第n个元素都消成0。将变换后的矩阵(1) C的第二列中第二行的元素 与其下面的此列的元素逐一进行比较，找到最大的元素(1) 2 k c，将第k行的元素与第2行的元素进行交换，然后通过行变换，将第2列中第3到第n个元素都消成0。以此方法将矩阵的左下部分全都消成0。最终形式如下： （A b）～ ()() 1 111 () 00 n n n n n nn g a a g a ?? ? * ? ? ? ?? ][ ]1 [1i b i d= +, ] ][ 1 [ ]2 [1i i B i a+ = +, ] ][ [ ]1 [1i i B i b= +, ]1 ][ [ ]1 [1+ = +i i B i c,q1[0]=0 , u1[0]=0 , 8, ,2,1 ]), [1 ][1 ][1 ( ][1 ][1??? = ? + - =i i q i a i b i c i q 9, ,2,1 ]), 1 [1 ][1 ][1 ( ]) [1 ][1 ][1 ( } [1??? = - ? + ? - =i i q i a i b i u i a i d i u x[9]=u1[9]；1, ,7,8 ], [1 ]1 [ ][1 ][??? = + + ? =i i u i x i q i x C源程序: #include #include #define L 9 void main() { float a[L][L]= /*三对角阵a*/ {{12.384120,-4.893077,0,0,0,0,0,0,0}, {-4.893077,25.398416,6.494097,0,0,0,0,0,0}, {0,6.494097,20.611500,8.243927,0,0,0,0,0}, {0,0,8.243927,23.422845,-13.880073,0,0,0,0}, {0,0,0,-13.880073,29.698272,4.534501,0,0,0}, {0,0,0,0,4.534501,16.006116,4.881435,0,0}, {0,0,0,0,0,4.881435,26.013316,-4.503634,0}, {0,0,0,0,0,0,-4.503634,21.254059,4.504498}, {0,0,0,0,0,0,0,4.504498,14.534124}}; float b[L]={2.1874369,33.992318,-25.173417,0.84671695, 1.784317,-86.612643,1.1101230,4.719345,-5.6784392}; float x[L]={0}; /*存放超松弛法求解结果*/ float xx[L]={0}; /*存放消去法求解结果*/ float d[L+1][L+1]; /*增广矩阵*/ float w=1.4; float t=0,tt=0; int i,j; for(j=0;j<=L;j++) /*开始迭代*/ { x[0]=(1-w)*x[0]+w*(-a[0][1]/a[0][0]*x[1]+b[0]/a[0][0]); for(i=1;i<=7;i++) x[i]=(1-w)*x[i]+w*(-a[i][i-1]/a[i][i]*x[i-1]-a[i][i+1]/a[i][i]*x[i +1]+b[i]/a[i][i]); /*求出x[1]到[x[7]的值*/ x[8]=(1-w)*x[8]+w*(-a[8][7]/a[8][8]+b[8]/a[8][8]); } cout<<"松弛法求x的结果:"< for(i=0;i<=8;i++) /*打印x的结果*/ cout< for(i=0;i<=8;i++) /*写出增广矩阵*/ { for(j=0;j<=8;j++) d[i][j]=a[i][j]; d[i][L]=b[i]; } for(i=0;i<=7;i++) { if(fabs(d[i][i]) for(j=0;j<=L;j++) {t=d[i][j];d[i][j]=d[i+1][j];d[i+1][j]=t;} else { for(j=0;j<=9;j++) d[i+1][j]=-d[i+1][i]/d[i][i]*d[i][j]+d[i+1][j]; } } xx[8]=d[8][9]/d[8][8]; for(i=7;i>=0;i--) /*求出xx[7]到xx[0]的值*/ { for(j=i+1;j<=8;j++) tt=tt+d[i][j]*xx[j]; xx[i]=(d[i][9]-tt)/d[i][i]; tt=0; /*为tt清零以进入下一个x值的求取*/ } cout<<"列主元消去法求x的结果:"< for(i=0;i<=L-1;i++) /*输出出结果*/ cout< } 运行结果: 松弛法求解x的值(2) x[0]=1.075097 x[1]=2.276997 x[2]=-2.855843 x[3]=2.292943 x[4]=2.112816 x[5]=-6.425753 x[6]=1.362390 x[7]=0.659479 x[8]=-0.700549 列主元消去法解得x的值(3) xx[0]=1.075799 xx[1]=2.275743 xx[2]=-2.855515 xx[3]=2.293099 xx[4]=2.112634 xx[5]=-6.423838 xx[6]=1.357923 xx[7]=0.634244 xx[8]=-0.587266 问题总结: 在本题中由于是三对角阵所以当判断行与行是否交换时,只需要比较此行与基下一行的首个非零元素既可,可减少比较次数，这是做本题的一个要点。 习题三： 已知函数值如下表: x 1 2 3 4 5 f(x) 0 0.69314718 1.0986123 1.3862944 1.6094378 x 6 7 8 9 10 f(x) 1.7917595 1.9459101 2.079445 2.1972246 2.3025851 f’(x) f’(1)=1 f’(10)=0.1 试用三次样条插值法求f(4.563)及f’(4.563)的近似值. 算法： 当()1,j j x x x -∈时, 33 22 111 11 111 1 11 ()()()()()()() 6666j j j j j j j j j j j j j j j j x x x x h x x h x x S x M M y M y M h h h h ---------------=++-+-其中 11j j j h x x --=- ，只要求出j M ，则得三次样条插值函数S （x ）。得到j M 所满足的方程: 1000 00 1111 1 111 116()6()(1,2,,1)61[()] 1j j j j j j j j j N N N N N N j j j j j j j j j y y d y h h y y y y d j N h h h h d y y y h h h h h h h h μλμ+----------? '=-???--=-=-?+? ? ?'=--?? ?==-=?++? 最后将其归结为成求解一个三对角阵的解。 其三对角阵的Doolittle 分解为： 1 1112221222111111111n n n n n n n n n b c a b c l A LU l b c a l a b γβγββγβ-----????????????? ???????????===? ?????? ????????????????? 可将求解方程组Ax=d 的问题化为 1,,,n L d LUx d L d Ux δδδδδδ?? =???===???=????? 即若记则由得 112111n n n d l l d δδ??????????????????=???????????? ?????? ， 综上所诉,可得求解方程Ax=d 的算法(称为追赶法)： 111 1111111111,,,,1,2,3,,1,,1,,2,1 i i i i i i i i i i i n i i i n i n i b d l b l c d l i n c x x x i n αβδββδδδδββ+++++++++? ====-??? =-=-??-?===-?? 2． 程序算法 由s(xi)=yi,I=1,2,…N s ’(x0)=y0’,s ’(xN)=yN ’可得 S(xi)= cj Ω3(xi-xj/h)=yi S ’(x0)=1/hcj Ω3’(x0-xj/h)=y ’0 S ’(xN)=1/hcj Ω3’(xN-xj/h)=y ’N 其中必须的方程 ： ??? ?? ? ?? ????? ? ?? ? ? ? ??=??????????????????? ????????????????????? ??---01551.1481551.1318338.1347667.12675460.11750557.106566268.93177664.85916738.6158839.40224141141141141141141141141141141421098765432101c 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 3． 程序清单 #include #include //根据题目所得的矩阵方程AX=b的系数矩阵A; float a[L][L]={{-2,-4,0,0,0,0,0,0,0,0,0,0}, {1,4,1,0,0,0,0,0,0,0,0,0}, {0,1,4,1,0,0,0,0,0,0,0,0}, {0,0,1,4,1,0,0,0,0,0,0,0}, {0,0,0,1,4,1,0,0,0,0,0,0}, {0,0,0,0,1,4,1,0,0,0,0,0}, {0,0,0,0,0,1,4,1,0,0,0,0}, {0,0,0,0,0,0,1,4,1,0,0,0}, {0,0,0,0,0,0,0,1,4,1,0,0}, {0,0,0,0,0,0,0,0,1,4,1,0}, {0,0,0,0,0,0,0,0,0,1,4,1}, {0,0,0,0,0,0,0,0,0,0,4,2}}; float x[L],c[L],sum=0; //矩阵方程AX=b中的b矩阵; float y[L]={1,0,0.69314718,1.0986123,1.3862944,1.6094378,1.7917595,1.9459101,2.079445 ,2.1972246,2.3025851,0.1}; main() { int i; float h=1.0,b[L],e[L],l[L],f[L],d[L],k[L],d1[L],d2[L],k1[L],k2[L],t,r; x[1]=1;x[0]=x[1]-h; for(i=1;i x[i+1]=1+i*h; b[0]=2*h*y[0]; b[L-1]=2*h*y[L-1]; for(i=1;i b[i]=6*y[i]; b[0]=b[0]-b[1]; b[L-1]=b[L-1]+b[L-2]; e[0]=a[0][0];f[0]=b[0]; for(i=0;i { l[i+1]=a[i+1][i]/e[i]; e[i+1]=a[i+1][i+1]-l[i+1]*a[i][i+1]; f[i+1]=b[i+1]-l[i+1]*f[i]; } c[L-1]=f[L-1]/e[L-1]; for(i=L-2;i>=0;i--) c[i]=(f[i]-a[i][i+1]*c[i+1])/e[i]; for(i=0;i k[i]=4.563-x[i]; for(i=0;i { if(fabs(k[i])>=2) d[i]=0; else if(fabs(k[i])<=1) d[i]=(1/2.0)*pow(fabs(k[i]),3)-k[i]*k[i]+(2/3.0); else if((fabs(k[i])<2)&&(fabs(k[i])>1)) d[i]=(-1/6.0)*pow(fabs(k[i]),3)+k[i]*k[i]-2*fabs(k[i])+(4/3.0); } for(i=0;i sum=sum+c[i]*d[i]; cout<<"输出f(4.563)的值："< cout< for(i=0;i k2[i]=3.5-x[i]; for(i=0;i { if(fabs(k2[i])>=1.5) d2[i]=0; else if(fabs(k2[i])<=(1/2.0)) d2[i]=-(k2[i]*k2[i])+(3/4.0); else if((fabs(k2[i])<=(3/2.0))&&(fabs(k2[i])>=(1/2.0))) d2[i]=(1/2.0)*k2[i]*k2[i]-(3/2.0)*fabs(k2[i])+(9/8.0); } sum=0; for(i=0;i sum=sum+c[i]*d2[i]; for(i=0;i k1[i]=2.5-x[i]; for(i=0;i { if(fabs(k1[i])>=(3/2)) d1[i]=0; else if(fabs(k1[i])<=1/2.0) d1[i]=-(k1[i]*k1[i])+(3/4.0); else if((fabs(k1[i])<=(3/2.0))&&(fabs(k1[i])>=(1/2.0))) d1[i]=(1/2.0)*k1[i]*k1[i]-(3/2.0)*fabs(k1[i])+(9/8.0); } t=0; for(i=0;i t=t+c[i]*d1[i]; r=sum-t; cout<<"输出f(4.563)倒数的值："< cout< } 运行结果如下: Please enter the value of x=4.563 y=1.518241 y ’=0.207437 习题四： 用Newton 法求方程 0142847=+-x x 在(0.1,1.9)的近似根(初始值取为区间端点,迭代6次或误差小于0.00001) 算法： 由定理可得Newton 迭代法满足下列条件： i. ()()f a f b 〈 0； ii. '' f 在区间[a,b]上不变号； iii. ()' 0f x ≠； iv. ()()f c f c b a '≤-，其中c 是a,b 中使'' min (),()f a f b 达到的一个。 则对任意初始近似值[]0,x a b ∈，迭代方程为： 1() (),0,1,2() k k k k k f x x x x k f x ?+==- =' C 源程序: #include #include { double x; /*设定初始值,此处设为0.9*/ cout<<"请输入x 的值x="; cin>>x; double x0,x1; do { x0=x; x=x0-(pow(x0,7)-28*pow(x0,4)+14)/(7*pow(x0,6)-112*pow(x0,3)); x1=x; } while((abs(x1-x0))>=0.00001); /*循环条件*/ cout< 用x0、x1为中间过度变量。 习题五： 用Romberg 算法求3 1.4213(57)sin x x x x dx +?（允许误差ε=0.000 01）。 算法： 1(0)12()(1)1111()(1)(1) 1[()()]21{[(21)]}2224,1,2,,,1,2,,141l 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[()()]2 b a T f a f b -= +，再把区间2等分，求出两个小梯形面积之和，记为(1)1T ，即(1)1[()()2()]42 b a b a T f a f b f -+=++ 由外推法便得(1)2(0) (1)(0)11(0) 11221 ()421411()2 T T T T T --== --，然后把区间再等分（即22等分），得复化梯形公式(2)1 T ，并由(1)1T 与(2)1 T 外推可得(2)(1)(1)112 441 T T T -=-， 2(1)(0)(0)223 2 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 +。当(0)(0)1l l T T +≈或(0)(0)1||l l T T +-<ε就停止计算，否则回到(1) 1l T +. 程序流程: (0)1(1)T (0)2(3)T (0)3(6)T (1)1(2)T (1)2(5)T (2)1(4)T C 源程序: #include # define Precision 0.00001 # define e 2.71828183 #define M 10 double function(double x)// { double s; s=pow(3,x)*pow(x,1.4)*(5*x+7)*sin(pow(x,2)); return s; } double Romberg(double a,double b,double f(double x)) { int m,n,k; double y[M],h,ep,p,xk,s,q; h=b-a; y[0]=h*(f(a)+f(b))/2.0; m=1; n=1; ep=Precision+1; while((ep>=Precision)&&(m { p=0.0; for(k=0;k { xk=a+(k+0.5)*h; p=p+f(xk); } p=(y[0]+h*p)/2.0; s=1.0; for(k=1;k<=m;k++) { s=4.0*s; q=(s*p-y[k-1])/(s-1.0); y[k-1]=p; p=q; } ep=fabs(q-y[m-1]); m=m+1; y[m-1]=q; n=n+n; h=h/2.0; } return q; } main() { double a,b,Result; cout<<"a=:"< cin>>a; cout<<"b=:"< cin>>b; Result=Romberg( a, b, function); cout<<"龙贝格积分结果："< return 0; } }运行结果如下: y=440.536011 问题总结: 本题在做的时候一定要注意r 的清零，否则结果会累积，导致运行结果出错。 习题六： 用定步长四阶Runge-Kutta 法求解 12 3323 1 23/1//10001000100(0)0(0)0(0)0 dy dt dy dt y dy dt y y y y y =??=??=--?? =??=?=?? h=0.000 5，打印(0.025),(0.045),(0.085),(0.1),(1,2,3)i i i i y y y y i =。 算法： 将n t 处附近的值的线形组合（有待定常数）再按Taylor 公式展开，然后确定待定常数，便是Runge-Kutta 法的思想方法。 112341213 2431(22)6(,)11(,) 2211(,)22 (,) n n n n n n n n n n Y Y K K K K K hF x Y K hF x h Y K K hF x h Y K K hF x h Y K +?=++++?? =??? =++?? ? =++??=++?? C 源程序： #include { int i,k,l,h=0.0005; double x; cout<<"输入x="; cin>>x; (int)l=x/h; double y[4][200],k1[4][200],k2[4][200],k3[4][200]; y[1][0]=y[2][0]=y[3][0]=0; for(i=0;i<=l;i++) { k1[1][i]=1,k1[2][i]=1,k1[3][i]=1,k1[4][i]=1; y[1][i+1]=y[1][i]+h*(k1[1][i]+2*k1[2][i]+2*k1[3][i]+k1[4][i]) /6; k2[1][i]=y[3][i],k2[2][i]=y[3][i],k2[3][i]=y[3][i],k2[4][i]=y[3][i]; y[2][i+1]=y[2][i]+h*(k2[1][i]+2*k2[2][i]+2*k2[3][i]+k2[4][i])/6; k3[1][i]=1000-1000*y[2][i]-100*y[3][i]; k3[2][i]=1000-1000*y[2][i]-100*(y[3][i]+0.5*h*k3[1][i]);