东南大学数值分析上机作业

习题1

17．（上机题）舍入误差与有效数 设S N =∑

j =2N

311⎫。 ，其精确值为1⎛ --⎪

2⎝2N N +1⎭j -1

2

1

111

，计算S N 的通用程序。 ++ +222

2-13-1N -1

11，计算S 的通用程序。 （2）编制按从小到大的顺序S N =1++ +N

N 2-1(N -1) 2-122-1

（1）编制按从大到小的顺序S N =

（3）按两种顺序分别计算S 102，S 104，S 106，并指出有效位数。（编制程序时用单精度） （4）通过本上机题，你明白了什么？

答：

1、从大到小顺序：

1. #include 2. int main() 3. {

4. float S=0,T,N,i=2; 5. printf(" 请输入N ：" ); 6. scanf("%f",&N); 7. while (i

9. T=1/(i*i-1); 10. i=++i; 11. S+=T; 12. }

13. printf("S=%f \n",S); 14. system("pause" ); 15. return 0; 16. }

2、从小到大顺序：

1. #include 2. int main() 3. {

4. float S=0,T,N; 5. printf(" 请输入N ：" ); 6. scanf("%f",&N); 7. while (N>=2) 8. {

9. T=1/(N*N-1); 10. N=--N; 11. S+=T; 12. }

13. printf("S=%f \n ",S); 14. system("pause" ); 15. return 0; 16. }

4、感想：

算法的不同导致计算结果的误差不同，一定要好好设计算法，尽量使算法精确。通过上面的例子，可以看出按不同的计算顺序所得到的结果是不同的，按从大到小的顺序计算的值

与精确值有较大的误差，而按从小到大的顺序计算的值与精确值吻合。计算机在进行数值计算时会出现“大数吃小数”的现象，导致计算结果的精度有所降低，使用计算机进行同号数的加法时，采用绝对值较小者先加的算法，得到的结果的相对误差较小。 5、附图：

图1从大到小算法及结果

图2从小到大算法及结果

习题2

20．（上机题）Newton 迭代法

（1）给定初值x 0及容许误差ε，编制Newton 法解方程f (x ) =0根的通用程序。

*（2）给定方程f (x ) =x /3-x =

0，易知其有三个根x 1=x 2=

0，x 3=

3**

1．由Newton 方法的局部收敛性可知存在δ>0，当x 0∈(-δ, δ) 时，Newton 迭代序

*

列收敛于根x 2。试确定尽可能大的δ。

2．试取若干初始值，观察当x 0∈(-∞, -1) ，(-1, -δ) ，(-δ, δ) ，(δ,1) ，(1,∞) 时Newton 序列是否收敛以及收敛于哪一个根。

（3）通过本上机题，你明白了什么？

答：

1、通用程序：

17. float x0,x1,d,eps=0.0000005;

1. #include 2. #include

3. float f(float x) 4. {

5. float f;

6. f=x*x*x/3-x; 7. return (f); 8. }

9. float df(float x) 10. {

11. float df; 12. df=x*x-1; 13. return (df); 14. }

15. int main() 16. {

18. int k=0;

19. printf(" 请输入初值：" ); 20. scanf("%f",&x0); 21. do 22. {

23. d=-f(x0)/df(x0); 24. x1=x0+d; 25. k++; 26. x0=x1; 27. }

28. while (fabs(d)>eps); 29. printf(" 迭代%d次后, 根

=%f\n",k,x0); 30. system("pause" ); 31. return 0; 32. }

2、确定尽可能大的δ：（由于函数是奇函数，故仅需在正半轴求δ）

4. {

1. #include 2. #include

3. float f(float x)

5. float f;

6. f=x*x*x/3-x; 7. return (f);

8. }

9. float df(float x) 10. {

11. float df; 12. df=x*x-1; 13. return (df); 14. }

15. int main() 16. {

17. float x0=0,x1,d;

18. float eps=0.000005,delta=0; 19. while (x0==0) 20. {

21. delta=delta+0.000001; 22. x0=x0+delta;

23. do 24. {

25. d=-f(x0)/df(x0); 26. x1=x0+d; 27. x0=x1; 28. }

29. while (fabs(d)>eps); 30. }

31. delta=delta-0.000001; 32. printf(" 尽可能大的值

=%f\n",delta); 33. system("pause" ); 34. return 0; 35. }

运算结果为：0.774596

3、测试是否在给定区间收敛：

在五个区间内分别取初值-100000，-0.8，0.5，0.9，110000测试是否收敛及收敛到哪个根，结果如下表：

4、感想：

选取的迭代初值不同，则迭代序列将有可能收敛于不同的根，牛顿法对迭代初值的要求较高。有些区间收敛范围相当大，有些却非常小。另外，我发现容许误差的大小对迭代结果的精度有很大的影响，如果想得到比较精确的结果，要将容许误差设置小一点。

5、附图：

图 1 利用牛顿法求方程的近似解

图 2 确定收敛范围

习题

3

39．（上机题）列主元Gauss 消去法

对于某电路的分析，归结为求解线性方程组RI=V，其中 R 为9阶矩阵，参见课本127页

VT=（−15，27，−23，0，−20，12，−7，7，0）

（1）编制解n 阶线性方程组Ax=b的列主元三角分解法的通用程序；

（2）用所编制的程序解线性方程组RI=V，并打印出解向量，保留五位有效数； （3）本编程之中，你提高了哪些编程能力？ 答:

1、通用程序

1. # include 2. # include 3. # define delta 1e-6 4. #define N 100 5. int main() 6. {

7. int i,j,t,r,n,u,c=0; 8. float p,L,max,s; 9. float X[N]; 10. float a[N][N+1];

11. printf(" 请输入系数矩阵的阶数：\n"); 12. scanf("%d",&n);

13. printf(" 输入的增广矩阵系数, 中间用空格隔开：\n"); 14. for (i=0;i

16. scanf("%f",&a[i][j]); 17. printf(" 你输入的增广矩阵为：\n"); 18. for (i=0;i

20. for (j=0;j

21. printf("%f\t",a[i][j]); 22. printf("\n"); 23. }

24. for (j=0;j

27. max=fabs(a[j][j]); 28. r=j; 29. }

30. for (i=j+1;imax) 32. {

33. max=a[i][j]; 34. r=i;

35. }

36. if (fabs(a[i][j])

40. p=a[j][t]; 41. a[j][t]=a[r][t]; 42. a[r][t]=p; 43. }

44. for (i=j+1;i

46. L=a[i][j]/a[j][j]; 47. for (t=j;t

48. a[i][t]=a[i][t]-L*a[j][t]; 49. } 50. }

51. printf(" 输出原方程的解为：\n"); 52. X[n-1]=a[n-1][n]/a[n-1][n-1]; 53. for (i=n-2;i>=0;i--) 54. {

55. s=a[i][n];

56. for (j=i+1;j

60. for (u=0;u

63. printf("%12f",a[u][j]); 64. c++;

65. if (c%(n+1)==0) 66. printf("\n"); 67. }

68. for (i=0;i

70. printf("x(%d)=%.4f\n",i+1,X[i]); 71. if (i==n-1) 72. printf("\n"); 73. }

74. system("pause" ); 75. return 0; 76. }

2、运行程序，方程的解为：

3、感想：

首先，通过本上机题，提高了我运用循环语句的能力；其次，是我对各种排序算法有了更深的理解；最后，是我对二维数组有了更好的把握，并且认识到MATLAB 在处理矩阵运算上的优势。在上机的过程中，加深了对Gauss 消去法的认识。 4、附图：

图 1 输入增广矩阵

习题4

37．（上机题）3次样条插值函数

（1）编制求第一型3次样条插值函数的通用程序；

端点条件为y010S(i+0.5)，i=0,1，„，9。 答：

（1）通用程序：

1. # include 2. # include 3. # include 4. # define delta 1e-6 5. # define N 100 6. int main() 7. {

8. int i,j,t,r,n,u,c=0; 9. float p,L,max,s,h,D1,D2; 10. float x[N],M[N],y[N]; 11. float a[N][N+1]={0}; 12. /*输入插值条件*/

13. printf(" 请输入插值节点的个数：\n"); 14. scanf("%d",&n);

15. printf(" 请从左到右输入插值节点：\n"); 16. for (i=0;i

18. scanf("%f",&x[i]); 19. }

20. printf(" 请输入插值节点对应的函数值：\n"); 21. for (i=0;i

23. scanf("%f",&y[i]); 24. }

25. printf(" 请输入插值初始条件：\n"); 26. printf("y'(0)="); 27. scanf("%f",&D1); 28. printf("y'(n)="); 29. scanf("%f",&D2); 30. /*产生弯矩方程组增广矩阵系数*/ 31. a[0][1]=1;

32. a[n-1][n-2]=1; 33. for (i=0;i

37. h=(x[i]-x[i-1])/(x[i+1]-x[i-1]); 38. a[i][i-1]=h; 39. a[i][i+1]=1-h; 40. }

41. a[0][n]=6*((y[n-2]-y[n-1])/(x[n-2]-x[n-1])-D2)/(x[n-2]-x[n-1]); 42. a[n-1][n]=6*((y[1]-y[0])/(x[1]-x[0])-D1)/(x[1]-x[0]); 43. for (i=1;i

45. a[i][n]=6*((y[i+1]-y[i])/(x[i+1]-x[i])-(y[i]-y[i-1])/(x[i]-x[i-1]))/

(x[i+1]-x[i-1]); 46. }

47. printf(" 请确认弯矩方程的增广矩阵为：\n"); 48. for (i=0;i

50. for (j=0;j

51. printf("%3.2f\t",a[i][j]); 52. printf("\n"); 53. }

54. for (j=0;j

57. max=fabs(a[j][j]); 58. r=j; 59. }

60. for (i=j+1;imax) 62. {

63. max=a[i][j]; 64. r=i; 65. }

66. if (max

67. printf(" 矩阵奇异\n"); 68. for (t=j;t

70. p=a[j][t]; 71. a[j][t]=a[r][t]; 72. a[r][t]=p; 73. }

74. for (i=j+1;i

75. {

76. L=a[i][j]/a[j][j];

77. for (t=j;t

78. a[i][t]=a[i][t]-L*a[j][t];

79. }

80. }

81. /*输出弯矩方程的计算结果*/

82. printf(" 整理后的增广矩阵为：\n");

83. M[n-1]=a[n-1][n]/a[n-1][n-1];

84. for (i=n-2;i>=0;i--)

85. {

86. s=a[i][n];

87. for (j=i+1;j

88. s=s-a[i][j]*M[j];

89. M[i]=s/a[i][i];

90. }

91. for (u=0;u

92. for (j=0;j

93. {

94. printf("%3.2f\t",a[u][j]);

95. c++;

96. if (c%(n+1)==0)

97. printf("\n");

98. }

99. printf(" 弯矩方程解向量为：\n");

100. for (i=0;i

101. {

102. printf("M(%d)=%f\n",i+1,M[i]);

103. if (i==n-1)

104. printf("\n");

105. }

106. /*计算插值函数*/

107. printf(" 样条函数如下：\n");

108. for (i=0;i

109. {

110. printf("S(x)=%f+%f(x-%.2f)+%f(x-%.2f)^2+%f(x-%.2f)^3\t%2.1f

" ,y[i],((y[i+1]-y[i])/(x[i+1]-x[i])-(M[i]/3+M[i+1]/6)*(x[i+1]-x[i])),x[i],M[i]/2,x[i],(M[i+1]-M[i])/6/(x[i+1]-x[i]),x[i],x[i],x[i+1]);

111. printf("\n");

112. }

113. printf("S(i+0.5)的值如下：\n");

114. for (i=0;i

115. {

116. printf("S(%d+0.5)=%f",i,y[i]+((y[i+1]-y[i])/(x[i+1]-x[i])-(M[i]/3+M

[i+1]/6)*(x[i+1]-x[i]))*(i+0.5-x[i])+M[i]/2*(i+0.5-x[i])*(i+0.5-x[i])+(M[i+1]-M[i])/6/(x[i+1]-x[i])*(i+0.5-x[i])*(i+0.5-x[i])*(i+0.5-x[i]));

117. printf("\n");

118. }

119. system("pause" );

120. return 0;

121. }

（2）打印出3次样条插值函数S(x),并打印出S(i+0.5)，（i=0，1，…，9）

1、样条函数如下：

2.510000+0.690000x +0.189038x 2−0.089039x 3 (0,1) 23 3.300000+0.800961(x−1) −0.078078(x−1) +0.017117(x−1) (1,2)

4.040000+0.696156(x−2) −0.026728(x−2) 2−0.009428(x−2) 3 (2,3) 4.700000+0.614416(x−3) −0.055011(x−3) 2−0.039405(x−3) 3 (3,4) 5.220000+0.386178(x−4) −0.173227(x−4) 2 +0.107049(x−4) 3 (4,5)S x =

23 5.540000+0.360870(x−5) +0.147919(x−5) −0.268789(x−5) (5,6)

5.780000−0.149659(x−6) −0.658448(x−6) 2 +0.428107(x−6) 3 (6,7) 5.400000−0.182234(x−7) +0.625873(x−7) 2−0.273639(x−7) 3 (7,8) 23 5.570000+0.248596(x−8) −0.195043(x−8) +0.076447(x−8) (8,9)

5.70000+0.087851(x−9) +0.034299(x−9) 2−0.022149(x−9) 3 (9,10)

（3）附图：

图 1 输入插值条件

图 2 生成弯矩方程的增广矩阵

图 3 输出插值函数及计算结果

习题6

23．（上机题）常微分方程初值问题数值解

常微分方程初值问题数值解

（1）编制RK4方法的通用程序；

（2）编制AB4方法的通用程序（由RK4提供初值）；

（3）编制AB4-AM4预测校正方法通用程序（由RK4提供初值）；

（4）编制带改进的AB4-AM4预测校正方法通用程序（由RK4提供初值）；

（5）对于初值问题

y ′=−x2y2 y 0 =3

取步长h=0.1，应用（1）~（4）中的四种方法进行计算，并将计算结果和精确解 y x =3/(1+x3) 作比较；

（6）通过本上机题，你能得到哪些结论？ 答：

（1）RK4通用程序： 1. #include

2. #include

3. #include

4. #define H 0.1

5. int main()

6. {

7. float x,y,r=3;

8. float k1,k2,k3,k4;

9. float fun1(float , float );

10. int i,n;

11. printf(" 请输入x,y 的初值：\n");

12. scanf("%f%f",&x,&y);

13. printf(" 请输入运算次数：\n");

14. scanf("%d",&n);

15. printf(" 运算结果为：\n");

16. for (i=0;i

17. {

18. printf("x%-2d=%.2f,\ty%-2d=%f,\ty(%-2d)=%f,\tr(%-2d)=%f",i,x,i,y,i,r

,i,r-y);

19. printf("\n");

20. k1=fun1(x,y);

21. k2=fun1(x+H/2,y+H*k1/2);

22. k3=fun1(x+H/2,y+H*k2/2);

23. k4=fun1(x+H,y+H*k3);

24. x=x+H;

25. y=y+H*(k1+2*k2+2*k3+k4)/6;

26. r=3/(1+x*x*x);

27. }

28. system("pause" );

29. return 0;

30. }

31. float fun1(float x,float y)

32. {

33. float f;

34. f=-x*x*y*y;

35. return f;

36. }

（2）AB4通用程序：

1. #include

2. #include

3. #include

4. #define H 0.1

5. int main()

6. {

7. float x,y,r=3;

8. float k1,k2,k3,k4,t[100]={0};

9. float fun1(float , float );

10. int i,j,n;

11. printf(" 请输入x,y 的初值：\n");

12. scanf("%f%f",&x,&y);

13. printf(" 请输入运算次数：\n");

14. scanf("%d",&n);

15. printf(" 运算结果为：\n");

16. for (i=4;i

17. {

18. if (i

19. {

20. for (j=0;j

21. {

22. printf("x%-2d=%.2f,\ty%-2d=%f,\ty(%-2d)=%f",j,x,j,y,j,r);

23. printf("\n");

24. k1=fun1(x,y);

25. k2=fun1(x+H/2,y+H*k1/2);

26. k3=fun1(x+H/2,y+H*k2/2);

27. k4=fun1(x+H,y+H*k3);

28. t[j]=y;

29. x=x+H;

30. y=y+H*(k1+2*k2+2*k3+k4)/6;

31. r=3/(1+x*x*x);

32. }

33. x=x-H;

34. }

35. t[i]=t[i-1]+H*(55*fun1(x,t[i-1])-59*fun1(x-H,t[i-2])+37*fun1(x-2*H,t

[i-3])-9*fun1(x-3*H,t[i-4]))/24;

36. x=x+H;

37. r=3/(1+x*x*x);

38. printf("x%-2d=%.2f,\ty%-2d=%f,\ty(%-2d)=%f,\tr(%-2d)=%f",i,x,i,t[i],

i,r,i,r-t[i]);

39. printf("\n");

40. }

41. system("pause" );

42. return 0;

43. }

44. float fun1(float x,float y)

45. {

46. float f;

47. f=-x*x*y*y;

48. return f;

49. }

（3）AB4-AM4通用程序：

1. #include

2. #include

3. #include

4. #define H 0.1

5. int main()

6. {

7. float x,y,r=3;

8. float k1,k2,k3,k4,t[100]={0},Y[100]={0};

9. float fun1(float , float );

10. int i,j,n;

11. printf(" 请输入x,y 的初值：\n");

12. scanf("%f%f",&x,&y);

13. printf(" 请输入运算次数：\n");

14. scanf("%d",&n);

15. printf(" 运算结果为：\n");

16. for (i=4;i

17. {

18. if (i

19. {

20. for (j=0;j

21. {

22. printf("x%-2d=%.2f,\ty%-2d=%f,\ty(%-2d)=%f",j,x,j,y,j,r);

23. printf("\n");

24. k1=fun1(x,y);

25. k2=fun1(x+H/2,y+H*k1/2);

26. k3=fun1(x+H/2,y+H*k2/2);

27. k4=fun1(x+H,y+H*k3);

28. Y[j]=y;

29. x=x+H;

30. y=y+H*(k1+2*k2+2*k3+k4)/6;

31. r=3/(1+x*x*x);

32. }

33. x=x-H;

34. }

35. t[i]=Y[i-1]+H*(55*fun1(x,Y[i-1])-59*fun1(x-H,Y[i-2])+37*fun1(x-2*H,Y

[i-3])-9*fun1(x-3*H,Y[i-4]))/24;

36. Y[i]=Y[i-1]+H*(9*fun1(x+H,t[i])+19*fun1(x,Y[i-1])-5*fun1(x-H,Y[i-2])

+fun1(x-2*H,Y[i-3]))/24;

37. x=x+H;

38. r=3/(1+x*x*x);

39. printf("x%-2d=%.2f,\ty%-2d=%f,\ty(%-2d)=%f,\tr(%-2d)=%f",i,x,i,Y[i],

i,r,i,r-Y[i]);

40. printf("\n");

41. }

42. system("pause" );

43. return 0;

44. }

45. float fun1(float x,float y)

46. {

47. float f;

48. f=-x*x*y*y;

49. return f;

50. }

（4）改进的AB4-AM4通用程序：

1. #include

2. #include

3. #include

4. #define H 0.1

5. int main()

6. {

7. float x,y,r=3;

8. float k1,k2,k3,k4,t[100]={0},Y[100]={0},R[100]={0};

9. float fun1(float , float );

10. int i,j,n;

11. printf(" 请输入x,y 的初值：\n");

12. scanf("%f%f",&x,&y);

13. printf(" 请输入运算次数：\n");

14. scanf("%d",&n);

15. printf(" 运算结果为：\n");

16. for (i=4;i

17. {

18. if (i

19. {

20. for (j=0;j

21. {

22. printf("x%-2d=%.2f,\ty%-2d=%f,\ty(%-2d)=%f",j,x,j,y,j,r);

23. printf("\n");

24. k1=fun1(x,y);

25. k2=fun1(x+H/2,y+H*k1/2);

26. k3=fun1(x+H/2,y+H*k2/2);

27. k4=fun1(x+H,y+H*k3);

28. R[j]=y;

29. x=x+H;

30. y=y+H*(k1+2*k2+2*k3+k4)/6;

31. r=3/(1+x*x*x);

32. }

33. x=x-H;

34. }

35. t[i]=R[i-1]+H*(55*fun1(x,R[i-1])-59*fun1(x-H,R[i-2])+37*fun1(x-2*H,R

[i-3])-9*fun1(x-3*H,R[i-4]))/24;

36. Y[i]=R[i-1]+H*(9*fun1(x+H,t[i])+19*fun1(x,R[i-1])-5*fun1(x-H,R[i-2])

+fun1(x-2*H,R[i-3]))/24;

37. R[i]=251.0/270*Y[i]+19.0/270*t[i];

38. x=x+H;

39. r=3/(1+x*x*x);

40. printf("x%-2d=%.2f,\ty%-2d=%f,\ty(%-2d)=%f,\tr(%-2d)=%f",i,x,i,R[i],

i,r,i,r-R[i]);

41. printf("\n");

42. }

43. system("pause" );

44. return 0;

45. }

46. float fun1(float x,float y)

47. {

48. float f;

49. f=-x*x*y*y;

50. return f;

51. }

（5）数据比较：

图 1 RK4解法

图 2

AB4解法

图 3 AB4-AM4解法

图 4 改进的AB4-AM4解法

（6）结论：

通过上述数据对比，可以发现，对于本题而言，RK4解法的精度最高，改进的AB4-AM4方法次之，AB4解法精度最低。