计算方法上机作业集合

1、第一次 沪a；b; f=Q(x)(x) +10; c=(a+b)/2: while abs(b-a)0.3i0 (-3) if 1(c) *f(b) fprintf C n x = of. f(x )=. of ti , x, f (s); x - 0. 0$0b8t(z)-0. 00067 Al 求得近似根为 0.09058 (2) 上机代码如图所示 x0=0 xO = 0 f=(x) (2-exp Ck) )/10 . riiiLa a.bU0-f (kO)0. 5*10* -3; sG=f (stO): x=sO： end fprintf ( Xq x5f,f (s)5f n 4x,f

2、tx): x= Ck 09064, f(x)=0. 09061 A 得近似根为0.09064; (3) 牛顿法上机代码如下 (z) fixp(x)+10*x-2 : E=esp i x +10;冬 ill数f 的导函数 y-e(i)x-f(x/g(x): hFQ; i=0 ; vhile abs(yCxO-xO)D.5*10(-3 sO=y (xO); 5=50 : end; fprintf C? n s = . 5W : x - 0, 09091 计算所得近似解为0.09091 第三次上机作业 上机作业181页第四题 线性方程组为 rl.1348 3.83261.16513.40171 xl

3、- r 9.5342 053011.7S752.533015435 工2 6.3941 3.4129 4.93178.76431.3142 x3 18.4231 1.23714.999810.67210.0147 x4 16.9237 I1 (1)顺序消元法 A二1.1348,3.8326,1.1651,3.4017;0.5301,1.7875,2.5330,1.5435;3.4129,4. 9317,8.7643,1.3142;1.2371,4.9998,10.6721,0.0147; b=9.5342;6.3941;18.4231;16.9237; 上机代码(函数部分)如下 fun cti

4、 on b = gaus( A,b )%用 b返回方程组的解 B=A,b; n=len gth(b); RA=ra nk(A); RB=ra nk(B); dif=RB-RA; if dif0 disp(此方程组无解); return end if RA=RB if RA=n format long; disp(此方程组有唯一解); for p=1: n-1 for k=p+1: n m=B (k,p)/B(p,p); B(k,p: n+1)=B(k,p: n+1)-m*B(p,p: n+1); end end%顺序消元形成上三角矩阵 b=B(1: n,n+1); A=B(1: n,1: n)

5、; b( n)=b( n)/A( n,n); for q=n-1:-1:1 b(q)=(b(q)-sum(A(q,q+1: n)*b(q+1: n)/A(q,q); end%回代求解 else disp(此方程组有无数组解); end end 上机运行结果为 A二1.1348,3.8326,1.1651,3.4017;0.5301,1.7875,2.5330,1.5435;3.4129,4. 9317,8.7643,1.3142;1.2371,4.9998,10.6721,0.0147; b二9.5342;6.3941;18.4231;16.9237; X二gaus(A,b) 此方程组有唯一解

6、 1.0000 1.0000 1.0000 1.0000 (2)列主元消元法 A=1.1348,3.8326,1.1651,3.4017;0.5301,1.7875,2.5330,1.5435;3.4129,4. 9317,8.7643,1.3142;1.2371,4.9998,10.6721,0.0147; b=9.5342;6.3941;18.4231;16.9237; 函数部分代码如下 function b = lieZhu(A,b )%用 b 返回方程组的解 B=A,b; RA=ra nk(A); RB=ra nk(B); n=len gth(b); dif=RB-RA; format

7、 Io ng; if dif0 disp(该方程组无解); return end if dif=0 if RA=n disp(该方程组有唯一解); c=zeros(1, n); for i=1: n-1 max=abs(A(i,i); m=i; for j=i+1: n if max A二1.1348,3.8326,1.1651,3.4017;0.5301,1.7875,2.53 30,1.5435;3.4129,4.9317,8.7643,1.3142;1.2371,4.99 98,10.6721,0.0147; b=9.5342;6.3941;18.4231;16.9237; X=lieZ

8、hu(A,b) 该方程组有唯一解 X = 1.0000 1.0002 0.999999999999999 0.999999999999999 根据两种方法运算结果，顺序消元法得到结果比列主元消元法得到的 结果精度更高。 (注：matlab使用的是2015b版本，不知道是保留小数位数太少，还 是程序原因，顺序消元输出结果总是等于准确解，请老师指正) 第四次上机作业 7.分析用下列迭代法解线性方程组 4 -1 0 -1 0 -1 4 -1 0 -1 0 -1 4 -1 0 -1 0 -1 4 -1 0 -1 0 -1 4 0 0 -1 0 -1 -1 4 0 -1 0 X3 的收敛性，并求出使o.

9、oool的近似解及相应 的迭代次数。 (1) 雅可比迭代法 解：上机编写的函数如下 function = Jacobi(X,b) %雅可比迭代法 D=diag(diag(X);%得到对角线元素组成的矩阵 B=in v(D)*(D-X); f=in v(D)*b; b(：,：)=0; x仁 B*b+f; num=1; while (n orm(x1-b)0.0001)%判断当前的解是否达到精度要求 b=x1; x1=B*b+f; num=nu m+1; end ; fprintf(求得的解为：n); disp(xl); fprintf(迭代次数：(次n ,num); end 上机运行，结果如下

10、A=4,-1,O,-1,O,0；-1,4,-1,O,-1,O;O,-1,4,-1,O,-1；-1,O,-1,4,-1 ,0；0,-1,0,-1,4,-1；0,0,-1,0,-1,4; b=0;5;-2;5;-2;6; Jacobi(A,b) 求得的解为： 0.999981765074381 1.99995018125674 0.999975090628368 1.99995018125674 0.999975090628368 1.99996353014876 迭代次数：28次 满足要求的解如输出结果所示，总共迭代了28次 (2) 高斯-赛德尔迭代法 上机程序如下所示 fun ctio n =

11、Gauss_Seidel( A,b ) %高斯赛德尔迭代法 D=diag(diag(A); L=D-tril(A); U=D-triu(A); B=i nv(D-L)*U; f=i nv(D-L)*b; b(：,：)=0； xO=B*b+f; num=1; while (norm(x0-b)0.0001) num=nu m+1; b=x0; x0=B*b+f; end ; fprintf(结果为 n); disp(x0); fprintf(迭代次数为：(次n ,num); end A=4,-1,0,-1,0,0;-1,4,-1,0,-1,0;0,-1,4,-1,0,-1;-1,0,-1,4,-

12、1 ,0;0,-1,0,-1,4,-1;0,0,-1,0,-1,4; b=0;5;-2;5;-2;6; Gauss_Seidel(A,b) 结果为 0.999960143810658 1.99995676152139 0.999963508299833 1.99996662162874 0.999972527179715 1.99998400886989 15次 迭代次数为：15次 满足精度要求的解如上述程序打印输出所示，迭代了 (3)SOF迭代法(w依次取 1.334,1.95,0.95) 上机代码如下 fun ctio n = SOR(A,b,w ) %S0!迭代法 D=diag(diag

13、(A); L=D-tril(A); U=D-triu(A); B=i nv(D-w*L)*(1-w)*D+w*U); f=w*i nv(D-w*L)*b; b(:,:)=0; x0=B*b+f; num=1; while (norm(x0-b)0.0001) num=nu m+1; b=x0; x0=B*b+f; end ; fprintf(结果为 n); disp(x0); fprintf(迭代次数为 %dn ,num); end 上机运行 A=4,-1,0,-1,0,0;-1,4,-1,0,-1,0;0,-1,4,-1,0,-1;-1 ,0,-1,4,-1,0;0,-1,0,-1,4,-1

14、;0,0,-1,0,-1,4; b=0;5;-2;5;-2;6; SOR(A,b,1.334) 结果为 1.009 1.99998698322858 1.068 2.053 0.999991858543476 2.69 迭代次数为13 SOR(A,b,1.95) 结果为 0.999984952088107 2.604 0.999959115182729 2.06 1.697 1.99997885113446 迭代次数为241 SOR(A,b,0.95) 结果为 0.999961518309351 1.99995706825231 0.999963054838453 1.999965805720

15、33 0.999971141727589 1.9999827300678 迭代次数为17 由以上输出得到w取值不同的情况下，得到的满足精度要求的结果, 迭代次数分别如输出所示 第五次上机作业 8.从函数表 x 0.0 0.1 0.2 0.3 0.401 0.5 f(x) 0.39894 0.39695 0.39142 0.38138 0.36812 0.35206 出发，用下列方法计算f(0.15),f(0.31) 及f(0.47)的近似值 (1) 分段线性插值 (2) 分段二次插值 (3) 全区间上拉格朗日插值 解：(1)线性插值 编写函数如下 fun ctio n R = div_li n

16、e( x0,y0,x) %线性插值 p=le ngth(xO); n=len gth(y0); m=le ngth(x); if(p=n)%x的个数与y的个数不等 error(数据输入有误，请重新输入); return; else fprintf(线性插值 n); for t=1:m z=x(t); if(zx0(p) fpri ntf(x%d不在所给自变量范围内，无法进行插值,t); con ti nue; else for i=1:p-1 if(z x0=0.0 0.1 0.195 0.3 0.401 0.5; y0=0.39894 0.39695 0.39142 0.38138 0.36

17、812 0.35206; x=0.15 0.31 0.47; div_line(x0,y0,x) 线性插值 ans = 0.394040.380070.35693 即结果为 f(0.15)0.39404,f(0.31)0.38007,f(0.47) 乂 0.35693 (2)分段二次插值 编写的函数如下 fun ctio n R = div2li ne(xO,yO,x) %分段二次插值 p=le ngth(xO); m=le ngth(yO); n=len gth(x); if(p=m) error(输入错误，请重新输入数据 ); end; for t=1: n if(x(t)x0(p) fp

18、rintf(x%d不在所给区间上,t); con ti nue; en d; tag=2; m=abs(x(t)-x0(1)+abs(x(t)-x0(2)+abs(x(t)-x 0(3); for i=3:p-1 sum=abs(x(t)-x0(i-1)+abs(x(t)-x0(i)+abs(x(t)-x0(i+1); if(sum x0=0.0 0.1 0.195 0.3 0.401 0.5; y0=0.39894 0.39695 0.39142 0.38138 0.36812 0.35206; x=0.15 0.31 0.47; div2li ne(x0,y0,x) ans = 0.394

19、48 0.38022 0.35725 即分段二次插值方法下， 0.35725 f(0.15) | = 0.39448,f(0.31)= 0.38022,f(0.47) (3) 上机编写的程序如下 fun cti on R = lagra nge(x0,y0,x) %全区间上拉格朗日插值 p=le ngth(yO); n=le ngth(xO);m=le ngth(x); %计算函数表和x的长度 if p = n error(数据输入有误，请重新输入)； %若函数表的x与y长度不一致则输入有误 else fprintf(拉格朗日插值n); for t=1:m %利用循环计算每个x的插值 s=0.

20、0; z=x(t); for k=1: n p=1； for i=1: n if i=k p=p*(z-x0(i)/(x0(k)-x0(i); end end s=s+y0(k)*p; end %根据拉格朗日插值公式求解y R(t)=s; %输岀插值结果 end end 上机运行结果为 x0=0.0 0.1 0.195 0.3 0.401 0.5; y0=0.39894 0.39695 0.39142 0.38138 0.36812 0.35206; x=0.15 0.31 0.47; lagra nge(x0,y0,x) 拉格朗日插值 ans = 0.394470.380220.35722

21、即分段二次插值方法下， f(0.15)0.39447,f(0.31)0.38022,f(0.47) 怎 0.35722 9.解：上机程序如下，为方便起见，将所有操作分在四个函数中进行 入口函数 function =spline( X,Y,xx,y1_0,y1_18 ) %输出自变量所对应的函数值 M=getM(X,Y,y1_0,y1_18);%得到 M s=xx; k=le ngth(xx); for a=1:k s(xx(a)=getExp(X,Y,M,xx(a)i%算自变量所在小区间对应曲线的表达式，并 根据表达式计算对应的函数值 fprin tf(s(%d)=%fn,xx(a),s(xx

22、(a); %输出打印函数值 end end 获取M fun ctio n M = getM(X,Y,y1_0,y1_1) %得到M n=le ngth(X); a=0*X;b=a;c=a;h=a;f=a; b=b+2; h(2: n) =X(2:n )-X(1: n-1);% h(1)不 用 a(2:n-1)=h(2: n-1)./(h(2: n-1)+h(3: n); c(2:n-1)=1-a(2: n-1); a(n )=1;c(1)=1; Yf(2: n)=Y(2: n)-Y(1: n-1); f(2: n-1)=6*(Yf(3: n)./h(3: n)-Yf(2: n-1)./h(2:

23、 n-1)./(h(2: n-1)+h(3: n); f(1)=6*(Yf(2)/h(2)-y1_0)/h(2); f(n )=6*(y1_1-Yf( n)/h( n)/h( n); M=CalM(n,a,b,c,f);%计算 M end 计算M fun ctio n f = CalM( n,a,b,c,f) %追赶法求解M eps=0.1e-15; %防止参数过小，是的计算误差过大 if abs(b(1)eps disp(除数为0,停止计算); return else c(1)=c(1)/b(1); end %追赶法：根据递推算式计算B for i=2: n-1 b(i)=b(i)-a(i)

24、*c(i-1); if abs(b(i)X(n 5) n1=n5; else n2=n5; end end % %计算y y=M( n1)*(X( n2)-xF3/(6*h( n2)+ M( n2)*(x-X( n1)F3/(6*h( n2); y=y+(Y( n1)-M( n1)*h( n2)*h( n2)/6)*(X( n2)-x)/h( n2); y=y+(Y( n2)-M( n2)*h( n2)*h( n2)/6)*(x-X( n1)/h( n2); % %结束 end 上机试运行，过程如下 X=0.52 3.1 8.0 17.95 28.65 39.62 50.65 78 104.6

25、 156.6 208.6 260.7 312.5 364.4 416.3 468 494 507 520; Y=5.28794 9.413.84 20.2 24.9 28.44 31.135 36.5 36.6 34.6 31.0 26.34 20.9 14.8 7.8 3.7 1.5 0.2; xx=2 4 6 12 16 30 60 110 180 280 400 515; y1_0=1.86548; y1_18=-0.046115;spl in e(X, Y,xx,y1_0,y1_18) s(2)=7.825123 s(4)=10.481311 s(6)=12.363477 s(12)=

26、16.575574 s(16)=19.091594 s(30)=25.386402 s(60)=32.804283 s(110)=36.647878 s(180)=35.917148 s(280)=29.368462 s(400)=16.799784 s(515)=0.542713 根据上述程序运行结果，可得所有自变量对应的函数值，如上输出结 果所示 第六次上机作业 10.已知一组实验数据 i i 2 3 4 C 6 7 8 9 J X i 3 4 5 6 7 8 C 1 0 y -i 5 4 2 1 1 2 3 4 0 试用最小二乘法求它的多项式拟合曲线，并求出最低点位置 解：试用matla

27、b绘图命令，将以上个点绘制在坐标图上，如下图所 10f 10 该函数图像形如二次函数的图像，现使用二次拟合 程序如下 function xmin,ymin = SecondFitting(x,y) %最小二乘法，二次拟合，形如a+bx+cxA2 x=x; y=y; m=le ngth(x); A(:,1)=on es(m,1); A(:,2)=x; A(:,3)=x.a2; cc=Ay; a=cc(1); b=cc(2); c=cc(3); fprintf(拟合的曲线为 %f%fx+%fxA2n,a,b,c); xx=1:0.01:10; yy二cc(1)+cc(2)*xx+cc(3)*xx.

28、A2; plot(x,y,r*,xx,yy); ymi n=min(yy); c1=a-y min; p=c b c1; xmi n=roots(p); fprintf(最低点的横纵坐标分别为); end 上机运行 x=1 3 4 5 6 7 8 9 10; y=10 5 4 2 1 1 2 3 4; Secon dFitt in g(x,y) 拟合的曲线为 13.459664-3.605309x+0.267571xA2 最低点的横纵坐标分别为ans = 6.74 6.7342 图像如下图所示 11 1234567B910 根据程序运行结果，得到拟合方程为 y=13.459664-3.6053

29、09x+0.267571x9,最低点坐标为(6.74 , 6.7342)。 计算方法第七次上机作业 2 2 龙 y 15.用龙贝格算法计算椭圆400+100=1的周长，使误差不超过 解：椭圆周长L计算公式为 L=4 上机程序如下所示 (注：把每次计算结果存储于一个二维数组之中，为了输出如书中所 示的表格形式，每一次结果的下标有一定的调整) %fun表示被积函数，a为区间下限，b为上限, fun cti on = Romberg( fun, a,b,e ) e为精度要求 %龙贝格算法求椭圆周长 T=zeros(); T(1,1)=(b-a)/2*(feval(fu n,a)+feval(fu n

30、,b); T(2,1)=0.5*T(1,1)+(b-a)/2*(feval(fu n, a+(b-a)/2); T(2,2)=(4*T(2,1)-T(1,1)/3; T(3,1)=0.5*T(2,1)+(b-a)/4*(sum(feval(fu n,a+(2*(1:2)*(b-a)/4); T(3,2)=(4*T(3,1)-T(2,1)/3; T(3,3)=(16*T(3,2)-T(2,2)/15; T(4,1)=0.5*T(3,1)+(b-a)/8*(sum(feval(fu n,a+(2*(1:4)*(b-a)/8); T(4,2)=(4*T(4,1)-T(3,1)/3; T(4,3)=(

31、16*T(4,2)-T(3,2)/15; T(4,4)=(64*T(4,3)-T(3,3)/63; T(5,1)=0.5*T(4,1)+(b-a)/16*(sum(feval(fu n,a+(2*(1:8)*(b-a)/16); T(5,2)=(4*T(5,1)-T(4,1)/3; T(5,3)=(16*T(5,2)-T(4,2)/15; T(5,4)=(64*T(5,3)-T(4,3)/63; %以上为求岀序列表中前五行的值 k=5; while(4*abs (T(k,4)-T(k-1,4)e) k=k+1; T(k,1)=0.5*T(k-1,1)+(b-a)/(2A(k-1)*(sum(f

32、eval(fu n,a+(2*(1:(2A(k-2) )*(b-a)/0(k-1); T(k,2)=(4*T(k,1)-T(k-1,1)/3; T(k,3)=(16*T(k,2)-T(k-1,2)/15; T(k,4)=(64*T(k,3)-T(k-1,3)/63; end disp(4*T);%区间为四分之一圆周，所以乘以4输岀为椭圆周长 上机运行，结果如下(复制数字结果，排版容易乱，截图如下) 运行命令 fun 二(x)power(300*(si n(x).A2)+100),0.5); Romberg(fu n, 0,pi/2,0.0001) f L3i=E,Cx pciFer i , 3M，x- ;s in (x)2J+10C : , Cl 61 ; RcnbcrE( fun，0, p i/2, (L (0 01) 9匚2右厂9別：旳詐 0 0 0 S-.制飾149“阳眸 0 0 W W0IL27E03 1D7.3IB906M4793 107.