计算方法上机作业集合(共28页).docx

精选优质文档-倾情为你奉上第一次&第二次上机作业上机作业：1.在Matlab上执行：>> 5.1-5-0.1和>> 1.5-1-0.5给出执行结果，并简要分析一下产生现象的原因。解：执行结果如下： 在Matlab中，小数值很难用二进制进行描述。由于计算精度的影响，相近两数相减会出现误差。2.（课本181页第一题）解：（1）n=0时，积分得I0=ln6-ln5，编写如下图代码从以上代码显示的结果可以看出，I20的近似值为0.7465（2）In=01xn5+xdx,可得01xn6dx01xn5+xdx01xn5dx,得 16(n+1)In15(n+1),则有1126I201105,取I20=1105,以此逆序估算I0。代码段及结果如下图所示 结果是从I19逆序输出至I0，所以得到I0的近似值为0.0227。（3）从I20估计的过程更为可靠。首先根据积分得表达式是可知，被积函数随着n的增大，其所围面积应当是逐步减小的，即积分值应是随着n的递增二单调减小的，（1）中输出的值不满足这一条件，（2）满足。设Sn表示In的近似值，Sn-In=(-5)n(S0-I0)(根据递推公式可以导出此式),可以看出，随着n的增大，误差也在增大，所以顺序估计时，算法不稳定性逐渐增大，逆序估计情况则刚好相反，误差不断减小，算法逐渐趋于稳定。2.（课本181页第二题）（1）上机代码如图所示求得近似根为0.09058（2）上机代码如图所示得近似根为0.09064；（3）牛顿法上机代码如下计算所得近似解为0.09091第三次上机作业上机作业181页第四题线性方程组为1.13483.83260.53011.78751.16513.40172.53301.54353.41294.93171.23714.99988.76431.67210.0147x1x2x3x4=9.53426.9237（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;上机代码（函数部分）如下function b = gaus( A,b )%用b返回方程组的解B=A,b;n=length(b);RA=rank(A);RB=rank(B);dif=RB-RA;if dif>0 disp('此方程组无解'); returnendif 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); 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('此方程组有无数组解'); endend上机运行结果为>> 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)此方程组有唯一解X = 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=rank(A);RB=rank(B);n=length(b);dif=RB-RA;format long;if dif>0 disp('该方程组无解'); returnendif 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<abs(A(j,i) max=abs(A(j,i); m=j; end end %求出每一次消元时绝对值最大的一行的行号 if m=i for k=i:n c(k)=A(i,k); A(i,k)=A(m,k); A(m,k)=c(k); end d1=b(i); b(i)=b(m); b(m)=d1;%函数值跟随方程一起换位置 end for k=i+1:n for j=i+1:n A(k,j)=A(k,j)-A(i,j)*A(k,i)/A(i,i); end b(k)=b(k)-b(i)*A(k,i)/A(i,i); A(k,i)=0; end end %完成消元操作，形成上三角矩阵 b(n)=b(n)/A(n,n); for i=n-1:-1:1 sum=0; for j=i+1:n sum=sum+A(i,j)*b(j); end b(i)=(b(i)-sum)/A(i,i) ;%回代求解其他未知数 end endelse disp('此方程组有无数组解');endend上机运行，结果为>> 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=lieZhu(A,b)该方程组有唯一解X = 1.0000 1.0002 0.9999 0.9999>>根据两种方法运算结果，顺序消元法得到结果比列主元消元法得到的结果精度更高。（注：matlab使用的是2015b版本，不知道是保留小数位数太少，还是程序原因，顺序消元输出结果总是等于准确解，请老师指正）第四次上机作业7.分析用下列迭代法解线性方程组 4-1-1 4 0-1-1 0 0 0-1 0 0-1-1 0 4-1-1 4 0 -1-1 0 0-10 0 0-1-1 0 4-1-1 4 x1x2x3x4x5x6 = 0 5-2 5-2 6 的收敛性，并求出使X(k+1)-X(k)20.0001的近似解及相应的迭代次数。（1） 雅可比迭代法解：上机编写的函数如下function = Jacobi(X,b)%雅可比迭代法D=diag(diag(X);%得到对角线元素组成的矩阵B=inv(D)*(D-X);f=inv(D)*b;b(:,:)=0;x1=B*b+f;num=1;while(norm(x1-b)>0.0001)%判断当前的解是否达到精度要求 b=x1; x1=B*b+f; num=num+1;end;fprintf('求得的解为：n');disp(x1);fprintf('迭代次数：%d次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,-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.4381 1.674 0.8368 1.674 0.8368 1.876迭代次数：28次满足要求的解如输出结果所示，总共迭代了28次（2） 高斯-赛德尔迭代法上机程序如下所示function =Gauss_Seidel( A,b )%高斯赛德尔迭代法D=diag(diag(A);L=D-tril(A);U=D-triu(A);B=inv(D-L)*U;f=inv(D-L)*b;b(:,:)=0;x0=B*b+f;num=1;while(norm(x0-b)>0.0001) num=num+1; b=x0; x0=B*b+f;end;fprintf('结果为n');disp(x0);fprintf('迭代次数为：%d次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,-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.0658 1.139 0.9833 1.874 0.9715 1.989迭代次数为：15次满足精度要求的解如上述程序打印输出所示，迭代了15次（3） SOR迭代法（w依次取1.334,1.95,0.95）上机代码如下function = SOR(A,b,w )%SOR迭代法¨D=diag(diag(A);L=D-tril(A);U=D-triu(A);B=inv(D-w*L)*(1-w)*D+w*U);f=w*inv(D-w*L)*b;b(:,:)=0;x0=B*b+f;num=1;while(norm(x0-b)>0.0001) num=num+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;0,0,-1,0,-1,4;>> b=0;5;-2;5;-2;6;>> SOR(A,b,1.334)结果为 1.009 1.858 1.068 2.053 0.3476 2.69迭代次数为13>> SOR(A,b,1.95)结果为 0.8107 2.604 0.2729 2.06 1.697 1.446迭代次数为241>> SOR(A,b,0.95)结果为 0.9351 1.231 0.8453 1.033 0.7589 1.78迭代次数为17由以上输出得到w取值不同的情况下，得到的满足精度要求的结果，迭代次数分别如输出所示第五次上机作业8.从函数表x0.00.10.20.30.4010.5f(x)0.398940.396950.391420.381380.368120.35206出发，用下列方法计算f(0.15),f(0.31)及f(0.47)的近似值（1） 分段线性插值（2） 分段二次插值（3） 全区间上拉格朗日插值解：（1）线性插值编写函数如下function R = div_line( x0,y0,x )%线性插值p=length(x0);n=length(y0);m=length(x);if(p=n)%x的个数与y的个数不等 error('数据输入有误，请重新输入'); return;else fprintf('线性插值n'); for t=1:m z=x(t); if(z<x0(1)|z>x0(p) fprintf('x%d不在所给自变量范围内，无法进行插值',t); continue; else for i=1:p-1 if(z<x0(i+1) break; end; end; R(t)=y0(i)*(x(t)-x0(i+1)/(x0(i)-x0(i+1)+y0(i+1)*(x(t)-x0(i)/(x0(i+1)-x0(i); end; end; 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;>> div_line(x0,y0,x)线性插值ans = 0.39404 0.38007 0.35693即结果为f(0.15)0.39404,f(0.31) 0.38007,f(0.47) 0.35693（2）分段二次插值编写的函数如下function R = div2line(x0,y0,x)%分段二次插值p=length(x0);m=length(y0);n=length(x);if(p=m) error('输入错误，请重新输入数据');end;for t=1:n if(x(t)<x0(1)|x(t)>x0(p) fprintf('x%d不在所给区间上',t); continue; end; tag=2; m=abs(x(t)-x0(1)+abs(x(t)-x0(2)+abs(x(t)-x0(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<m) m=sum; tag=i; end end; fprintf('tag=%dn',tag);R(t)=y0(tag-1)*(x(t)-x0(tag)*(x(t)-x0(tag+1)/(x0(tag-1)-x0(tag)*(x0(tag-1)-x0(tag+1)+y0(tag)*(x(t)-x0(tag-1)*(x(t)-x0(tag+1)/(x0(tag)-x0(tag-1)*(x0(tag)-x0(tag+1)+y0(tag+1)*(x(t)-x0(tag-1)*(x(t)-x0(tag)/(x0(tag+1)-x0(tag-1)*(x0(tag+1)-x0(tag);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;>>div2line(x0,y0,x)ans = 0.39448 0.38022 0.35725即分段二次插值方法下，f(0.15)0.39448,f(0.31) 0.38022,f(0.47) 0.35725（3）上机编写的程序如下function R = lagrange(x0,y0,x)%全区间上拉格朗日插值p=length(y0);n=length(x0);m=length(x); %计算函数表和x的长度if p = n error('数据输入有误，请重新输入'); %若函数表的x与y长度不一致则输入有误 else fprintf('拉格朗日插值n'); for t=1:m %利用循环计算每个x的插值 s=0.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 %根据拉格朗日插值公式求解yR(t)=s; %输出插值结果 endend上机运行结果为>> 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;>> lagrange(x0,y0,x)拉格朗日插值ans = 0.39447 0.38022 0.35722即分段二次插值方法下，f(0.15)0.39447,f(0.31) 0.38022,f(0.47) 0.357229.解：上机程序如下，为方便起见，将所有操作分在四个函数中进行入口函数function =spline( X,Y,xx,y1_0,y1_18 )%输出自变量所对应的函数值M=getM(X,Y,y1_0,y1_18);%先得到Ms=xx;k=length(xx);for a=1:k s(xx(a)=getExp(X,Y,M,xx(a);%计算自变量所在小区间对应曲线的表达式，并根据表达式计算对应的函数值 fprintf('s(%d)=%fn',xx(a),s(xx(a); %输出打印函数值endend获取Mfunction M = getM(X,Y,y1_0,y1_1)%得到Mn=length(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: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);%计算Mend计算Mfunction f = CalM(n,a,b,c,f)% 追赶法求解Meps=0.1e-15; %防止参数过小，是的计算误差过大if abs(b(1)<eps disp('除数为0，停止计算'); returnelse c(1)=c(1)/b(1);end%追赶法：根据递推算式计算for i=2:n-1 b(i)=b(i)-a(i)*c(i-1); if abs(b(i)<eps disp('除数为0，停止计算'); return else c(i)=c(i)/b(i); endendb(n)=b(n)-a(n)*c(n-1);%追赶法:根据递推算式计算f(1)=f(1)/b(1);for i=2:n f(i)=(f(i)-a(i)*f(i-1)/b(i);end%以下求解Ux=y, x的值存入ffor i=n-1:-1:1 f(i)=f(i)-c(i)*f(i+1);endreturnend得到自变量所在区间的表达式，并求自变量对应的函数值function y = getExp(X,Y,M,x)%根据X、Y、M计算表达式，并根据表达式计算对应的函数值n=length(X);h(2:n)=X(2:n)-X(1:n-1); %判断x落在哪个小区间n1=1;n2=n;while n2=n1+1 n5=fix(n1+n2)/2); if x>X(n5) n1=n5; else n2=n5; endend%计算yy=M(n1)*(X(n2)-x)3/(6*h(n2)+ M(n2)*(x-X(n1)3/(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 156.6 208.6 260.7 312.5 364.4 416.3 468 494 507 520;>> Y=5.28794 9.4 13.84 20.2 24.9 28.44 31.1 35 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.;spline(X,Y,xx,y1_0,y1_18)s(2)=7.s(4)=10.s(6)=12.s(12)=16.s(16)=19.s(30)=25.s(60)=32.s(110)=36.s(180)=35.s(280)=29.s(400)=16.s(515)=0.根据上述程序运行结果，可得所有自变量对应的函数值，如上输出结果所示第六次上机作业10.已知一组实验数据i123456789xi1345678910yi1054211234试用最小二乘法求它的多项式拟合曲线，并求出最低点位置。解：试用matlab绘图命令，将以上个点绘制在坐标图上，如下图所示该函数图像形如二次函数的图像，现使用二次拟合程序如下function xmin,ymin = SecondFitting(x,y)%最小二乘法，二次拟合,形如a+bx+cx2x=x'y=y'm=length(x);A(:,1)=ones(m,1);A(:,2)=x;A(:,3)=x.2;cc=Ay;a=cc(1);b=cc(2);c=cc(3);fprintf('拟合的曲线为%f%fx+%fx2n',a,b,c);xx=1:0.01:10;yy=cc(1)+cc(2)*xx+cc(3)*xx.2;plot(x,y,'r*',xx,yy);ymin=min(yy);c1=a-ymin;p=c b c1;xmin=roots(p);fprintf('最低点的横纵坐标分别为');end上机运行>> x=1 3 4 5 6 7 8 9 10;>> y=10 5 4 2 1 1 2 3 4;>> SecondFitting(x,y)拟合的曲线为13.-3.x+0.x2最低点的横纵坐标分别为ans = 6.74 6.7342图像如下图所示根据程序运行结果，得到拟合方程为y=13.-3.x+0.x2,最低点坐标为（6.74，6.7342）。计算方法第七次上机作业15.用龙贝格算法计算椭圆x2400+y2100 =1的周长，使误差不超过10-4.解：椭圆周长L计算公式为L=402300(sin)2+100 d上机程序如下所示（注：把每次计算结果存储于一个二维数组之中，为了输出如书中所示的表格形式，每一次结果的下标有一定的调整）function = Romberg( fun,a,b,e )%fun表示被积函数，a为区间下限，b为上限，e为精度要求%龙贝格算法求椭圆周长T=zeros();T(1,1)=(b-a)/2*(feval(fun,a)+feval(fun,b);T(2,1)=0.5*T(1,1)+(b-a)/2*(feval(fun,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(fun,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(fun,a+(2*(1:4)*(b-a)/8);T(4,2)=(4*T(4,1)-T(3,1)/3;T(4,3)=(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(fun,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)/(2(k-1)*(sum(feval(fun,a+(2*(1:(2(k-2)*(b-a)/(2(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;enddisp(4*T);%区间为四分之一圆周，所以乘以4输出为椭圆周长上机运行，结果如下（复制数字结果，排版容易乱，截图如下）运行命令>> fun=(x)power(300*(sin(x).2)+100),0.5);>> Romberg(fun,0,pi/2,0.0001)由输出结果可得，满足精度要求的椭圆的周长计算结果为96.57，大约为96.8846，序列表长宽均为23计算方法第八次上机作业17.用下列方法求dydx=23xy-2,x0,1y0=1的数值解（取h=0.1）,并将计算结果与准确解y=31+x2进行比较:(1) 欧拉方法(2) 改进欧拉方法(3) 经典R-K法解：本题f(x,y)= 23xy-2,h=0.1,上机程序如下function A = Equation(fun,a,b,y0,h)%f表示函数,a,b为区间端%y0为0处的函数值，h为步长%欧拉方法x=a:h:b;n=(b-a)/h;y=0*x;y(1)=y0;A(1,1)=x(1);A(1,2)=y(1);for k=1:n y(k+1)=y(k)+h*feval(fun,x(k),y(k); A(k+1,1)=x(k+1); A(k+1,2)=y(k+1);end%改进欧拉方法¨A(1,3)=y(1);for k=1:n s1=feval(fun,x(k),y(k); ys=y(k)+h*s1; s2=feval(fun,x(k+1),ys); y(k+1)=y(k)+0.5*h*(s1+s2); A(k+1,3)=y(k+1);end%R-K方法¨A(1,4)=y(1);for k=1:n K1=feval(fun,x(k),y(k); K2=feval(fun,x(k)+0.5*h,y(k)+0.5*h*K1); K3=feval(fun,x(k)+0.5*h,y(k)+0.5*h*K2); K4=feval(fun,x(k)+h,y(k)+h*K3); y(k+1)=y(k)+h*(K1+2*K2+2*K3+K4)/6; A(k+1,4)=y(k+1);end%精确解for k=0:n x=a+k*h; y(k+1)=(1+x2)(1/3); A(k+1,5)=y(k+1);end end 上机运行，结果如下>> f=(x,y)2*x*power(y,-2)/3;>> Equation(f,0,1,1,0.1)清晰版从左往右，依次是欧拉方法，改进欧拉方法，R-K方法计算的结果专心-专注-专业