西安交大2014级计算方法b大作业.doc

计算方法B上机作业

目录
第二章习题 1
题目2.4 1
题目分析 1
算法组织 2
结果分析 2
题目2.5 2
题目分析 2
算法组织 3
结果分析 4
第三章习题 5
题目 5
题目分析 5
算法组织 6
结果分析 6
第四章习题 6
题目 6
题目分析 6
算法组织 6
结果分析 6
第五章题目 6
题目 6
题目分析 6
算法组织 6
结果分析 6
第六章习题 6
题目 6
算法组织： 6
结果分析： 6
附录 6

第二章习题
题目2.4
求20阶三对角方程式Tx=f的解，其中
T= f=

题目分析
T=

T=LU
L=, U=
b1=u1,c1=r1
ak=lkuk-1,bk=lkrk-1+uk,ck=rk,k=2,3,...,n
由此可得计算公式
u1=b1,r1=c1
Lk=ak/uk-1,uk=bk-lkrk-1,rk=ck,k=2,3,...,n
由此可解三对角方程组
Tx=d
算法组织
采用追赶法
1. b1→u1，d1→y1
2. for k=2,3,....,n
2.1 ak/uk-1→lk,
2.2 bk-lk*ck-1→uk
2.3 dk-lk*yk-1→yk.
3. yn/un→xn
4. for k=n-1,n-2,....,1
4.1 (yk-ck*xk+1)/uk→xk
结果分析
x=[1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000]T;
理论解x=[1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1]T;迭代满足要求。
题目2.5
用Jacobi迭代和Gauss-Seidel迭代解电路电流方程组，使各部分电流的误差均小于10-3

题目分析
Jacobi迭代
矩阵A可分为三部分写成:A=D-E-F,其中D,-E,-F分别是A的对角线元、下三角部分元和上三角部分元所构成的矩阵。
（D-E-F）X=b可得
X=D-1(E+F)x+b
Gauss-Seidel迭代
矩阵A可分为三部分写成:A=D-E-F,其中D,-E,-F分别是A的对角线元、下三角部分元和上三角部分元所构成的矩阵。
（D-E）x=Fx+b可得
x=(D-E)-1Fx+(D-E)-1b
算法组织
Jacobi迭代
1.1 for i=1,2...,n
1.1.1()
1.2 Max
1.3 If e< then stop
1.4 for i=1,2,...,n
1.4.1
2. 输出错误信息：“不收敛”；stop
Gauss-Seidel迭代
1. for i=1,2,…,n
1.1
2. for 1,2....
2.1 for i=1,2...
2.1.1 ()
2.2 If max then stop
2.3 for i=1,2...
2.3.1
3. 输出错误信息：“不收敛”；stop
结果分析
Jacobi迭代
k=1，x=0.3571，0.0000，0.0000，0.0000，0.0000
k=2，x=0.3571，0.0282，0.0000，0.0000，0.0000
k=3，x=0.3602，0.0282，0.0113，0.0000，0.0047
k=4，x=0.3602，0.0320，0.0113，0.0038，0.0047
k=5，x=0.3606，0.0320，0.0151，0.0038，0.0053
k=6，x=0.3606，0.0331，0.0151，0.0050，0.0053
k=7，x=0.3607，0.0331，0.0163，0.0050，0.0055
k=8，收敛，x=0.3607，0.0335，0.0163，0.0054，0.0055
Gauss-Seidel迭代
k=1，x=0.3571，0.0282，0.0113，0.0038，0.0047
k=2，x=0.3602，0.0320，0.0151，0.0050，0.0053
k=3，x=0.3606，0.0331，0.0163，0.0054，0.0055
k=4，x=0.3607，0.0335，0.0166，0.0055，0.0056
经过计算，Jacobi迭代需要8次迭代才能得到结果，而Gauss-Seidel迭代仅需要四次迭代就能得到结果。这说明Gauss-Seidel迭代的收敛速度大于Jacobi迭代。
第三章习题
题目
Runge函数，，对n=5,10,20
1.计算在点，处的值；
2.求插值数据点，的Newton插值多项式N(x)和样条插值多项式s(x);
3.对n=5,20计算相应的函数值，和.
4.计算，.
题目分析
牛顿插值多项式：

其中，

计算出插值多项式中的各项系数，从而求出牛顿插值多项式。
样条插值多项式：
式中，

若用差商记号，并记

令，即有

采用自然样条插值，即有，代入解三对角方程组，即可求得
利用求得的牛顿插值多项式和样条差值多项式，分别进行计算，然后进行结果比较。
算法组织
1. 求得时各点的函数值，构造插值数据表。
2. 利用算法DD，求得插值多项式中的系数，，。
2.1. 算法DD（x,y,n）
2.1.1. for k=1,2,...,n
2.1.2 for i=n, n-1,...,k
2.1.2.1 (yi-yi-1)/(xi-xi-k)yi
3. 利用算法SPLINEM，解出系数，从而求得样条分段插值多项式。
3.1. for i=0,1,2,.n
1.1 yiMi
3.2. for k=1,2
3.2.1 for i=n, n-1,..,k
3.2.1.1 (Mi-Mi-1)/(xi-xi-1)Mi
3.3. x1-x0h1
3.4. for i=1,2,...,n-1
3.4.1 xi+1-xihi+1
3.4.2 hi+1/(hi+hi+1)ci;1-ciai;2bi
3.4.3 6Mi+1Mi
3.5. d0M0;dnMn;c0;
3.6. 解出Mi
4. 已知，，求得。利用算法NI求得，利用算法FINDK确定点所在区间，利用算法EVASPLINE求得。
4.1 算法NI
4.4.1.yny
4.4.2.for i=n-1,n-2,...1,0
4.4.2.1 y(x-xi)+yiy
4.2算法FINDK（x,y,n,,k），可找出所在区间
4.2.1. 1k
4.2.2. for i=1,2...,n-1
4.2.2.1 If xi then ik; goto 3
Else i+1k
4.2.3. 结束
4.3算法EVASPLINE（x,y,n,,,d0,un,dn）
4.3.1. xk-xk-1h; xk-; ;

5. 将所求数值进行比较，求得，
结果分析
1. 求得各点的函数值
n=5时，
插值节点xi
xi = -1.0000 -0.6000 -0.2000 0.2000 0.6000 1.0000
插值节点xi处的函数值f(xi)
fx = 0.0385 0.1000 0.5000 0.5000 0.1000 0.0385
n=10时，
插值节点xi
xi =-1.0000 -0.8000 -0.6000 -0.4000 -0.2000 0 0.2000 0.4000 0.6000 0.8000 1.0000
插值节点xi处的函数值f(xi)
fx = 0.0385 0.0588 0.1000 0.2000 0.5000 1.0000 0.5000 0.2000 0.1000 0.0588 0.0385
n=20时，
插值节点xi
xi =-1.0000 -0.9000 -0.8000 -0.7000 -0.6000 -0.5000 -0.4000 -0.3000 -0.2000 -0.1000 0 0.1000 0.2000 0.3000 0.4000
0.5000 0.6000 0.7000 0.8000 0.9000 1.0000
插值节点xi处的函数值f(xi)
fx = 0.0385 0.0471 0.0588 0.0755 0.1000 0.1379 0.2000 0.3077 0.5000 0.8000 1.0000 0.8000 0.5000 0.3077 0.2000 0.1379 0.1000 0.0755 0.0588 0.0471 0.0385
2. 求插值数据点的牛顿插值多项式和样条插值多项式
详细表达式可见附录
3. 求得个点相对应的函数值
当n=5时，
=0.0400 0.0416 0.0433 0.0451 0.0471 0.0491 0.0513 0.0536 0.0561 0.0588 0.0617 0.0648 0.0681 0.0716
0.0755 0.0796 0.0841 0.0890 0.0943 0.1000 0.1063 0.1131 0.1206 0.1289 0.1379 0.1479 0.1590 0.1712 0.1848 0.2000 0.2169 0.2358 0.2571 0.2809 0.3077 0.3378 0.3717 0.4098 0.4525 0.5000 0.5525 0.6098 0.6711 0.7353 0.8000 0.8621 0.9174 0.9615 0.9901 1.0000 0.9901 0.9615 0.9174 0.8621 0.8000 0.7353 0.6711 0.6098 0.5525 0.5000 0.4525 0.4098 0.3717 0.3378 0.3077 0.2809 0.2571 0.2358 0.2169 0.2000 0.1848 0.1712 0.1590 0.1479 0.1379 0.1289 0.1206 0.1131 0.1063 0.1000 0.0943 0.0890 0.0841 0.0796 0.0755 0.0716 0.0681 0.0648 0.0617 0.0588 0.0561 0.0536 0.0513 0.0491 0.0471 0.0451 0.0433 0.0416
0.0400
=0.0137 -0.0069 -0.0236 -0.0366 -0.0460 -0.0522 -0.0553 -0.0555 -0.0530 0.0481 -0.0408 -0.0314 -0.0200 -0.0069
0.0078 0.0240 0.0414 0.0600 0.0796 0.1000 0.1211 0.1427 0.1648 0.1872 0.2097 0.2323 0.2549 0.2773
0.2994 0.3212 0.3424 0.3632 0.3833 0.4027 0.4213 0.4390 0.4558 0.4716 0.4864 0.5000 0.5125 0.5238
0.5338 0.5426 0.5501 0.5563 0.5611 0.5645 0.5666 0.5673 0.5666 0.5645 0.5611 0.5563 0.5501 0.5426
0.5338 0.5238 0.5125 0.5000 0.4864 0.4716 0.4558 0.4390 0.4213 0.4027 0.3833 0.3632 0.3424 0.3212
0.2994 0.2773 0.2549 0.2323 0.2097 0.1872 0.1648 0.1427 0.1211 0.1000 0.0796 0.0600 0.0414 0.0240
0.0078 -0.0069 -0.0200 -0.0314 -0.0408 -0.0481 -0.0530 -0.0555 -0.0553 -0.0522 -0.0460 -0.0366 -0.0236 -0.0069
0.0137
=0.0360 00337 0.0315 0.0296 0.0280 0.0269 0.0262 0.0261 0.0266 0.0279 0.0301 0.0331 0.0371 0.0422
0.0485 0.0560 0.0648 0.0750 0.0867 0.1000 0.1149 0.1313 0.1490 0.1679 0.1878 0.2085 0.2300 0.2520
0.2744 0.2970 0.3196 0.3422 0.3645 0.3864 0.4077 0.4282 0.4479 0.4665 0.4840 0.5000 0.5145 0.5275
0.5390 0.5490 0.5574 0.5643 0.5696 0.5735 0.5758 0.5765 0.5758 0.5735 0.5696 0.5643 0.5574 0.5490
0.5390 0.5275 0.5145 0.5000 0.4840 0.4665 0.4479 0.4282 0.4077 0.3864 0.3645 0.3422 0.3196 0.2970
0.2744 0.2520 0.2300 0.2085 0.1878 0.1679 0.1490 0.1313 0.1149 0.1000 0.0867 0.0750 0.0648 0.0560
0.0485 0.0422 0.0371 0.0331 0.0301 0.0279 0.0266 0.0261 0.0262 0.0269 0.0280 0.0296 0.0315 0.0337
0.0360
当n=20时，
=- -58.2381 -50.8644 -28.6626 -10.3345 0.0471 3.9451 4.0691 2.6744 1.1353 0.0588 -0.4459 -0.5136 -0.3459 -0.1143
0.0755 0.1823 0.2101 0.1858 0.1408 0.1000 0.0770 0.0750 00896 0.1130 0.1379 0.1590 0.1742 0.1843
0.1918 0.2000 0.2115 0.2281 0.2501 0.2770 0.3077 0.3411 0.3765 0.4143 0.4550 0.5000 0.5503 0.6065
0.6682 0.7337 0.8000 0.8632 0.9189 0.9626 0.9904 1.0000 0.9904 0.9626 0.9189 0.8632 0.8000 0.7337
0.6682 0.6065 0.5503 0.5000 0.4550 0.4143 0.3765 0.3411 0.3077 0.2770 0.2501 0.2281 0.2115 0.2000
0.1918 0.1843 0.1742 0.1590 0.1379 0.1130 0.0896 0.0750 0.0770 0.1000 0.1408 0.1858 0.2101 01823
0.0755 -0.1143 -0.3459 -0.5136 -0.4459 0.0588 1.1353 2.6744 4.0691 3.9451 0.0471 -10.3345 -28.6626 -50.8644
-58.2381
Sk =0.0401 0.0417 0.0434 0.0452 0.0471 0.0491 0.0513 0.0536 0.0561 0.0588 0.0617 0.0648 0.0681 0.0716
0.0755 0.0796 0.0841 0.0890 0.0942 0.1000 0.1063 0.1131 0.1206 0.1289 0.1379 0.1479 0.1589 0.1712
0.1848 0.2000 0.2170 0.2359 0.2572 0.2810 0.3077 0.3375 0.3711 0.4089 0.4517 0.5000 0.5541 0.6127
0.6742 0.7372 0.8000 0.8606 0.9152 0.9594 0.9891 1.0000 0.9891 0.9594 0.9152 0.8606 0.8000 0.7372
0.6742 0.6127 0.5541 0.5000 0.4517 0.4089 0.3711 0.3375 0.3077 0.2810 0.2572 0.2359 0.2170 0.2000
0.1848 0.1712 0.1589 0.1479 0.1379 0.1289 0.1206 0.1131 0.1063 0.1000 0.0942 0.0890 0.0841 0.0796
0.0755 0.0716 0.0681 0.0648 0.0617 0.0588 0.0561 0.0536 0.0513 0.0491 0.0471 0.0452 0.0434 0.0417
0.0401
4. 当n=5时，=0.4327，=0.4235

当n=10时，= 1.9156，= 0.0220

当n=20时，=58.2781，=0.0031

由以上结果可知，随着插值数据点的数目增大，牛顿插值公式的误差增大，特别是两端的误差。这告诉我们，用高次插值多项式可能产生极大的误差，因此在实际应用中，要采用相对较小的n值。
而样条插值的误差随着数据点的增加而减小，可以避免大的误差，同时插值曲线比较光滑。
第四章习题
题目
用n=4的复化Simpson公式求积分方程

的近似解.然后用三次样条插值函数S(x)近似y(x),并对,(k=0,1,…,20,h=1/20),比较S(x)和真解在这些点上的值.
题目分析
令

即有

记

令

即，

即，
解方程可得Y。
算法组织
1. 构造矩阵
2. for i=,……9 （构造系数矩阵）
2.1 for j=1,2,……9
2.1.1 if j=1,9
2.1.2 if j=3.5.7，
2.1.3 if j=2,4,8,
3. 调用消元法程序解方程 ,得到Y。
3.1 for k=1,3,...,n-1
3.1.1 找满足的下标uk
3.1.2 If =0 then输出错误信息；stop
3.1.3 for j=1,2,..n
3.1.3.1
3.1.4
3.1.5 for i=k+1,k+2,..,n
3.1.5.1
3.1.5.2 for j=k+1,k+2,..n
3.1.5.2.1
3.1.5.3
3.2.[回带过程]
3.2.1
3.2.2 for k=n-1, n-2,...1
3.2.2.1
4. 调用样条差值函数近似求得S(x)
4.1. for i=0,1,2,..,n
1.1 yiMi
4.2. for k=1,2

4.2.1 for i=n, n-1,..,k
4.2.1.1 (Mi-Mi-1)/(xi-xi-1)Mi
4.3. x1-x0h1
4.4. for i=1,2,...,n-1
4.4.1 xi+1-xihi+1
4.4.2 hi+1/(hi+hi+1)ci;1-ciai;2bi
4.4.3 6Mi+1Mi
4.5. d0M0;dnMn;c0;
4.6. 解出Mi
4.7算法EVASPLINE（x,y,n,,,d0, un, dn）求得S(x)
4.7.1. 调用FINDK
4.7.2. xk-xk-1h; xk-; ;

5. 比较S(x)与y(x)在点上的值。
结果分析
求得Y= [1.0001，0.7902，0.6401，0.5290，0.4445，0.3788，0.3266，0.2845，0.2500]T
采用样条插值得
S(x)= 6.98*x^3 - 1.79*x + 1.0 0.00<=x<=0.13
4.24*x^2 - 2.32*x - 4.32*x^3 + 1.02 0.13<=x<=0.25
1.23*x^2 - 1.57*x - 0.305*x^3 + 0.959 0.25<=x<=0.38
1.83*x^2 - 1.79*x - 0.839*x^3 + 0.988 0.38<=x<=0.50
1.11*x^2 - 1.43*x - 0.363*x^3 + 0.92 0.50<=x<=0.63
1.09*x^2 - 1.42*x - 0.349*x^3 + 0.925 0.63<=x<=0.75
0.372*x^2 - 0.88*x - 0.0309*x^3 0.791 0.75<=x<=0.88
2.33*x^2 - 2.59*x - 0.776*x^3 + 1.29 0.88<=x<=1.00
将点代入，得到误差图像

求得其最大相对误差为0.0160
第五章题目
题目
5. 用Newton迭代法求下列方程组的解
a. (取)
b. (取)
题目分析
本题使用newton迭代法求非线性方程组解。
设在所讨论的区域内关于自变量各分量的二阶偏导数连续，利用多元函数的Taylor公式，在处展开并略去二次以上的项，可得近似公式

从而得线性方程组

得到

算法组织
1. 建立函数F（X）
2. 建立函数的一阶导数F’(X)
3. 算法NEWTON2
3.1 for k=0,1,2……，Nmax
3.1.1
3.1.2 if then
3.1.2.1 输出“已足够小”
3.1.2.2 ；STOP
3.1.3
3.1.4 解方程组
3.1.5 if then
3.1.5.1 输出“已足够小”；stop
3.1.6
3.2 输出“在Nmax次迭代后，仍不收敛”；stop
结果分析
a 由程序得
k =0，x =1.2392，y =16720
k =1，x =1.2338，y =1.6603
k =2，x =1.2343，y =1.6616
k =3，||x-x0||已足够小，x =1.2343，y =1.6616
经过四步迭代收敛
b 由程序得
k=0，x=1.0386，y=0.4717
k=1，||x-x0||已足够小，x=1.0386，y=0.4717
经过两步收敛。
第六章习题
题目
用Runge-Kutta方法(6-54)计算初值问题.

取步长h=0.25
算法组织：
采用四阶（经典）龙格-库塔公式：

结果分析：
t
h=0.25
0.25
0.2398

0.50
0.4240
0.75
0.5174
1.0
0.5068

附录
第三章习题牛顿插值多项式和样条插值多项式的表达式
当n=5时，
N(x)=1/26-(x+1)*(((((33*x)/9007199254740992-703687441776641287/585467951558164480)*(x-1/5) + 25/13)*(x + 1/5)-5/52)*(x + 3/5) - 2/13)
样条插值的表达式：在第1个区间(-1.0000<=x<=-0.6000)内，S(x)的表达式为：
y =x/76 + (425*(x + 1)^3)/247 + 51/988
在第2个区间(-0.6000<=x<=-0.2000)内，S(x)的表达式为：
y = (131*x)/152 - (3855*x^2)/988 - (3275*x^3)/988 + 707/520
在第3个区间(-0.2000<=x<=0.2000)内，S(x)的表达式为：
y = - (775*x)/104 - (945*x^2)/494 - 3519/9880
在第4个区间(0.2000<=x<=0.6000)内，S(x)的表达式为：
y = (3275*x^3)/988 (3855*x^2)/988 + (36399*x)/3952 - 3781/1040

在第5个区间(0.6000<=x<=1.0000)内，S(x)的表达式为：
y =12145082295111409541/4449556431842050048 - (425*(x - 1)^3)/247 - (12145082295111425843*x)/4449556431842050048
n=10时，
N(x)= 1/26 - ((((x + 2/5)*((x*(((x - 2/5)*(((7773696463348413*x)/35184372088832 - 69963268170135707/175921860444160)*(x - 3/5) + 5908009312144791/35184372088832) - 18750/221)*(x - 1/5) + 15625/884) + 5625/884)*(x + 1/5) - 2375/884) - 175/221)*(x + 3/5) - 115/442)*(x + 4/5) - 45/442)*(x + 1)
样条插值的表达式：在第1个区间(-1.0000<=x<=-0.8000)内，S(x)的表达式为：
y =(141951938076098285*x)/468374361246531584 (54675*(x + 1)^3)/160004 + 159966336585580269/468374361246531584
在第2个区间(-0.8000<=x<=-0.6000)内，S(x)的表达式为：
y =(27242268629192565*x)/36028797018963968 - (54675*(x + 3/5)^3)/160004 + (3800*(x + 4/5)^3)/3077 + 25363413836023409/36028797018963968
在第3个区间(-0.6000<=x<=-0.4000)内，S(x)的表达式为：
y = (145499871745901944681*x)/55430304213676064768 - (3800*(x + 2/5)^3)/3077 + (2332137800463225*(x + 3/5)^3)/1125899906842624 + 252935492656638812009/138575760534190161920
在第4个区间(-0.4000<=x<=-0.2000)内，S(x)的表达式为：
y = (15096357045089445*x^3)/1125899906842624 + (19514911134385269*x^2)/1125899906842624 + (558436177390937502207*x)/34643940133547540480 + 115287874206230894431/21652462583467212800
在第5个区间(-0.2000<=x<=0.0000)内，S(x)的表达式为：
y =- (537604469093067959*x)/11258999068426240 - (25721847370083*x^2)/1099511627776 - (30663557178580495*x^3)/562949953421312 - 179105966515547739/28147497671065600
在第6个区间(0.0000<=x<=0.2000)内，S(x)的表达式为：
y =(30663557178580505*x^3)/562949953421312 - (25721847370083*x^2)/1099511627776 - (152168429151608429*x)/2814749767106560 - 89552983257773857/14073748835532800
在第7个区间(0.2000<=x<=0.4000)内，S(x)的表达式为：
y =- (90578142270536885*x^3)/6755399441055744 + (39029822268770605*x^2)/2251799813685248 + (5696372587309794739*x)/11258999068426240 - 20173191042207372827/168884986026393600
在第8个区间(0.4000<=x<=0.6000)内，S(x)的表达式为：
y =(3800*(x - 2/5)^3)/3077 - (88928733325263785*x)/70368744177664 - (13992826802779255*(x - 3/5)^3)/6755399441055744 + 41540526145056859/70368744177664
在第9个区间(0.6000<=x<=0.8000)内，S(x)的表达式为：
y = (68409786027632455*x)/35184372088832 + (54675*(x - 3/5)^3)/160004 - (3800*(x - 4/5)^3)/3077 - 46954228541630179/35184372088832
在第10个区间(0.8000<=x<=1.0000)内，S(x)的表达式为：
y =69962787255772577/35184372088832 - (54675*(x - 1)^3)/160004 - (38868241859560495*x)/17