计算方法实验指导书学生用

1、实验一 插值方法一、实验目的和要求（1） 理解插值的基本概念，掌握各种插值方法，包括（2） 通过实验进一步理解并掌握各种插值方法。二、实验器材PC 机、Turbo c 或 C+、WINDOWS 系统三、实验内容插值、插值、分段插值。任务（一）1、 实验内容插值法：用插值多项式求函数近似值。设函数 y=f(x)在给定的两两互异节点 x0，x1，xn 上的函数值为 y0，y1，yn，求作一个次数n 的多n项式 Pn(x)=a0+a1x+anx 使它满足：Pn(xi)=nyi(i=0,1,2,n)。p (x) = y l (x) 其中l (j为插值公式。nk kkx - xk =0j =0 j kk

2、j2、 实验步骤#include #define N 3void main()float XN+1,YN+1,x;double L,LI; int i,j;printf(Please input Xi,Yi:n); for(i=0;i=N;i+)scanf(%f,%f,&Xi,&Yi); printf(Please input x:n); scanf(%f,&x);L=0;for(i=0;i=N;i+)LI=1;for(j=0;j=N;j+)if(j=i) continue;LI=(LI*(x-Xj)/(Xi-Xj); L=L+LI*Yi;printf(nN=%dnL(%f)=%fn,N,x,

3、L);3、 实验结果与分析Please input Xi,Yi: 0.56160,0.827410.56280,0.826590.56401,0.825770.56521,0.82495Please input x: 0.5635N=3L(0.563500)=0.826116任务（二）1、实验内容插值方法：造出差，并用插值多项式求函数近似值。f (x) 在 x0，x1，xk 上的k 阶差商 f - f f (x j )kf k = j =0k -1kx - xk(x j - xi )i =0 i j0kk - f x0 ,f x0 ,k -1 =- xk -1xkNn (x) = f (x0

4、) + f x0 ,0 ) +L+ f n (-n-1 )=f(x0)+(x-x0)fx0,x1+(x-x1)fx0,x1,x2+(x-x2).fx0,.,xn-1+(x-xn-1)fx0,.,xn.构造差x0 x1 x2 x3x4如下: y0y1 y2 y3y4fx0,x1fx0,x2fx0,x3fx0,x4fx0,x1,x2fx0,x1,x3fx0,x1,x4fx0,x1,x2,x3fx0,x1,x2,x4fx0,x1,x2,x3,x42实验步骤#include #define N 4void main()float XN+1,YN+1,FN+1;float x1,x2,N1,N2; in

5、t i,k;printf(please input Xi,Yin); for(i=0;i=N;i+)scanf(%f,%f,&Xi,&Yi);printf(0t%.2ft%.5fn,X0,Y0); for(k=1;k=N;k+)F0=Yk;printf(%dt%.2ft%.5f ,k,Xk,Yk); for(i=0;i=0;i-)N1=(x1-Xi)*N1+Yi;N2=(x2-Xi)*N2+Yi;printf(N(%.3f)=%fn,x1,N1);printf(N(%.3f)=%fn,x2,N2);3实验结果与分析please input 0.4,0.410750.55,0.578150.65

6、,0.696750.8,0.888110.9,1.02652Xi,Yi:01234please0.400.550.650.800.90input0.410750.578150.696750.888111.026521.1160001.1440001.1934001.2315400.2800020.3096000.3301150.1973190.2004510.031312x1,x2:0.596,0.895N(0.596)=0.631918 N(0.895)=1.019368本程序给定插值点序列，构造插值多项式。输入要计算的函数点 x 并计算函数值，利用插值公式，当增加一个节点时，只需在后面多计

7、算一项，而前面的计算仍有用；另一方面插值公式的各项系数恰好又是各阶差商，而各阶差商可用差商公式来计算。四、出现的问题及解决方案Newton 插值插商的计算,把原插商公式变形可简化运算: f - f f (x j )kf k = j =0k -1kx - xk(x j - xi )i =0 i j0kk - f x0 ,f x0 ,k -1 =- xk -1xk五、结论通过实验进一步理解并掌握了各种插值方法。实验二 数值一、实验目的及要求（1）进一步理解数值的概念。（2）掌握数值的基本方法、梯形求积公式、复化求积公式。（3）掌握自动变步长的数值二、实验器材算法。PC 机、Turbo c 或 C+

8、、WINDOWS 系统三、实验任务任务（一）用复化梯形法计算。1、实验内容把区间a,bn 等分，分点 xk = a + kh （k=0,1,n），h = (b - a)求积公式计算被积函数 f (x) 在其上的，即：n ，对每个小区间xk , xk +1 用梯形f (x)dx 1 (x2xk +1- x ) f (x ) + f (x)k +1kk +1kxk则：n-1I =f (x)dx = n-1f (x)dx 1 (xbxk +1- x ) f (x ) +f (x)k +1kkk +1k =0 2axkk =0n-1= h f (a) + 2 f (x ) + f (b) = Tkn2

9、k =1Tn 称为复化求积公式，下标 n 表示将区间 n 等分。2、实验步骤11dx/计算1 + x20#include #include #define F(x) 1/(1+x*x) void main()float A,B,H,X,T; int N,i;A=0;B=1;printf(please input N:n); scanf(%d,&N);H=(B-A)/N; X=A;T=0;for(i=1;i=N-1;i+)X=X+H;T=T+F(X); T=(F(A)+2*T+F(B)*H/2;printf(T(%d)=%fn,N,T);3 实验结果与分析please input N:128T(

10、128)=0.785396please input N:256T(256)=0.785398please input N:512T(512)=0.785398本程序运行三次,仅当 n=29=512,n=28=256 时有六位数相同。任务（二）用自动变步长梯形法计算。1、实验内容将区间a，bn 等分，分点 xk = a + kh ， h = (b - a)值。n ，（k=0,1,n）， Tn 表示用复化梯形法求得的= 1 (2+ 1 h ，该小间二分前后的两个积2小区间x , x，记该区间的中点为 xkk +1k + 1k2分值（用梯形法求）分别记为Tk1 和Tk 2 ，则：T = h f (x

11、 ) + f (x)k +1k1k2T = h f (x ) + f (x) + h f (x) + f (x ) = h f (x ) + 2 f (x) + f (x)k +1k +1k 2kk + 1k + 1kk + 1444222= 1 T+ h2n-1n-11+ h2 k 2k1故： Tf (x)= T，Tf (x) 。k 2k1k + 12k + 1222k =0k =0h n-11k =0即： T2n = 2 Tn +f (x)k +122这一公式是递推形式的，式中的 h 是二分前的步长即 h = (b - a) n 。2、实验步骤#include #include #defi

12、ne F(x) 1/(1+x*x) void main()float A,B,E,H,S,X,T1,T2;int N,k; A=0; B=1;printf(please input E:n); scanf(%e,&E);H=B-A;N=1;T1=(F(A)+F(B)*H/2;printf(n%ft,T1); while(1)X=A+H/2; S=0;for(k=0;k=N-1;k+)S=S+F(X); X=X+H; T2=(T1+S*H)/2;printf(%ft,T2); N=N*2;if(fabs(T2-T1)E) break; H=H/2;T1=T2;printf(nT(%d)=%fn,

13、N,T2);3、实验结果与分析please input 1e-60.7500000.785357E:0.7750000.7853880.7827940.7853960.7847470.7853970.7852350.785398T(512)=0.785398本程序实质上等分了 9 次，区间被分为 29=512 等分，与复化梯形法实质一样。四、出现的问题及解决方案，使用(T2-T1)E 来精度，实质上等分了 9 次，区间被分为 29=512 等分。在自动变步长梯形法五、结论掌握了数值的基本方法、梯形求积公式、复化求积公式，自动变步长的数值算法。实验三 常微分方程的数值解法一、实验目的及要求（1）

14、理解常微分方程初值问题数值解的概念。（2）掌握几种简单的数值解法：方法，改进的方法。、线性多步法。（3）了解常微分方程初值问题数值解法中二、实验器材PC 机、Turbo c 或 C+、WINDOWS 系统三、实验内容任务（一）用改进法解常微分方程的初值问题。1、实验内容 y p= yn + hf (xn , yn )-y= y + hf (x , y )yn+1nnn= y + hf (x, y )或n+1ph-cn yn+1= yn + 2 f (xn , yn ) + f (xn+1 , yn+1 )1 y=( y + y )n+1pc2这即为改进的1、 实验步骤法公式。y = -xy 2

15、y(0) = 2(0 x 5)/用改进法解#include void main()float X,Y,Xn,Yp,Yc,Yt,H; int N,i;printf(please input X,Y,N,Xn:n); scanf(%f,%f,%d,%f,&X,&Y,&N,&Xn); H=(Xn-X)/N;printf(N=%d,H=%fn,N,H); printf(ittxttyttytn); printf(%dtt%ftt%ftt%fn,0,X,Y,Y); for(i=1;i=N;i+)Yp=Y-H*X*Y*Y; X=X+H;Yc=Y-H*X*Yp*Yp; Y=(Yp+Yc)/2; Yt=2/(

16、1+X*X);printf(%dtt%ftt%ftt%fn,i,X,Y,Yt);3实验结果与分析please input X,Y,N,Xn: 0,2,20,5N=20,H=0.250000i 012345x 0.0000000.2500000.5000000.7500001.0000001.250000y 2.0000001.8750001.5938911.2823901.0096210.793188yt 2.0000001.8823531.6000001.2800001.0000000.780488678910111213141516171819201.5000001.7500002.000

17、0002.2500002.5000002.7500003.0000003.2500003.5000003.7500004.0000004.2500004.5000004.7500005.0000000.6281510.5037300.4096670.3378650.2823570.2388570.2043000.1764900.1538360.1351750.1196420.1065920.0955300.0860800.0779480.6153850.4923080.4000000.3298970.2758620.2335770.2000000.1729730.1509430.1327800

18、.1176470.1049180.0941180.0848810.076923输入 x0=0，y0=2，n=20，xn=5输出近似解 y(i)，准确解 yt(i),本题准确y=2/(1+x2)。任务（二）用1、 实验内容校正公式，求解初值问题。- 校正公式：h24= y +(55y - 59y+ 3= f (x, y0y0y0) 预报：n +1n+n+11n -1n+1nnh24校正： y= y + 19y - 5 y=0y(9 yf (x, y) n +1n +1nnn -1n+1n+1n+1公式： y+ k)4kkk k2、实验步骤/用- 校正公式解y#include void main(

19、)float a,b,h,k1,k2,k3,k4;float x101,y0101,f0101,y101,f101,yc101;int n,i;printf(please input a,b,n:n); scanf(%f,%f,%d,&a,&b,&n);h=(b-a)/n;printf(a=%.2f,b=%.2f,n=%d,h=%.2fn,a,b,n,h); printf(itxittyitttycin);x0=a;y0=2;f0=-x0*y0*y0;yc0=2/(x0*x0+1);printf(%dt%.2ftt%ftt%f,0,x0,y0,yc0); for(i=1;i=3;i+)xi=x

20、i-1+h; k1=fi-1;k2=-(xi-1+h/2)*(yi-1+h*k1/2)*(yi-1+h*k1/2);k3=-(xi-1+h/2)*(yi-1+h*k2/2)*(yi-1+h*k2/2);k4=-xi*(yi-1+h*k3)*(yi-1+h*k3); yi=yi-1+(k1+2*(k2+k3)+k4)*h/6; fi=-xi*yi*yi;yci=2/(xi*xi+1); printf(n%dt%.2ftt%ftt%f,i,xi,yi,yci);for(i=4;i=n;i+)xi=xi-1+h;y0i=yi-1+(55*fi-1-59*fi-2+37*fi-3-9*fi-4)*h/

21、24; f0i=-xi*y0i*y0i;yi=yi-1+(9*f0i+19*fi-1-5*fi-2+fi-3)*h/24; fi=-xi*yi*yi;yci=2/(xi*xi+1); printf(n%dt%.2ftt%ftt%f,i,xi,yi,yci);printf(n);3、实验结果与分析please input a,b,n: 0,5,20a=0.00,b=5.00,n=20,h=0.25i 012345678910xi 0.000.250.500.751.001.251.501.752.002.252.50yi 2.0000001.8823081.5998961.2799480.996

22、4360.7806570.6180330.4945020.4015730.3310350.276622yci 2.0000001.8823531.6000001.2800001.0000000.7804880.6153850.4923080.4000000.3298970.275862111213141516171819202.753.003.253.503.754.004.254.504.755.000.2340820.2003440.1732120.1511130.1329030.1177390.1049870.0941710.0849230.076956yci。本题准确0.2335770

23、.2000000.1729730.1509430.1327800.1176470.1049180.0941180.0848810.076923y=2/(1+x2)。yi 为近似解，准确四、出现的问题及解决方案在-校公式计算值 y1，y2，y3。，计算 yi+1 时要用到 fi，fi-1，fi-2，fi-3，因此，开始必须借助于五、结论掌握了常微分方程数值解法：方法，改进的方法，-校正公式。实验四 线性方程组的数值解法一、实验目的及要求（1）理解线性方程组的消去法。（2） 掌握常用的几种迭代公式。（3） 掌握矩阵分解法。（4） 掌握迭代收敛性和误差估计。二、实验器材PC 机、Turbo c 或

24、C+、WINDOWS 系统三、实验任务任务（一）用消去法求解线性方程组。1、实验内容n对于线性代数方程组 AX=b 即： aij x j = bi(i = 1,2,L, n)j =1Gauss 消元法的计算步骤分为消元和两个过程。 a11a12 a22 0M0a13 a23LLLOLa1nb1(2) a11a12LLOLa1nb1 a(2)(2)(2)00M0a2nb2aab消元：21222n2a(3)a(3)(3)bM MMM33M03nM3aMaab n1n n2nna(n)(n)bn nn第一步： 第1行 - ai1 + 第i行, i = 2,L, na11a1n a11b1 a11a1

25、2LLOLa12LLOLa1nb1(2) aa(2)a(2)aab0b21222n2 22M2nM2MMMMM M aa(2)a(2)(2)aab0b n1n n n2nnn 2nn- a(2)第二步： 第2行i2 + 第i行, i = 3,L, na(2)22 a11a12 a22 0M0LLLa13 a23LLLOLb1a1nb1(2) a11a12LLOLa1nb1(2)(2)(2)(2)00M0a13a2nb2a(2)a(2)0M0b22M2nM2Ma(3)a(3)(3)b33M3nM3Ma(2)a(2)(2)bn 2nnna(3)a(3)(3)bn n3a1nnn a11a12a(2

26、)a(2)a(2)b(2) 00M0220M0232n2(3)，化为a(3)a(3)b经过 n-133M03n3MOMa(n)(n)Lbn nn顺序消元的计算公式为：a (k ) = a (k -1) - a (k -1) a (k -1) / a (k -1)ijijikkjkk(i = k + 1, k + 2,L, n; j = k, k + 1,L, n)(k ) = b(k -1) - a (k -1) b(k -1) / a(k -1)biiikkkk(k = 1,2,L, n -1)x= b(n-1) / a(n-1)nnnn的计算公式为： xn- 其中i = n - 1, n

27、- 2,L,1= (b(i-1)a(i-1) x) / a(i-1)iiijjiij =i+12、实验步骤#include #include #include #define N 3void main()float AN+1N+2,S,L;int i,j,k;printf(please input Aij:n); for(i=1;i=N;i+)for(j=1;j=N+1;j+) scanf(%f,&Aij);printf( (A,b) 1 n); printf(n);for(i=1;i=N;i+)for(j=1;j=N+1;j+) printf(%.2f ,Aij);printf(n);for

28、(k=1;k=N-1;k+)if(fabs(Akk)1e-6)printf(fail!); exit(0);for(i=k+1;i=N;i+)L=Aik/Akk;printf( L%d%d=%f,i,k,L); for(j=k;j=N+1;j+)Aij=Aij-L*Akj;printf(n (A,b) %dn,k+1);printf(n);for(i=1;i=N;i+)for(j=1;j=1;i-)S=0;for(j=i+1;j=N;j+)S=S+Aij*AjN+1;AiN+1=(AiN+1-S)/Aii;printf(nThe Solution of Equations is for(i=1

29、;i=N;i+)printf(%f,AiN+1);X=);3、实验结果与分析please input Aij:2.55.38.12.39.61.7-5.1 3.71.5 3.8-4.35.5(A,b) 1. 2.500 2.3005.3009.600-5.1003.7001.5003.8008.1001.700L21=2.120000 (A,b) 2.-4.3005.500L31=3.2400002.500-0.000-0.0002.3004.724-5.752-5.10012.31212.2243.700-4.044-6.488L32=-1.217612 (A,b) 3.2.500-0.00

30、0-0.0002.3004.724-0.000-5.10012.31227.2153.700-4.044-11.412The solution of Equations is X= 0.4067050.236817 -0.419325本程序在执行时首先读入方程组的增广矩阵 An(n+1),其元素按行输入，通过消元后，从最后一列 Ain+1中输出方程组的。任务（二）用矩阵分解法求解线性方程组。1、实验内容给定 n 阶非奇异矩阵 A，两个 n 阶矩阵1lr11r12 rLLOr1n1MOrL = 21R = 2n 22 MML 1ln1ln2rnn 使 A=LR。则有：LRX=b,令 RX=Y，有

31、：LY=b ，RX=Y。由式 LY=b 求出Y，由式 RX=Y 再求出解 X。分解计算公式如下：( j = 1,2,L, n)(i = 2,3,L, n)= a1 jr1 j(1)(2)li1 = ai1 / a11k -1rkj = akj - lkqrqj( j = k, k +1,L, n)(3)q=1k -1- liq rqk / rkk= aiklik(4)q=1解方程组 LY=b,(5)计算公式为: y1 = b1 yn- ly= b(k=2,3,n)(7)kkkq qq=k +1解方程组 RX=Y,(6)计算公式为:= ynxnxn- l x ) / r(k=n-1,n-2,1)

32、(8)= ( ykkkq qkkq=k +1a11a1222a1323a33MLLLOLaa21设计为紧凑格式： a3132MMaann 2n32、实验步骤#include #define N 4void main()float AN+1N+2;int i,j,k,q; printf(please input for(i=1;i=N;i+)Aij:n);for(j=1;j=N+1;j+)scanf(%f,&Aij);printf(A,b)n);printf(.n);for(i=1;i=N;i+)for(j=1;j=N+1;j+)printf( %.1f ,Aij); printf(n);for

33、(i=2;i=N;i+)/*按公式(2)计算 li1*/Ai1=Ai1/A11;for(k=2;k=N-1;k+)for(j=k;j=N+1;j+) for(q=1;q=k-1;q+)Akj=Akj-Akq*Aqj;/*按公式(3)计算 rkj*/for(i=k+1;i=N;i+)for(q=1;q=k-1;q+)/*按公式(4)计算 lik*/Aik=Aik-Aiq*Aqk;Aik=Aik/Akk;for(j=N;j=N+1;j+) for(q=1;q=1;i-)for(q=i+1;q=N;q+)AiN+1=AiN+1-Aiq*AqN+1;AiN+1=AiN+1/Aii;printf(nX=

34、); for(i=1;i=N;i+)printf( %f ,AiN+1);printf( n );/*按公式(8)计算 X*/3、实验结果与分析please input Aij 1 2 -12 8 275 4 7 -2-3 7 9 56 -12 -8 (A,b)4113 49.1.05.0-3.06.02.04.07.0-12.0-12.07.09.0-8.08.0-2.05.03.027.04.011.049.0X=3.000000-2.0000011.000000 5.000000四、出现的问题及解决方案，为检查 akk(k)=0 是否成立，设置了一个检查口，一旦遇到 akk(k)=0，便

35、使其打印出 fail，以在消去法示简单消去法的失败。五、结论理解线性方程组的消去法和矩阵分解法。实验五 方程求根的数值解法一、实验目的及要求（1）理解非线性方程求根的数值解法，掌握二分法。（2）掌握二、实验器材。PC 机、Turbo c 或 C+、WINDOWS 系统三、实验任务任务（一）用二分法求解非线性方程的根。1、 实验内容用二分法求实根的近似值的基本思路，就是逐步将含有根的区间二分，通过函数值的符号，逐步对半缩小有根区间，直到区间缩小到容许误差范围之内，然后取区间的中点为根的近似值。2、实验步骤/求方程 2x3-4x2+3x-6=0 在(-10,+10)的根#include #include#define F(x) x*(2*x-4)*x+3)-6void main()float x0,x1,x2,fx0,fx1,fx2; doprintf(please input x1,x2:n); scanf(%f,%f,&x1,&x2); fx1=F(x1);fx2=F(x2);while(fx1*fx20); dox0