西安交通大学计算方法上机作业
计算方法上机作业
1.对以下和式计算：01421181
84858616n n S n n n n ∞
=??
=--- ?++++??∑
，要求： (1)若只需保留11个有效数字，该如何进行计算； (2)若要保留30个有效数字，则又将如何进行计算；
（1）解题思想和算法实现：
根据保留有效位数的要求，可以由公式
得出计算精度要求。只需要很
少内存，时间复杂度和d 呈线性，不需要高浮点支持。先根据while 语句求出符合精度要求的n 值的大小，然后利用for 语句对这n 项进行求和，输出计算结果及n 值大小即可。
（2）matlab 源程序：
保留11位有效数字时； clear clc
format long n=0;
sum=1/(16^n)*(4/(8*n+1)-2/(8*n+4)-1/(8*n+5)-1/(8*n+6)); while sum>=5*10^(-11); n=n+1;
sum=1/(16^n)*(4/(8*n+1)-2/(8*n+4)-1/(8*n+5)-1/(8*n+6)); end
fori=0:n-1;
sum=sum+1/(16^i)*(4/(8*i+1)-2/(8*i+4)-1/(8*i+5)-1/(8*i+6)); end
vpa(sum,11) n
保留30位有效数字时； clear clc
format long n=0;
sum=1/(16^n)*(4/(8*n+1)-2/(8*n+4)-1/(8*n+5)-1/(8*n+6)); whilesum>=5*10^(-30); n=n+1;
sum=1/(16^n)*(4/(8*n+1)-2/(8*n+4)-1/(8*n+5)-1/(8*n+6)); end
fori=0:n-1;
sum=sum+1/(16^i)*(4/(8*i+1)-2/(8*i+4)-1/(8*i+5)-1/(8*i+6)); end
vpa(sum,30) n
（3）实验结果分析
图1.1 保留11位有效数字的n值及计算结果图图1.2 保留30位有效数字的n值及计算结果图
由计算结果可知，通过合理的误差控制，分别通过7次和22次循环，可以实现题目所要求的精确度。
2.某通信公司在一次施工中，需要在水面宽度为20米的河沟底部沿直线走向铺设一条沟底光缆。在铺设光缆之前需要对沟底的地形进行初步探测，从而估计所需光缆的长度，为工程预算提供依据。已探测
(1)请用合适的曲线拟合所测数据点；
(2)预测所需光缆长度的近似值，并作出铺设河底光缆的曲线图；（1）解题思想和算法原理
给定区间[a, b]一个分划
⊿：a=x0 若函数S(x)满足下列条件： 1)S(x)在每个区间[x i, x j]上是不高于3次的多项式。 2)S(x)及其2阶导数在[a, b]上连续。则称S(x)使关于分划⊿的三次样条函数。（2）matlab源程序： clc,clear x=0:1:20; y=[9.01 8.96 7.96 7.97 8.02 9.05 10.13 11.18 12.26 13.28 13.32 12.61 11.29 10.22 9.15 7.9 7.95 8.86 9.81 10.80 10.93]; n=length(x); l(1)=0;m(n)=0; h=diff(x); df=diff(y)./diff(x); d(1)=0;d(n)=0; for j=2:n-1 l(j)=h(j)/(h(j-1)+h(j)); m(j)=h(j-1)/(h(j-1)+h(j)); d(j)=6*(df(j)-df(j-1))/(h(j-1)+h(j)); end m=m(2:end); u=diag(m,-1);r=diag(l,1);a=diag(2*ones(1,n)); A=u+r+a; M=inv(A)*d'; syms g for j=1:n-1 s(j)=M(j)*(x(j+1)-g)^3/(6*h(j))+M(j+1)*((g-x(j))^3/(6*h(j)))+(y(j)-M(j)*h(j)^2/6)*(x(j+1) -g)/h(j)+(y(j+1)-M(j+1)*h(j)^2/6)*(g-x(j))/h(j); end s r=0; for j=1:n-1 df=diff(s(j),g); warning off all; q=int(sqrt(1+df.^2),g,j-1,j); r=r+q; end L=vpa(r,8); disp('the length of the label is L='); disp(L); for j=1:n-1 S(j,:)=sym2poly(s(j)); end for j=1:n-1 x1=x(j):0.1:x(j+1); y1=polyval(S(j,:),x1); if j==1 y2=y1; else fori=1:11 k=(j-1)*10+i; y2(k)=y1(i); end end end x2=x(1):0.1:x(n); plot(x,y,'o') grid hold on plot(x2,y2,'r') （3）实验结果分析 图2.1 铺设河底电缆长度 图2.2 铺设河底光缆的曲线图 由三次样条插值得出的函数曲线的长度和即铺设河底电缆的长度为 26.498514。为了提高插值精度，用三次样条插值可以增加插值节点的办法来满 足要求，而且在给定节点数的条件下，三次样条插值的精度要优于多项式插值以及线性分段插值，虽然舍弃了降低误差这个优点，但是其曲线的光滑性要好一些。 3.假定某天的气温变化记录如下表所示，试用数据拟合的方法找出这 （1）解题思想和数学原理： 对于具体实验时，通常不是先给出函数的解析式，再进行实验，而是通过实验的观察和测量给出离散的一些点，再来求出具体的函数解析式。又因为测量误差的存在，实际真实的解析式曲线并不一定通过测量给出的所有点。最小二乘法，形成法方程是求解这一问题的很好的方法，本实验运用这一方法实现对给定数据的拟合。 对于给定的测量数据(x i ,f i )(i=1,2,…，n )，设函数分布为 ∑==m j j j x a x y 0 ) ()(? 特别的，取 ) (x j ?为多项式 j j x x =)(? (j=0, 1,…，m ) 则根据最小二乘法原理，可以构造泛函 ∑∑==-=n i m j i j j i m x a f a a a H 1 10) )((),,,(? 令 0=??k a H (k=0, 1,…，m ) 则可以得到法方程 ???? ????? ???=????????????????????????),(),(),(),(),(),(),(),(),(),(),(),(1010101111000100m m m m m m m m f f f a a a ????????????????????? 求该解方程组，则可以得到解 m a a a ,,,10 ，因此可得到数据的最小二乘解 ∑=≈m j j j x a x f 0 ) ()(? （2）matlab 源程序： x=[0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24]; %给出题目数据（时间） y=[15 14 14 14 14 15 16 18 20 20 23 25 28 31 32 31 29 27 25 24 22 20 18 17 16]; %给出题目数据（温度） plot(x,y, 'm*') %画出各个离散数据点 hold on for n=2:4; %2、3、4代表拟合函数最高项次数 alltemp=25; % alltemp 代表数据点总共有25个 A=zeros(n+1,n+1); %定义初始正规方程组的系数矩阵A C=ones(n+1,1); %定义初始正规方程组的系数向量C D=zeros(n+1,1); %定义初始正规方程组的向量D fori=1:n+1 for j=1:n+1 for k1=1: alltemp A(i,j)=A(i,j)+(x(k1).^(i-1+j-1)); end end for k2=1: alltemp D(i,1)=D(i,1)+(x(k2).^(i-1)).*(y(k2)); end end %以上为计算出正规方程组矩阵A 、D 的所有元素的程序 tol=1.0e-12; maxit=1000; C=bicg(A,D,tol,maxit); %使用bicg 迭代算出正规方程组的系数向量C p=0; %误差分量 E=0; %误差总量 if n==2 b=0:24; f=C(1)+C(2).*b+C(3).*(b.^2); p=y(b+1)-f; for v=1:25 E=E+(p(v)).^2; end plot(b,f, 'r-') end %以上是对2阶拟合函数的图形处理与误差计算 if n==3 b=0:24; f=C(1)+C(2).*b+C(3).*(b.^2)+C(4).*(b.^3); p=y(b+1)-f; for v=1:25 E=E+(p(v)).^2; end plot(b,f, 'g-') end %以上是对3阶拟合函数的图形处理与误差计算 if n==4 b=0:24; f=C(1)+C(2).*b+C(3).*(b.^2)+C(4).*(b.^3)+C(5).*(b.^4); p=y(b+1)-f; for v=1:25 E=E+(p(v)).^2; end plot(b,f, 'b-') end %以上是对4阶拟合函数的图形处理与误差计算 C,E end n=2; %重新对n赋值，进行指数函数拟合 A=zeros(n+1,n+1); %重新对A矩阵赋初值 C=zeros(n+1,1); %重新对C向量赋初值 D=zeros(n+1,1); %重新对D向量赋初值 fori=1:n+1 for j=1:n+1 for k=1: alltemp A(i,j)=A(i,j)+(x(k).^(i-1+j-1)); end end for l=1: alltemp D(i,1)=D(i,1)+((x(l).^(i-1)).*(log(y(l)))); end end %计算出A矩阵、D向量各元素数值 C=bicg(A,D,tol,maxit); %利用bicg迭代求解系数 b=0:24; p=0; E=0; f=exp(C(1)+C(2).*b+C(3).*(b.^2)); p=y(b+1)-f; for v=1:25 E=E+(p(v)).^2; end plot(b,f, 'c-') %对指数拟合函数进行图形处理和误差计算 b=-C(3); c=C(2)/(2*b); a=exp(b*(c^2)+C(1)); %算出题设要求的指数拟合函数的各个系数 a,b,c,E grid on legend('测量数据','二次函数','三次函数','四次函数','指数拟合','Location','NorthWest') hold off %hold on与hold off联合使用，表示将各个曲线画在同一个图中 图3.1 二次曲线拟合系数与2范数误差 图3.2 三次曲线拟合系数与2范数误差 图3.3 四次曲线拟合系数与2范数误差 图3.4 指数曲线拟合系数与2范数误差 图3.5 数据原始点与拟合曲线对比图 （3）实验结果分析： 根据国家有关规定，用的是02、08、14、20时四个观测时次的数据做平均，最有代表性。从图中可以看出并不是多项式次数越高越好，随着次数的增高，曲线所呈现出的给定点处和实际的吻合度越好，但对于其他地方的吻合度降低了。 4.设计算法，求出非线性方程52 -+=的所有根，并使误差不 x x 645200 超过4 10-。 解：（1）解题思想和算法实现： 对于一个非线性方程的数值解法很多。在此介绍两种最常见的方法：二分法和Newton 法。首先要研究函数的形态，确定根的数量和大致区间的位置。 对于二分法，其数学实质就是说对于给定的待求解的方程f(x)，其在[a,b ]上连续，f(a)f(b)<0，且f(x)在[a,b ]内仅有一个实根x *，取区间中点c ，若，则c 恰为其根，否则根据f(a)f(c)<0是否成立判断根在区间[a,c ]和[c,b ]中的哪一个，从而得出新区间，仍称为[a,b ]。重复运行计算，直至满足精度为止。这就是二分法的计算思想。 Newton 法通常预先要给出一个猜测初值x 0，然后根据其迭代公式 )() ('1k k k k x f x f x x - =+ 产生逼近解x *的迭代数列{x k }，这就是Newton 法的思想。当x 0接近x *时收敛很快，但是当x 0选择不好时，可能会发散，因此初值的选取很重要。另外，若将该迭代公式改进为 )()('1k k k k x f x f r x x -=+ 其中r 为要求的方程的根的重数，这就是改进的Newton 法，当求解已知重数的方程的根时，在同种条件下其收敛速度要比Newton 法快的多。 本题中用matlab 画曲线52y 64520x x =-+如下： 图4.1 y=6*x.^5-45*x.^2+20的曲线 由曲线可以看出，该方程有三个实根，根所在区间依次为： [-1,-0.5] [0.5,1] [1.5,2] 采用Newton 迭代法，依据迭代式： 1f x x f +=-(k)(k ) (k) '(k)(x ) (x ) ，k=0,1,2,… (8-1) 该方程的迭代式为： 521464520 3090x x x x x x +-+=--(k)(k)(k ) (k) (k)(k) （8-2） 分别选用迭代初值01x =-()、01x =()、02x =()进行迭代，分别求得满足精度的三个实根。 （2）matlab 源程序： clc;clear x=-2:0.1:2; y=6*x.^5-45*x.^2+20; plot(x,y,'g')%画相应的曲线 grid title('y=6*x.^5-45*x.^2+20曲线') format long roots([6 0 0 -45 0 20])%根的真实解 clear -2 -1.5-1-0.500.51 1.52 -400 -350-300-250-200-150-100-500 50y=6*x.5-45*x.2+20曲线 x0=input('请输入迭代初值x0='); format long i=1; e=10^-4;%精度 x=x0-(6*x0.^5-45*x0.^2+20)/(30*x0.^4-90*x0); while abs(x-x0)>e i=i+1; x0=x; x=x0-(6*x0.^5-45*x0.^2+20)/(30*x0.^4-90*x0);%迭代 end display('方程的根') x display('迭代的次数') i （3）实验结果分析： 图4.2 运行结果 对于Newton迭代法，三个初值x0都使得迭代收敛，这是非常重要的。考虑Newton法迭代的收敛性条件： （1）存在一个区间，满足。由曲线和所选的三个区间可知这一条件满足。 （2）f(x)是[a,b]上的单调函数，即对一切不变号。经验证所选的三个区间满足这一条件。 （3）f(x)的凹向在[a,b]上不变，即在[a,b]上不改变符号。经验证所选的三个区间满足这一条件。 它一般是超线性收敛的。 5.编写程序实现大规模方程组的列主元高斯消去法程序，并对所附的方程组进行求解。针对本专业中所碰到的实际问题，提炼一个使用方 程组进行求解的例子，并对求解过程进行分析、求解。 解:（1）算法原理 由于一般线性方程在使用Gauss 消去法求解时，从求解的过程中可以看到， 若)1(-k kk a =0，则必须进行行交换，才能使消去过程进行下去。有的时候即使 ≠-)1(k kk a 0，但是其绝对值非常小，由于机器舍入误差的影响，消去过程也会出现不稳定得现象，导致结果不正确。因此有必要进行列主元技术，以最大可能的消除这种现象。这一技术要寻找行r ，使得 ) 1()1(max ||->-=k ik k i k rk a a 并将第r 行和第k 行的元素进行交换，以使得当前的) 1(-k kk a 的数值比0要大的多。 这种列主元的消去法的主要步骤如下： 1. 消元过程 对k =1,2,…,n -1,进行如下步骤。 1) 选主元，记ik k i rk a a >=max || ||rk a 很小，这说明方程的系数矩阵严重病态，给出警告，提示结果可能不对。 2) 交换增广阵A 的r ，k 两行的元素。 kj rj a a ? (j=k,…,n +1) 3）计算消元 kk kj ik ij ij a a a a a /-=(i=k+1,…,n ; j =k +1,……,n +1) 2. 回代过程 对k = n , n -1,…,1,进行如下计算 ) /(1 1,∑-=+- =n k j kk j kj n k k a x a a x 至此，完成了整个方程组的求解。 （2）matlab 源程序： %非压缩,dat51.dat 、dat53.dat clear;clc fp=fopen('dat53.dat','rb'); id=fread(fp,1,'int32'); ver=fread(fp,1,'int32'); N=fread(fp,1,'int32'); q=fread(fp,1,'int32'); fori=1:N A(i,:)=fread(fp,N,'float'); end fori=1:N d(i)=fread(fp,1,'float'); end %正向消去过程 fori=1:N-q for k=1:p ll=A(i+k,i)/A(i,i); for j=i:i+q A(i+k,j)=A(i+k,j)-ll*A(i,j); end d(i+k)=d(i+k)-ll*d(i); end end fori=N-q+1:N for k=1:N-i ll=A(i+k,i)/A(i,i); for j=i:N A(i+k,j)=A(i+k,j)-ll*A(i,j); end d(i+k)=d(i+k)-ll*d(i); end end %回代过程 x(N)=d(N)/A(N,N); fori=N-1:-1:1 S=0; ifi+q>N cv=N;%cv-critical value else cv=i+q; end for j=i+1:cv S=S+A(i,j)*x(j); end x(i)=(d(i)-S)/A(i,i); end x %压缩,dat52.dat、dat54.dat clear;clc fp=fopen('dat54.dat','rb'); ver=fread(fp,1,'int32'); N=fread(fp,1,'int32'); q=fread(fp,1,'int32'); p=fread(fp,1,'int32'); fori=1:N A(i,:)=fread(fp,p+q+1,'float'); end fori=1:N d(i)=fread(fp,1,'float'); end %正向消去过程 fori=1:N ifi+p cv=p;%cv-critical value else cv=N-i; end for k=1:cv ll=A(i+k,p+1-k)/A(i,p+1); for j=p+1:p+q+1 A(i+k,j-k)=A(i+k,j-k)-ll*A(i,j); end d(i+k)=d(i+k)-ll*d(i); end end %回代过程 x(N)=d(N)/A(N,p+1); fori=N-1:-1:1; S=0; ii=i+1; jj=p+2; while (ii<=N)&&(jj<=p+q+1) S=S+A(i,jj)*x(ii); ii=ii+1; jj=jj+1; end x(i)=(d(i)-S)/A(i,p+1); end x （3）实验结果： 非压缩矩阵求解结果（部分） 压缩矩阵求解结果（部分） （4）分析心得： 采用Gauss消去法时，如果在消元时对角线上的元素始终较大，那么本方法不需要进行列主元计算，计算结果一般就可以达到要求，否则必须进行列主元这一步，以减少机器误差带来的影响，使方法得出的结果正确。 （5）实例 化学反应方程式严格遵守质量守恒定律，书写化学反应方程式写出反应物和 生成物后，往往左右两边各原子数目不相等，不满足质量守恒定律，这就需要通 过计算配平来解决。对于化学反应方程式 X1KMnO4+x2MnSO4+x3H2O=x4MnO2+x5K2SO4+x6H2SO4，可按以下方式配平：上述化学反应式中包含5种不同的原子（钾、锰、氧、硫、氢），按照原子 守恒有: K: x1=2x5 Mn: x1+x2=x4 O: 4x1+4x2+x3=2x4+4x5+4x6 S: x2=x5+x6 H: 3x2=2x6 上述方程组中有5个方程，6个未知量，因此有无数解。可先令x6=1，于是形成方程组。 使用Gauss消去法求解此方程组得: x =[1.0000,1.5000,1.0000,2.5000,0.5000,1.0000]；由于化学方程式通常取最简的正整数，因此将x乘2得配平的化学方程式： 2KMnO4+3MnSO4+2H2O=5MnO2+K2SO4+2H2SO4 程序如下： clear;clc; a=[1 0 0 0 -2 0; 1 1 0 -1 0 0; 4 4 1 -2 -4 -4; 0 1 0 0 -1 -1; 0 0 2 0 0 -2; 0 0 0 0 0 -1]; b=[0 0 0 0 0 -1]'; n=length(b); for k=1:n-1 if a(k,k)==0 disp('error'); return; end fori=k+1:n a(i,k)=a(i,k)/a(k,k); for j=k+1:n a(i,j)=a(i,j)-a(i,k)*a(k,j); end b(i)=b(i)-a(i,k)*b(k); end end x(n)=b(n)/a(n,n); for k=n-1:-1:1 S=b(k); for j=k+1:n S=S-a(k,j)*x(j); end x(k)=S/a(k,k); end x 西安交大统计学考试试卷 一、单项选择题（每小题2 分，共20 分） 1.在企业统计中，下列统计标志中属于数量标志的是（C） A、文化程度 B、职业 C、月工资 D、行业 2.下列属于相对数的综合指标有（B ） A、国民收入 B、人均国民收入 C、国内生产净值 D、设备台数 3.有三个企业的年利润额分别是5000 万元、8000 万元和3900 万元，则这句话中有（B）个变量？ A、0 个 B、两个 C、1 个 D、3 个 4.下列变量中属于连续型变量的是（A ） A、身高 B、产品件数 C、企业人数 D、产品品种 5.下列各项中，属于时点指标的有（A ） A、库存额 B、总收入 C、平均收入 D、人均收入 6.典型调查是（B ）确定调查单位的 A、随机 B、主观 C、随意 D 盲目 7.总体标准差未知时总体均值的假设检验要用到（A ）： A、Z 统计量 B、t 统计量 C、统计量 D、X 统计量 8.把样本总体中全部单位数的集合称为（A ） A、样本 B、小总体 C、样本容量 D、总体容量 9.概率的取值范围是p（D ） A、大于1 B、大于－1 C、小于1 D、在0 与1 之间 10.算术平均数的离差之和等于（A ） A、零 B、1 C、－1 D、2 二、多项选择题（每小题2 分，共10 分。每题全部答对才给分，否则不计分） 1.数据的计量尺度包括（ABCD ）： A、定类尺度 B、定序尺度 C、定距尺度 D、定比尺度 E、测量尺度 2.下列属于连续型变量的有（BE ）： A、工人人数 B、商品销售额 C、商品库存额 D、商品库存量 E、总产值 3.测量变量离中趋势的指标有（ABE ） A、极差 B、平均差 C、几何平均数 D、众数 E、标准差 4.在工业企业的设备调查中（BDE ） A、工业企业是调查对象 B、工业企业的所有设备是调查对象 C、每台设备是 填报单位D、每台设备是调查单位E、每个工业企业是填报单位 5.下列平均数中，容易受数列中极端值影响的平均数有（ABC ） A、算术平均数 B、调和平均数 C、几何平均数 D、中位数 E、众数 三、判断题（在正确答案后写“对”，在错误答案后写“错”。每小题1 分，共10 分） 1、“性别”是品质标志。（对） 数值分析上机实验报告 选题：曲线拟合的最小二乘法 指导老师： 专业： 学号： 姓名： 课题八曲线拟合的最小二乘法 一、问题提出 从随机的数据中找出其规律性，给出其近似表达式的问题，在生产实践和科学实验中大量存在，通常利用数据的最小二乘法求得拟合曲线。 在某冶炼过程中，根据统计数据的含碳量与时间关系，试求含碳量y 与时间t 的拟合曲线。 二、要求 1、用最小二乘法进行曲线拟合； 2、近似解析表达式为()33221t a t a t a t ++=?； 3、打印出拟合函数()t ?，并打印出()j t ?与()j t y 的误差，12,,2,1 =j ； 4、另外选取一个近似表达式，尝试拟合效果的比较； 5、*绘制出曲线拟合图*。 三、目的和意义 1、掌握曲线拟合的最小二乘法； 2、最小二乘法亦可用于解超定线代数方程组； 3、探索拟合函数的选择与拟合精度间的关系。 四、计算公式 对于给定的测量数据(x i ,f i )(i=1,2,…，n )，设函数分布为 ∑==m j j j x a x y 0)()(? 特别的，取)(x j ?为多项式 j j x x =)(? (j=0, 1,…，m ) 则根据最小二乘法原理，可以构造泛函 ∑∑==-=n i m j i j j i m x a f a a a H 1 10))((),,,(? 令 0=??k a H (k=0, 1,…，m ) 则可以得到法方程 ???? ??????? ?=????????????????????????),(),(),(),(),(),(),(),(),(),(),(),(1010101111000100m m m m m m m m f f f a a a ????????????????????? 求该解方程组，则可以得到解m a a a ,,,10 ，因此可得到数据的最小二乘解 ∑=≈m j j j x a x f 0)()(? 曲线拟合：实际工作中，变量间未必都有线性关系，如服药后血药浓度与时间的关系；疾病疗效与疗程长短的关系；毒物剂量与致死率的关系等常呈曲线关系。曲线拟合是指选择适当的曲线类型来拟合观测数据，并用拟合的曲线方程分析两变量间的关系。 五、结构程序设计 在程序结构方面主要是按照顺序结构进行设计，在进行曲线的拟合时，为了进行比较，在程序设计中，直接调用了最小二乘法的拟合函数polyfit ，并且依次调用了plot 、figure 、hold on 函数进行图象的绘制，最后调用了一个绝对值函数abs 用于计算拟合函数与原有数据的误差，进行拟合效果的比较。 1.计算以下和式：01421181 84858616n n S n n n n ∞ =?? =--- ?++++??∑ ，要求： (1)若保留11个有效数字，给出计算结果，并评价计算的算法； (2)若要保留30个有效数字，则又将如何进行计算。 （1）题目分析 该题是对无穷级数求和，因此在使用matlab 进行累加时需要一个累加的终止条件。这里令?? ? ??+-+-+-+= 681581482184161n n n n a n n ，则 ()()1.016 1 6855844864816114851384128698161 681581482184161148113811282984161111<< ? ??? ????? ??++++++???? ????? ??++++++=??? ????? ??+-+-+-+??? ????? ??+-+-+-+=+++n n n n n n n n n n n n n n n n a a n n n n n n 故近似取其误差为1+≈k a ε，并且有m -1m -111021 21 ?=?=≈+βεk a ， （2）算法依据 使用matlab 编程时用digits 函数和vpa 函数来控制位数。 （3）Matlab 运行程序 %%保留11位有效数字 k1=11; s1=0;%用于存储这一步计算值 for n=0:50 a=(1/16^n)*(4/(8*n+1)-2/(8*n+4)-1/(8*n+5)-1/(8*n+6)); n1=n-1; if a<=0.5*10^(1-k1) break end end; for i=0:1:n1 t=(1/16^i)*(4/(8*i+1)-2/(8*i+4)-1/(8*i+5)-1/(8*i+6)); s1=s1+t; end s11=vpa(s1,k1); disp('保留11位有效数字的结果为：');disp(s11); disp('此时n 值为：');disp(n1); %%保留30位有效数字 clear all; k2=30; 计算方法上机报告 姓名： 学号： 班级： 上课班级： 说明： 本次上机实验使用的编程语言是Matlab 语言，编译环境为MATLAB 7.11.0，运行平台为Windows 7。 1. 对以下和式计算： ∑ ∞ ? ?? ??+-+-+-+=0681581482184161n n n n S n ，要求： ① 若只需保留11个有效数字，该如何进行计算； ② 若要保留30个有效数字，则又将如何进行计算； （1） 算法思想 1、根据精度要求估计所加的项数，可以使用后验误差估计，通项为： 1421114 16818485861681 n n n a n n n n n ε??= ---<< ?+++++??； 2、为了保证计算结果的准确性，写程序时，从后向前计算； 3、使用Matlab 时，可以使用以下函数控制位数： digits(位数)或vpa(变量，精度为数) （2）算法结构 1. ;0=s ?? ? ??+-+-+-+= 681581482184161n n n n t n ； 2. for 0,1,2,,n i =??? if 10m t -≤ end; 3. for ,1,2,,0n i i i =--??? ;s s t =+ （3）Matlab源程序 clear; %清除工作空间变量 clc; %清除命令窗口命令 m=input('请输入有效数字的位数m='); %输入有效数字的位数 s=0; for n=0:50 t=(1/16^n)*(4/(8*n+1)-2/(8*n+4)-1/(8*n+5)-1/(8*n+6)); if t<=10^(-m) %判断通项与精度的关系break; end end; fprintf('需要将n值加到n=%d\n',n-1); %需要将n值加到的数值 for i=n-1:-1:0 t=(1/16^i)*(4/(8*i+1)-2/(8*i+4)-1/(8*i+5)-1/(8*i+6)); s=s+t; %求和运算 end s=vpa(s,m) %控制s的精度 （4）结果与分析 当保留11位有效数字时，需要将n值加到n=7， s =3.1415926536； 当保留30位有效数字时，需要将n值加到n=22, s =3.14159265358979323846264338328。 通过上面的实验结果可以看出，通过从后往前计算，这种算法很好的保证了计算结果要求保留的准确数字位数的要求。 第一章 一、题目 设∑ =-= N N j S 2 j 2 1 1，其精确值为)11 123(21+--N N 。 1) 编制按从大到小的顺序1 1 13112122 2-+??+-+-=N S N ，计算S N 的通用程序。 2) 编制按从小到大的顺序1 21 1)1(111222-+ ??+--+-= N N S N ，计算S N 的通用程序。 3) 按两种顺序分别计算64210,10,10S S S ,并指出有效位数。（编制程序时用单精度） 4) 通过本次上机题，你明白了什么？ 二、通用程序 N=input('Please Input an N (N>1):'); AccurateValue=single((0-1/(N+1)-1/N+3/2)/2); Sn1=single(0); for a=2:N; Sn1=Sn1+1/(a^2-1); end Sn2=single(0); for a=2:N; Sn2=Sn2+1/((N-a+2)^2-1); end fprintf('The value of Sn (N=%d)\n',N); fprintf('Accurate Calculation %f\n',AccurateValue); fprintf('Caculate from large to small %f\n',Sn1); fprintf('Caculate from small to large %f\n',Sn2); disp('____________________________________________________') 三、结果 从结果可以看出有效位数是6位。 感想：可以得出，算法对误差的传播有一定的影响，在计算时选一种好的算法可以使结果更为精确。从以上的结果可以看到从大到小的顺序导致大数吃小数的现象，容易产生较大的误差，求和运算从小数到大数所得到的结果才比较准确。 2.用下列方法求方程e^x+10x-2=0的近似根，要求误差不超过5*10的负4次方，并比较计算量 （1）二分法 （局部，大图不太看得清，故后面两小题都用局部截图） （2）迭代法 （3）牛顿法 顺序消元法 #include for(k=p+1;k 数值计算方法I 上机实验考试题（两题任选一题） １．小型火箭初始质量为900千克，其中包括600千克燃料。火箭竖直向上发射时燃料以15千克/秒的速率燃烧掉，由此产生30000牛顿的恒定推力．当燃料用尽时引擎关闭。设火箭上升的整个过程中，空气阻力与速度平方成正比，比例系数为0.4（千克/米）．重力加速度取9.8米/秒2. A. 建立火箭升空过程的数学模型（微分方程）； B. 求引擎关闭瞬间火箭的高度、速度、加速度，及火箭到达最高点的时间和高度． ２．小型火箭初始质量为1200千克，其中包括900千克燃料。火箭竖直向上发射时燃料以15千克/秒的速率燃烧掉，由此产生40000牛顿的恒定推力．当燃料用尽时引擎关闭。设火箭上升的整个过程中，空气阻力与速度平方成正比，比例系数记作k ，火箭升空过程的数学模型为 0)0(,0,01222==≤≤-+?? ? ??-==t dt dx x t t mg T dt dx k dt x d m 其中)(t x 为火箭在时刻t 的高度，m =1200-15t 为火箭在时刻t 的质量，T （=30000牛顿）为推力，g (=9.8米/秒2)为重力加速度, t 1 (=900/15=60秒)为引擎关闭时刻． 今测得一组数据如下（t ~时间（秒），x ~高度（米），v ~速度（米/秒））： 现有两种估计比例系数k 的方法： 1．用每一个数据(t,x,v )计算一个k 的估计值（共11个），再用它们来估计k 。 2．用这组数据拟合一个k ． 请你分别用这两种方法给出k 的估计值，对方法进行评价，并且回答，能否认为空气阻力系数k=0.5（说明理由）． 《数值计算方法》上机实验报告华北电力大学 实验名称数值il?算方法》上机实验课程名称数值计算方法专业班级：电力实08学生姓名：李超然学号:200801001008 成绩: 指导教师:郝育黔老师实验日期：2010年04月华北电力大学实验报告数值计算方法上机实验报吿一. 各算法的算法原理及计算机程序框图1、牛顿法求解非线性方程 *对于非线性方程，若已知根的一个近似值，将在处展开成一阶 xxfx ()0, fx ()xkk 泰勒公式 "f 0 / 2 八八,fxfxfxxxxx 0 0 0 0 0 kkkk2! 忽略高次项，有 ，fxfxfxxx 0 ()()()，，, kkk 右端是直线方程，用这个直线方程来近似非线性方程。将非线性方程的 **根代入，即fx ()0, X ，* fxfxxx 0 0 0 0, ,, kkk fx 0 fx 0 0, 解出 fX 0 *k XX,, k' fx 0 k 水将右端取为，则是比更接近于的近似值，即xxxxk, Ik, Ik fx ()k 八XX, Ikk* fx()k 这就是牛顿迭代公式。 ,2,计算机程序框图:，见, ,3,输入变量、输出变量说明: X输入变量:迭代初值，迭代精度，迭代最大次数，\0 输出变量:当前迭代次数，当前迭代值xkl ,4,具体算例及求解结果: 2/16 华北电力大学实验报吿 开始 读入 l>k /fx()0?,0 fx 0 Oxx，，01* fx ()0 XX,,，?10 kk, ,1，kN, ?xx, 10 输出迭代输出X输出奇异标志1失败标志 ,3,输入变量、输出变量说明: 结束 例：导出计算的牛顿迭代公式，并il ?算。（课本P39例2-16） 115cc （0）, 求解结果: 10. 750000 10.723837 10. 723805 10. 723805 2、列主元素消去法求解线性方程组，1,算法原理: 高斯消去法是利用现行方程组初等变换中的一种变换，即用一个不为零的数乘 -个 方程后加只另一个方程，使方程组变成同解的上三角方程组，然后再自下而上 对上三角 3/16 华北电力大学实验报告方程组求解。 列选主元是当高斯消元到第步时，从列的以下（包括）的各元素中选出绝 aakkkkkk 对值最大的，然后通过行交换将其交换到的位置上。交换系数矩阵中的 两行（包括常ekk 数项），只相当于两个方程的位置交换了，因此，列选主元不影响求解的结 ,2,计算机程序框图:，见下页, 输入变量:系数矩阵元素，常向量元素baiji 输出变量:解向量元素bbb,,12n 计算方法第四次上机报告 2.用欧拉方法解初值 y’=10x(1-y) 0<=x<=1 Y(0)=0 取步长h=0.1,保留5位有效数字，并与准确解相比较 分析：该题目考察欧拉方法解初值问题 程序如下： function Heun(a,b,y0,n) h=(b-a)/n;x=a:h:b; y=y0*ones(1,n+1); for j=2:n+1 yp=y(j-1)+h*f(x(j-1),y(j-1)); yc=y(j-1)+h*f(x(j),yp); y(j)=1/2*(yp+yc); end for k=1:n+1 fprintf('x[%d]=%f\ty[%d]=%f\n',k-1,x(k),k-1,y(k)); end function z=f(xx,yy) z=10*xx*(1-yy); 运行结果： >> Heun(0,1,0,10) x[0]=0.000000 y[0]=0.000000 x[1]=0.100000 y[1]=0.050000 x[2]=0.200000 y[2]=0.183000 x[3]=0.300000 y[3]=0.362740 x[4]=0.400000 y[4]=0.547545 x[5]=0.500000 y[5]=0.705905 x[6]=0.600000 y[6]=0.823543 x[7]=0.700000 y[7]=0.901184 x[8]=0.800000 y[8]=0.947627 x[9]=0.900000 y[9]=0.973290 x[10]=1.000000 y[10]=0.986645 >> 分析： 该结果与准确结果相比比较接近，但是有一定的误差。 6.用四阶龙格—库塔公式解第三题中的初值问题，取步长h=0.2，保留五位有效数字。 题目目的分析： 该题考查四阶龙格-库塔方法和改进欧拉方法求解精确度问题。 程序： 改进欧拉法： function Heun(a,b,y0,n) h=(b-a)/n;x=a:h:b; y=y0*ones(1,n+1); for j=2:n+1 yp=y(j-1)+h*f(x(j-1),y(j-1)); yc=y(j-1)+h*f(x(j),yp); y(j)=1/2*(yp+yc); end for k=1:n+1 fprintf('x[%d]=%f\ty[%d]=%f\n',k-1,x(k),k-1,y(k)); end 计算方法 一、填空题 1.假定x ≤1，用泰勒多项式?+??+++=! !212n x x x e n x ,计算e x 的值，若要求截断误差不超过0.005，则n=_5___ 2. 解 方 程 03432 3=-+x － x x 的牛顿迭代公式 )463/()343(121121311+--+--=------k k k k k k k x x x x x x x 3.一阶常微分方程初值问题 ?????= ='y x y y x f y 0 0)() ,(,其改进的欧拉方法格式为)],(),([21 1 1 y x y x y y i i i i i i f f h +++++= 4.解三对角线方程组的计算方法称为追赶法或回代法 5. 数值求解初值问题的四阶龙格——库塔公式的局部截断误差为o(h 5 ) 6.在ALGOL 中，简单算术表达式y x 3 + 的写法为x+y ↑3 7.循环语句分为离散型循环,步长型循环,当型循环. 8.函数)(x f 在[a,b]上的一次（线性）插值函数= )(x l )()(b f a b a x a f b a b x --+-- 9.在实际进行插值时插值时，将插值范围分为若干段，然后在每个分段上使用低阶插值————如线性插值和抛物插值，这就是所谓分段插值法 10、数值计算中，误差主要来源于模型误差、观测误差、截断误差和舍入误差。 11、电子计算机的结构大体上可分为输入设备 、 存储器、运算器、控制器、 输出设备 五个主要部分。 12、算式2 cos sin 2x x x +在ALGOL 中写为))2cos()(sin(2↑+↑x x x 。 13、ALGOL 算法语言的基本符号分为 字母 、 数字 、 逻辑值、 定义符四大 计算方法实验报告 班级： 学号： 姓名： 成绩： 1 舍入误差及稳定性 一、实验目的 （1）通过上机编程，复习巩固以前所学程序设计语言及上机操作指令； （2）通过上机计算，了解舍入误差所引起的数值不稳定性 二、实验内容 1、用两种不同的顺序计算10000 21n n -=∑，分析其误差的变化 2、已知连分数() 1 01223//(.../)n n a f b b a b a a b =+ +++，利用下面的算法计算f ： 1 1 ,i n n i i i a d b d b d ++==+ (1,2,...,0 i n n =-- 0f d = 写一程序，读入011,,,...,,,...,,n n n b b b a a 计算并打印f 3、给出一个有效的算法和一个无效的算法计算积分 1 041 n n x y dx x =+? (0,1,...,1 n = 4、设2 2 11N N j S j == -∑ ，已知其精确值为1311221N N ?? -- ?+?? （1）编制按从大到小的顺序计算N S 的程序 （2）编制按从小到大的顺序计算N S 的程序 （3）按两种顺序分别计算10001000030000,,,S S S 并指出有效位数 三、实验步骤、程序设计、实验结果及分析 1、用两种不同的顺序计算10000 2 1n n -=∑，分析其误差的变化 （1）实验步骤： 分别从1~10000和从10000~1两种顺序进行计算，应包含的头文件有stdio.h 和math.h （2）程序设计： a.顺序计算 #include 实验报告一 题目：非线性方程求解 摘要：非线性方程的解析解通常很难给出，因此线性方程的数值解法就尤为重要。本实验采用两种常见的求解方法二分法和Newton法及改进的Newton法。 前言：（目的和意义） 掌握二分法与Newton法的基本原理和应用。 数学原理： 对于一个非线性方程的数值解法很多。在此介绍两种最常见的方法：二分法和Newton法。 对于二分法，其数学实质就是说对于给定的待求解的方程f(x)，其在[a,b]上连续，f(a)f(b)<0，且f(x)在[a,b]内仅有一个实根x*，取区间中点c，若，则c恰为其根，否则根据f(a)f(c)<0是否成立判断根在区间[a,c]和[c,b]中的哪一个，从而得出新区间，仍称为[a,b]。重复运行计算，直至满足精度为止。这就是二分法的计算思想。 Newton法通常预先要给出一个猜测初值x0，然后根据其迭代公式 产生逼近解x*的迭代数列{x k}，这就是Newton法的思想。当x0接近x*时收敛很快，但是当x0选择不好时，可能会发散，因此初值的选取很重要。另外，若将该迭代公式改进为 其中r为要求的方程的根的重数，这就是改进的Newton法，当求解已知重数的方程的根时，在同种条件下其收敛速度要比Newton法快的多。 程序设计： 本实验采用Matlab的M文件编写。其中待求解的方程写成function的方式，如下 function y=f(x); y=-x*x-sin(x); 写成如上形式即可，下面给出主程序。 二分法源程序： clear %%%给定求解区间 b=1.5; a=0; %%%误差 R=1; k=0;%迭代次数初值 while (R>5e-6) ; c=(a+b)/2; if f12(a)*f12(c)>0; a=c; else b=c; end R=b-a;%求出误差 k=k+1; end x=c%给出解 Newton法及改进的Newton法源程序：clear %%%% 输入函数 f=input('请输入需要求解函数>>','s') %%%求解f(x)的导数 df=diff(f); 1． 设?+=1 05dx x x I n n ， （1） 由递推公式n I I n n 1 51+-=-，从0I 的几个近似值出发，计算20I ； 解：易得：0I =ln6-ln5=0.1823, 程序为： I=0.182; for n=1:20 I=(-5)*I+1/n; end I 输出结果为：20I = -3.0666e+010 （2） 粗糙估计20I ，用n I I n n 51 5111+- =--，计算0I ； 因为 0095.05 6 0079.01020 201 020 ≈<<≈??dx x I dx x 所以取0087.0)0095.00079.0(2 1 20=+= I 程序为：I=0.0087; for n=1:20 I=(-1/5)*I+1/(5*n); end I 0I = 0.0083 （3） 分析结果的可靠性及产生此现象的原因(重点分析原因)。 首先分析两种递推式的误差；设第一递推式中开始时的误差为000I I E '-=，递推过程的舍入误差不计。并记n n n I I E '-=，则有01)5(5E E E n n n -==-=-Λ。因为=20E 20020)5(I E >>-，所此递推式不可靠。而在第二种递推式中n n E E E )5 1(5110-==-=Λ，误差在缩小， 所以此递推式是可靠的。出现以上运行结果的主要原因是在构造递推式过程中，考虑误差是否得到控制， 即算法是否数值稳定。 2． 求方程0210=-+x e x 的近似根，要求4 1105-+?<-k k x x ，并比较计算量。 （1） 在[0，1]上用二分法； 程序：a=0;b=1.0; while abs(b-a)>5*1e-4 c=(b+a)/2; 数值计算方法上机实习题 1． 设?+=1 05dx x x I n n ， （1） 由递推公式n I I n n 1 51+ -=-，从I 0=0.1824, 0=0.1823I 出发，计算20I ； （2） 20=0I ，20=10000I , 用n I I n n 51 5111+- =--，计算0I ； （3） 分析结果的可靠性及产生此现象的原因(重点分析原因)。 答：第一个算法可得出 e 0=|I 0?I 0 ?| e n =|I n ?I n ?|=5n |e 0| 易知第一个算法每一步计算都把误差放大了5倍，n 次计算后更是放大了5n 倍，可靠性低。 第二个算法可得出 e n =|I n ?I n ?| e 0=(15 )n |e n | 可以看出第二个算法每一步计算就把误差缩小5倍，n 次后缩小了5n 倍，可靠性高。 2． 求方程0210=-+x e x 的近似根，要求41105-+?<-k k x x ，并比较计算量。 （1） 在[0，1]上用二分法； 计算根与步数程序： fplot(@(x) exp(x)+10*x-2,[0,1]); grid on; syms x; f=exp(x)+10*x-2; [root,n]=EFF3(f,0,1); fprintf('root=%6.8f ,n=%d \n',root,n); 计算结果显示： root=0.09057617 ,n=11 （2） 取初值00=x ，并用迭代10 21 x k e x -=+； （3） 加速迭代的结果; （4） 取初值00 x ，并用牛顿迭代法； 计算方法（B）实验报告 姓名： 学号： 学院： 专业： 实验一 三对角方程组Tx f =的求解 一、 实验目的 掌握三对角方程组Tx f =求解的方法。 二、 实验内容 求三对角方程组Tx f =的解，其中： 4 -1 -1 4 -1 -1 4 1 -1 4T ????????=?? ?? ???? ， 3223f ?? ? ? ?= ? ? ??? 三、 算法组织 设系数矩阵为三对角矩阵 11222333111 b c a b c a b c a b c b n n n n T ---???????? =?????? ?????? 则方程组Tx f =称为三对角方程组。 设矩阵T 非奇异，T 可分解为T=LU ，其中L 为下三角矩阵，U 为单位上三角矩阵，记 1 1 212 313 1 1 1111 ,11n n n n n r l r l r L U l r l μμμμμ---???? ? ? ? ? ? ?== ? ? ? ? ? ? ? ? ? ?? ? ? ? 可先依次求出,L U 中的元素后，令Ux y =，先求解下三角方程组Ly f =得出 y ，再求解上三角方程组Ux y =。 追赶法的算法组织如下： 1.输入三对角矩阵T 和右端向量f ； 2.将Tx f =压缩为四个一维数组{}{}{}{}i i i i a b c d 、、、，{}{}{}i i i a b c 、、是T 的三对角线性方程组的三个对角，{}i d 是右端向量。将分解矩阵压缩为三个一维数组 {}{}{}i i i l r μ、、。 3.对T 做Crout 分解（也可以用Doolittle 分解）导出追赶法的计算步骤如下： 1111,b r c μ== for 2i n = 111, , ,i i i i i i i i i i i i i l a b a r r c y d l y μμ---==-==- end 4.回代求解x /n n n x y μ= for 11i n =- 1()/i i i i i x y c x μ+=- end 5. 停止，输出结果。 四、 MATLAB 程序 MATLAB 程序见附件1. 五、 结果及分析 实验结果为： (1.0000 1.0000 1.0000 1.0000)T x = 第一次&第二次上机作业 上机作业： 1.在Matlab上执行：>> 5.1-5-0.1和>> 1.5-1-0.5 给出执行结果，并简要分析一下产生现象的原因。 解：执行结果如下： 在Matlab中，小数值很难用二进制进行描述。由于计算精度的影响，相近两数相减会出现误差。 2.（课本181页第一题） 解：（1）n=0时，积分得I0=ln6-ln5，编写如下图代码 从以上代码显示的结果可以看出，I 20的近似值为0.7465 （2）I I =∫I I 5+I 10dx,可得∫I I 610dx ≤∫I I 5+I 10dx ≤∫I I 510dx,得 16(I +1)≤I I ≤15(I +1),则有1126≤I 20≤1105, 取I 20=1 105 ,以此逆序估算I 0。代码段及结果如下图所示 （3）从I20估计的过程更为可靠。首先根据积分得表达式是可知，被积函数随着n的增大，其所围面积应当是逐步减小的，即积分值应是随着n的递增二单调减小的，（1）中输出的值不满足这一条件，（2）满足。设I I表示I I的近似值，I I-I I=(?5)I(I0?I0)(根据递推公式可以导出此式),可以看出，随着n的增大，误差也在增大，所以顺序估计时，算法不稳定性逐渐增大，逆序估计情况则刚好相反，误差不断减小，算法逐渐趋于稳定。 2.（课本181页第二题） （1）上机代码如图所示 求得近似根为0.09058 （2）上机代码如图所示 得近似根为0.09064； （3）牛顿法上机代码如下 计算所得近似解为0.09091 第三次上机作业上机作业181页第四题 线性方程组为 [1.13483.8326 0.53011.7875 1.16513.4017 2.53301.5435 3.4129 4.9317 1.23714.9998 8.76431.3142 10.67210.0147 ][ I1 I2 I3 I4 ]=[ 9.5342 6.3941 18.4231 16.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); 计算方法A 上机大作业 1. 共轭梯度法求解线性方程组 算法原理：由定理3.4.1可知系数矩阵A 是对称正定矩阵的线性方程组Ax=b 的解与求解二次函数1()2 T T f x x Ax b x =-极小点具有等价性，所以可以利用共轭梯度法求解1()2 T T f x x Ax b x = -的极小点来达到求解Ax=b 的目的。 共轭梯度法在形式上具有迭代法的特征，在给定初始值情况下，根据迭代公式： (1)()()k k k k x x d α+=+ 产生的迭代序列(1)(2)(3)x x x ，，，... 在无舍入误差假定下，最多经过n 次迭代，就可求得()f x 的最小值，也就是方程Ax=b 的解。 首先导出最佳步长k α的计算式。 假设迭代点()k x 和搜索方向()k d 已经给定，便可以通过()()()() k k f x d φαα=+的极小化 ()()min ()()k k f x d φαα=+ 来求得，根据多元复合函数的求导法则得: ()()()'()()k k T k f x d d φαα=?+ 令'()0φα=，得到: ()() ()()k T k k k T k r d d Ad α=，其中()()k k r b Ax =- 然后确定搜索方向()k d 。给定初始向量(0)x 后，由于负梯度方向是函数下降最快的方向，故第一次迭代取搜索方向(0) (0)(0)(0)()d r f x b Ax ==-?=-。令 (1)(0)00x x d α=+ 其中(0)(0)0(0)(0) T T r d d Ad α=。第二次迭代时，从(1) x 出发的搜索方向不再取(1)r ，而是选取(1) (1)(0)0d r d β=+，使得(1)d 与(0)d 是关于矩阵A 的共轭向量，由此可 求得参数0β: 数值分析上机题(matlab版)(东南大学) 数值分析上机报告 第一章 一、题目 精确值为)1 1 123(21+--N N 。 1) 编制按从大到小的顺序 1 1 131121222-+??+-+-= N S N ，计算S N 的通用程序。 2) 编制按从小到大的顺序 1 21 1)1(111222-+??+--+-= N N S N ，计算S N 的通用程序。 3) 按两种顺序分别计算6 42 10,10, 10S S S ,并指出有效位 数。（编制程序时用单精度） 4) 通过本次上机题，你明白了什么？ 二、通用程序 clear N=input('Please Input an N (N>1):'); AccurateValue=single((0-1/(N+1)-1/N+3/2)/2); Sn1=single(0); for a=2:N; Sn1=Sn1+1/(a^2-1); end Sn2=single(0); for a=2:N; Sn2=Sn2+1/((N-a+2)^2-1); end fprintf('The value of Sn using different algorithms (N=%d)\n',N); disp('____________________________________________________') fprintf('Accurate Calculation %f\n',AccurateValue); fprintf('Caculate from large to small %f\n',Sn1); fprintf('Caculate from small to large %f\n',Sn2); 《计算方法》上机指导书 实验1 MATLAB 基本命令 1．掌握MATLAB 的程序设计 实验内容：对以下问题，编写M 文件。 (1) 生成一个5×5矩阵，编程求其最大值及其所处的位置。 (2) 编程求∑=20 1!n n 。 (3) 一球从100米高度自由落下，每次落地后反跳回原高度的一半，再落下。求它在 第10次落地时，共经过多少米？第10次反弹有多高？ 2．掌握MATLAB 的绘图命令 实验内容：对于自变量x 的取值属于[0，3π]，在同一图形窗口画出如下图形。 （1）1sin()cos()y x x =?； （2）21 2sin()cos()3 y x x =-； 实验2 插值方法与数值积分 1. 研究人口数据的插值与预测 实验内容：下表给出了从1940年到1990年的美国人口，用插值方法推测1930年、1965年、2010年人口的近似值。 美国人口数据 1930年美国的人口大约是123,203千人，你认为你得到的1965年和2010年的人口数字精确度如何？ 2．最小二乘法拟合经验公式 实验内容：某类疾病发病率为y ‰和年龄段x (每五年为一段，例如0~5岁为第一段，6~10岁为第二段……)之间有形如bx ae y =的经验关系，观测得到的数据表如下 （1）用最小二乘法确定模型bx ae y =中的参数a 和b 。 （2）利用MATLAB 画出离散数据及拟合函数bx ae y =图形。 3.复化求积公式 实验内容：对于定积分? +=1 02 4dx x x I 。 （1）分别取利用复化梯形公式计算，并与真值比较。再画出计算误差与n 之间的曲线。 （2）取[0，1]上的9个点，分别用复化梯形公式和复化辛普森公式计算，并比较精度。 实验一 矩阵的分解 一、实验目的 掌握矩阵的分解原理和一般方法，学会利用矩阵分解直接求解线性方程组。 二、实验内容 求矩阵() 2020 =ij A α?的T LDL 分解与Cholesky 分解，其中 ,min(,),ij i i j i j i j α=?=? ≠? 。 三、问题分析 1． Cholesky 分解 Cholesky 分解是针对被分解矩阵为对称正定的情况给出的。 分解步骤如下： 11g =1111/y b g =，1111i i g g α= 2i n = ； DO 2j n = jj g = IF 0jj g < STOP ，JUMP TO (5) DO 1i j n =+ 1 1j ij ik kj k ij jj g g g g α-=??- ? ? ?=∑ ji ij g g = 1 1j i ik k k i jj b g y y g -=??- ? ? ?=∑ END DO END DO 2． T LDL 分解 T LDL 分解是针对Cholesky 分解中的开平方运算进行的改进。 分解步骤如下： 11i i r α=，1111/i i r r r =，11y b = 1i n = DO 2i n = DO j i n = 1 1i ij ij ik kj k r l r α-=??=- ??? ∑ /ji ij ii l r r = 1 1i i i ik k k y b l b -=??=- ??? ∑ END DO END DO 四、matlab 求解 分别写出T LDL 分解和Cholesky 分解的函数程序gaijinsqrt.m 和.cholesky m ，调用格 式如下： 1. [index,x,r]=gaijinsqrt(A,b) 参数说明： A 和b 分别是线性代数方程组Ax =b 的系数矩阵和右端向量；输出x 为解向量。 [index,x,g]=Cholesky(A,b) 参数说明： A 和b 分别是线性代数方程组Ax =b 的系数矩阵和右端向量；输出x 为解向量。 然后写出主程序2homework .m 如下： %生成矩阵A A=zeros(20,20); for i=1:20 for j=1:20 if i~=j if i>j A(i,j)=j; else A(i,j)=i; end统计西安交大期末考试试题(含答案)
数值分析上机作业
西安交通大学计算方法B上机试题
计算方法上机作业
数值分析上机题目详解
计算方法上机题答案
数值计算方法I上机实验考试题
《数值计算方法》上机实验报告
计算方法上机作业
计算方法试题库讲解
计算方法上机实习题大作业(实验报告).
(完整版)哈工大-数值分析上机实验报告
(完整版)数值计算方法上机实习题答案
数值计算方法上机实习题
西交大计算方法上机报告
计算方法上机作业集合
西交计算方法A上机大作业
数值分析上机题(matlab版)(东南大学)
计算方法作业2
西安交通大学计算方法A实验报告