常州大学数值分析作业(共六章)
《数值分析》作业参考答案-2
《数值分析》作业参考答案一. 选择题1． A ； 2.B ； 3.B ； 4.D; 5.C; 6.D; 7.C; 8.B; 9.D; 10.C ； 11.B ； 12.A ； 13.A; 14.C; 15.A; 16.B; 17.D; 18.A 19.D,20.C,21.A,22.D,23.C,24.C. 二. 填空题1. 3，3，3 ；2. 1，2/3 ；3. 100！2^100 ；4. (-1 ,1)；5.)())((102010n x x x x x x ---Λ ； 6. xxx x x g sin 1cos )(+--=，2；7. (-1, 1)； 8. x ； 9. 4 ； 10. 5，9 ； 11. )(211nn n x cx x +=+, 2； 12. 31x x ++ ； 13. 10/9, 4； 14. 10, 55, 550; 15. 3ab b -+ 16. 312-x 17. 3118．431,431,21i i +-- 19.x x +22 20. 3b a a -+ . 三.1. 12)(2++=x x x p 2. )()(x f x p = 3.12292.512,916,910====a B C A , 代数精确度为5 3.证明：||)'(||||'||)'(1-⋅=A A A A A A cond设}'m ax {的特征值的模A A =λ，})'m ax {(11的特征值的模--=A A β，则 上式=2212))((||||||||A cond AA =⋅=⋅-βλ4.1. (12分)⎪⎪⎪⎪⎪⎭⎫⎝⎛---=1111433221L ,⎪⎪⎪⎪⎪⎭⎫ ⎝⎛---=4534231112U ,⎪⎪⎪⎪⎪⎭⎫ ⎝⎛-=0101X 2. (8分)Seidel 收敛，因为A 实正定对称阵. 迭代格式⎪⎪⎩⎪⎪⎨⎧+=++-=+=+=+++++++2/)1(2/)2(2/)(2/)2()1(3)1(4)(4)1(2)1(3)(3)1(1)1(2)(2)1(1k k k k k k k k k k x x x x x x x x x x 5. 14)(22+-=x x p π，余项6*54|)2(|61)(cos 322ππ≤-≤-x x x p x 6.6. 证明：当0=A 时，结论显然成立；当0≠A 时，因222||||||||||||||X A Y AX Y T<，故 2220,0||||||||||||||sup A Y X AX Y T Y X ≤≠≠；又A A T是实对称矩阵，故存在正交阵),,,(21n p p p P Λ=使得⎪⎪⎪⎭⎫ ⎝⎛==n TT D AP A P λλO 1， i P 是特征值i λ对应的特征向量。
常州大学数值分析07-08试卷A及参考答案
江苏工业学院2007~2008学年第 2 学期硕士生考试试题参考解答一、（10分）叙述防止误差的几个基本原则，并举一例说明其在数值计算中的应用。
答：防止误差的几个基本原则主要有： 1） 防止大数“吃”小数；2） 避免除数绝对值远远小于被除数绝对值的除法； 3） 避免相近数相减；4） 避免使用不稳定的算法；5） 注意简化计算步骤，减少运算次数； ………… 5 分 例如：当x 充分大时，即1x >>时，计算可以用表达式来计算，以避免相近数相减。
………… 5 分二、(15分)（1）叙述Lagrange 插值或Newton 插值方法的方法思想。
(2) 设(1)0,(2)3,(3)10f f f ===, 试求)(x f 的二次Newton 插值多项式。
解：（1）拉格朗日插值、牛顿插值的方法思想分别如下： 对于给定的节点(,),0,1,2,,i i x y i n = 拉格朗日插值通过引入满足如下条件的基函数1,(),0,i j j il x j i=⎧=⎨≠⎩ 构造如下形式的插值多项式()()nn i i i P x l x y ==∑其中0()()()nj i j ijj ix x l x x x =≠-=-∏。
………… 4 分牛顿插值方法是通过构造如下形式的多项式01020101()()()()()()n n n N x a a x x a x x x x a x x x x -=+-+--++--其中,0,1,2,,i a i n =通过Newton 差商公式得到，且仅与0,1,,,i x x x 有关，由此可以保证在增加节点时， 原先的计算量能够被充分利用。
………… 4 分 (2) 根据列表函数可得差商表如下：0 0 0 3 3 0 10 7 2)(x f 的二次Newton 插值多项式为()3(1)2(1)(2)P x x x x =-+--即2()231P x x x =-+。
数值分析课后习题及答案
第一章 绪论（12） 第二章 插值法（40-42）2、当2,1,1-=x 时，4,3,0)(-=x f ，求)(x f 的二次插值多项式。
[解]372365)1(34)23(21)12)(12()1)(1(4)21)(11()2)(1()3()21)(11()2)(1(0))(())(())(())(())(())(()(2221202102210120120102102-+=-++--=+-+-⨯+------⨯-+-+-+⨯=----+----+----=x x x x x x x x x x x x x x x x x x x y x x x x x x x x y x x x x x x x x y x L 。
3、给出x x f ln )(=的数值表用线性插值及二次插值计算54.0ln 的近似值。
X 0.4 0.5 0.6 0.7 0.8 x ln -0.916291 -0.693147 -0.510826 -0.357765 -0.223144[解]若取5.00=x ，6.01=x ，则693147.0)5.0()(00-===f x f y ，510826.0)6.0()(11-===f x f y ，则604752.182321.1)5.0(10826.5)6.0(93147.65.06.05.0510826.06.05.06.0693147.0)(010110101-=---=--⨯---⨯-=--+--=x x x x x x x x x y x x x x y x L ，从而6202186.0604752.19845334.0604752.154.082321.1)54.0(1-=-=-⨯=L 。
若取4.00=x ，5.01=x ，6.02=x ，则916291.0)4.0()(00-===f x f y ，693147.0)5.0()(11-===f x f y ，510826.0)6.0()(22-===f x f y ，则 217097.2068475.404115.2)2.09.0(5413.25)24.0(3147.69)3.01.1(81455.45)5.06.0)(4.06.0()5.0)(4.0()510826.0()6.05.0)(4.05.0()6.0)(4.0()693147.0()6.04.0)(5.04.0()6.0)(5.0(916291.0))(())(())(())(())(())(()(22221202102210120120102102-+-=+--+-⨯++-⨯-=----⨯-+----⨯-+----⨯-=----+----+----=x x x x x x x x x x x x x x x x x x x x x x y x x x x x x x x y x x x x x x x x y x L ，从而61531984.0217097.21969765.259519934.0217097.254.0068475.454.004115.2)54.0(22-=-+-=-⨯+⨯-=L补充题：1、令00=x ，11=x ，写出x e x y -=)(的一次插值多项式)(1x L ，并估计插值余项。
常州大学数值分析作业
常州大学数值分析作业1.解：（1）x = [ 3*π/8 π/2];Y = cos(x);x0 = π/3;[A,Y] = lagrange(x,y,x0);P1 = vpa(poly2sym(A),3)结果如下：P1 = 1.53*x - 0.974Y = 0.5102（2）x = [π/4 3*π/8π/2];Y = cos(x);[A,Y] = lagrange(x,y,x0);P2=vpa(poly2sym(A),3)结果如下：P2 = 1.18*x^2 - 0.455*x - 0.189Y = 0.4973（3）x = [0 π/8 π/4 3*π/8 π/2];Y = cos(x);[A,Y]=lagrange(x,y,x0);结果如下：P3 = x^4 + 0.00282*x^3 - 0.514*x^2 + 0.0232*x + 0.0287 Y = 0.50017.function [T]=aitken(x,y,x0,T0)If nargin == 3T0=[];endn0=size(T0,1);m=max(size(x));n=n0+m;T=zeros(n,n+1);T(1:n0,1:n0+1)=T0;T(n0+1:n,1)=x;T(n0+1:n, 2)=y; ifn0==0i0=2;elsei0=n0+1;Endx=[0 1];y=[0.5 1.25];x0=2.8;T0=aitken(x,y,x0);T=T0;x=[3.0,4.0]';y=[3.5,2.75]';x0=2.8;T=aitken(x,y,x0,T0);n=max(size(x))+size(T0,1);for i=1:nfor j=1:i+1fprintf('%10.4f',T(i,j));endfprintf('\n');EndReturn0.0000 0.5000 0 0 01.0000 1.25002.6000 0 03.0000 3.5000 3.3000 3.2300 04.0000 2.7500 2.0750 2.2850 3.419016.function [C,D,Y]=newpoly(x0,y0,x)if nargin < 2 | nargin> 3error( 'Incorrect Number of Inputs'); endif length(x0)~=length(y0)error('The length of x0 must be equal to it of y0');end n=length(x0); D=zeros(n,n); D(:,1)=y0'; for j=2:n%计算差商表for k=j:nIf abs(x0(k)-x0(k-j+1))<eps< bdsfid="127" p=""></eps<> error('DividedbyZero,therearetwonodesarethes ame');endD(k,j)=(D(k,j-1)-D(k-1,j-1))/(x0(k)-x0(k-j+1));EndEndC=D(n,n);For k=(n-1):-1:1C=conv(C,poly(x0(k)));m=length(C);C(m)=C(m)+D(k,k);endIf nargin==3Y=polyval(C,x);endC=fliplr(C);Returnx = [0 1 2 3 4 ];y = [0.5,1.25,2.75,3.5,2.75];[A,Y]=lagrange(x,y,x0)x0 = [0 1 2 3 4 ];y0 = [0.5,1.25,2.75,3.5,2.75];[C,D,X]=newpoly(x0,y0,x)plot(x,Y,'b-',x0,X,'r:')A = 0.5000 -0.3125 1.4687 -0.4375 0.0313Y = 0.5000 1.2500 2.75003.5000 2.7500C = 0.0313 -0.4375 1.4688 -0.3125 0.5000D = 0.5000 0 00 01.2500 0.7500 0 0 02.7500 1.5000 0.3750 0 03.5000 0.7500 -0.3750 -0.2500 02.7500 -0.7500 -0.7500 -0.1250 0.0313X = 0.5000 1.2500 2.7500 3.5000 2.7500fl(x)=0.5*x^4 - 0.312*x^3 + 1.47*x^2 - 0.438*x + 0.0312fn(x)=0.5*x^4 - 0.312*x^3 + 1.47*x^2 - 0.438*x + 0.03126. 解：对题中函数进行变形：原式→y/x = a* exp(b*x) →ln(y/x) = ln(a) + b*exp(x) 化为线性形式计算：>> a = [1 2 3 4 5];>> b = [1.222 2.984 5.466 8.902 13.592];>> x = exp(a);>> y = log(b)-log(a);>> n = 1; >> [C]=lspoly(x,y,n);>> y = vpa(poly2sym(C),3)结果如下：y = 0.00464*x + 0.384写成题中拟合函数的形式即为：y = 1.4679*x*exp(0.00464*x)7.function [a0,a1,a2]=ec2(h,w)S=log(s)';N=length(h);A=zeros(N,3);for i=1:5A(i,1)=1;A(i,2)=log(h(i));A(i,3)=log(w(i));endc=inv(A'*A)*(A'*S); a0=exp(c(1)); a1=c(2); a2=c(3);return%给出数据h=[175 172 183 164 156]; w=[80 90 80 70 65];s=[1000 900 1200 750 800]; [a0, a1,a2]=ec2(h,w,s)输出结果为：a0 =1.614815742043648e-04a1 =3.383163094165866a2 =-0.4191650115826638.x= lsqcurvefit(fun,x0,xdata,ydata)x=lsqcurvefit(fun,x0,xdata,ydata,lb,ub)x=lsqcurvefit(fun,x0,xdata,ydata,lb, ub,options)[x,resnorm] = lsqcurvefit(…)[x,resnorm,residual]=lsqcurvefit(…)[x,resnorm,residual,exitflag]= lsqcu rvefit(…)[x,resnorm,residual,exitflag,output] = lsqcurvefit(…)[x,resnorm,residual,exitflag,output,l ambda] = lsqcurvefit(…) [x,resnorm,residual,exitflag,output,l ambda,jacobian] =lsqcurvefit(…)function F = myfun(x,xdata)F=(x(1).*xdata).*(exp(x(2).*xdat a));xdata=[1,2,3,4,5];ydata=[1.222,2.984,5.466,8.902,13. 592]; x0=[0,0];[x,resnorm]=lsqcurvefit(@myfun,x0, xdata,ydata)输出结果为：Local minimum found.Optimization completed because t he size of the gradient is less than t he default value of the function toler ance.x =0.999958348976391 0.2000141328 12834aaresnorm = 8.067930437509675e -7。
数值分析作业(完整版)
的逆阵 A ，用左除命令 A \ E 检验你的结果。
clc clear close all A=[1 1 1 1 1;1 2 3 4 5;1 3 6 10 15;1 4 10 20 35;1 5 15 35 70]; fprintf('对上述矩阵进行列主元素分解：\n') for i=1:1:r-1 [mx,ro]=max(abs(A(i:r,i))); % 寻找a阵第i列的最大值 [A(i,:),A(ro+i-1,:)]=exchange(A(i,:),A(ro+i-1,:)); % 进行行与行交换 for j=i+1:1:r A(j,:)=A(j,:)-A(j,i)/A(i,i)*A(i,:); end A End %--矩阵A的逆阵 A1=inv(A) %--左除验证 E=eye(5); A2=A\E % 5x5单位阵 % A阵的逆矩阵 % 输出每次交换后的A
第一章
1、计算积分 I n
Code： clc clear close all n=9; %--梯形积分法 x=0:0.01:1; y=(x.^n).*exp(x-1); In = trapz(x,y); In2=vpa(In,6) % 6位有效数字 %--高精度积分法 F = @(x1)(x1.^n).*exp(x1-1); s = quad(F,0,1); s1=vpa(s,6)
0
0, 0, 0, 0, 0 。
T
if abs(er(:,i-1))<=e fprintf('在迭代 %d 次之后，满足精度要求，x向量的值如下：\n',i); fprintf('x1=%.5f, x2=%.5f, x3=%.5f, x4=%.5f, x5=%.5f\n',x(1,i),x(2,i),x(3,i),x(4,i),x(5,i)); break end end %--绘图 figure(1) plot(1:1:i,x(1,:),'b',1:1:i,x(2,:),'k',1:1:i,x(3,:),'g',1:1:i,x(4,:), 'r',1:1:i,x(5,:),'c') legend('x1','x2','x3','x4','x5') grid on title('Jacobi迭代法——x值随迭代次数变化曲线') figure(2) plot(1:1:i-1,er(1,:),'b',1:1:i-1,er(2,:),'k',1:1:i-1,er(3,:),'g',1:1: i-1,er(4,:),'r',1:1:i-1,er(5,:),'c') legend('△x1','△x2','△x3','△x4','△x5') grid on title('Jacobi迭代法——△x值随迭代次数变化曲线') %% fprintf('\n-------------Gauss-Seidel迭代法---------------------\n'); U=-(A1-D); L=-(A2-D); DL_1=inv(D-L); M1=DL_1*U; b2=DL_1*b; x1(:,1)=M1*x0+b2; for j=2:1:100 x1(:,j)=M1*x1(:,j-1)+b2; er1(:,j-1)=x1(:,j)-x1(:,j-1); if abs(er1(:,j-1))<=e fprintf('在迭代 %d 次之后，满足精度要求，x向量的值如下：\n',j); fprintf('x1=%.5f, x2=%.5f, x3=%.5f, x4=%.5f, x5=%.5f\n',x1(1,j),x1(2,j),x1(3,j),x1(4,j),x1(5,j)); break end end %--绘图 figure(3) plot(1:1:j,x1(1,:),'b',1:1:j,x1(2,:),'k',1:1:j,x1(3,:),'g',1:1:j,x1(4 ,:),'r',1:1:j,x1(5,:),'c') legend('x1','x2','x3','x4','x5')
常州大学数值分析第三章
第三章作业1.设节点x 0=0,x 1=π/8,x2=π/4,x3=3π/8,x4=π/2,试适当选取上述节点，用拉格朗日插值法分别构造cosx 在区间[0，π/2]上的一次、二次、四次差值多项式P 1（x ）,P 2(x)和P 4（x)，并分别计算P 1（π/3）,P 2(π/3)和P 4（π/3). 解： x0 x1 x2 x3 x4 xπ/8 π/4 3π/8 π/2 y=cosx 10.9238790.7071060.382683（1）选择x0=0,x4=π/2的节点y0=cosx0=1,y4=cosx4=0,可得)()()()()(010110101x x x x y x x x x y x P --+--=，即 333333.0)3/(1636620.0)(11≈+-≈πP x x P（2）选择x0=0,x2=π/4,x4=π/2的节点y0=cosx0=1,y2=cosx2=0.707106,y4=cosx4=0,可得))(())(())(())(())(())(()(1202102210120120102101x x x x x x x x y x x x x x x x x y x x x x x x x x y x P ----+----+----=,即145968.1)3/(1511124.5482067.1)(222≈++-≈πP x x x P（3）选择x0=0,,x1=π/8,x2=π/4,x3=3π/8，x4=π/2的节点y0=cosx0=1,y1=cosx1=0.923879,y2=cosx2=0.707106,y3=cosx3=0.382683,y4=cosx4=0可得)()(4,044∏∑≠==--=ij j ji j i i x x x x y x P ， 得P3(x)=1+0.0031x-0.51542x +0.02423x +0.02844x4(3) 0.5001P π=/7.解：选取0123=0=1=2=3x x x x ，，，为节点 >> T0=[0.0 0.5];x=[1 2 3]';y=[1.25 2.75 3.5]';x0=2.8;T=aitken(x,y,x0,T0) T =0.0000 0.5000 0.0000 0.0000 0.0000 1.0000 1.2500 2.6000 0.0000 0.0000 2.0000 2.7500 3.6500 4.4900 0.0000 3.0000 3.5000 3.3000 3.2300 3.4820 161）拉格朗日差值．选取函数],[),sin()cos(ππ-∈+=x x x yx0=-pi:0.5*pi:pi; y0=cos(x0);x=-pi:0.05*pi:pi;if length(x0)~=length(y0)error('The length of x0 must be equal to it of y0'); endw=length(x0); n=w-1;L=zeros(w,w); for k=1:n+1 V=1;for j=1:n+1 if k~=j ifabs(x0(k)-x0(j))<epserror('Divided by Zero,there are two nodes are the same'); endV=conv(V,poly(x0(j)))/(x0(k)-x0(j));end end L(k,:)=V;endC=y0*L;Y=polyval(C,x); y=cos(x)+sin(x); r=y-Y;plot(x,Y,'r--',x,y,'b-',x,r,'k-.','LineWidth',2);legend('Lagrange polynomial','Theoriginal f(x) ','Error',0)2）牛顿差值．选取函数]5,5[,112-∈+=x x y >> x0=-5:1:5;>>y0=1./(1+x0.*x0); >>x=-5:0.1:5;>>if length(x0)~=length(y0)>>error('The length of x0 must be equal to it of y0'); end>>n=length(x0); >>D=zeros(n,n); >>D(:,1)=y0'; >>for j=2:n for k=j:nif abs(x0(k)-x0(k-j+1))<epserror('Divided by Zero,there are two nodes are the same'); end>>D(k,j)=(D(k,j-1)-D(k-1,j-1))/(x0(k )-x0(k-j+1)); end endC=D(n,n);for k=(n-1):-1:1C=conv(C,poly(x0(k))); m=length(C);C(m)=C(m)+D(k,k); endY=polyval(C,x); y=1./(1+x.*x+20*x); r=y-Y;>plot(x,Y ,'r--',x,y,'b-',x,r,'k-.','LineWidth',2); >legend('Newton polynomial','The original f(x)','Error',0)。
(完整版)数值分析第一次作业
问题1：20.给定数据如下表：试求三次样条插值S(x)，并满足条件 (1)S`(0.25)=1.0000，S`(0.53)=0.6868； （2）S ’’(0.25)=S ’’(0.53)=0。
分析：本问题是已知五个点，由这五个点求一三次样条插值函数。
边界条件有两种，（1）是已知一阶倒数，（2）是已知自然边界条件。
对于第一种边界(已知边界的一阶倒数值)，可写出下面的矩阵方程。
⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡432104321034322110d M M M M M 200020000020022d d d d λμμλμλμλ其中μj =j1-j 1-j h h h +，λi=j1-j j h h h +，dj=6f[x j-1,x j ,x j+1]， μn =1，λ0=1对于第一种边界条件d 0=0h 6（f[x 0,x 1]-f 0`），d n =1-n h 6（f`n-f `[x n-1,x n ]） 解：由matlab 计算得：由此得矩阵形式的线性方程组为：⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡ 2.1150-2.4286-3.2667-4.3143-5.5200-M M M M M 25714.00001204286.000004000.026000.0006429.023571.0001243210解得 M 0=-2.0286;M 1=-1.4627;M 2= -1.0333; M 3= -0.8058; M 4=-0.6546S(x)=⎪⎪⎩⎪⎪⎨⎧∈-+-+-∈-+-+-∈-+-+-∈-+-+-]53.0,45.0[x 5.40x 9.1087x 35.03956.8.450-x 1.3637-x .5301.67881- ]45.0,39.0[x 9.30x 11.188x 54.010.418793.0-x 2.2384-x .450(2.87040-]39.0,30.0[x 03.0x 6.9544x 9.30 6.107503.0-x 1.9136-x .3902.708779-]30.0,25.0[x 5.20x 10.9662x 0.3010.01695.20-x 4.8758-x .3006.76209-33333333），（）（）（）（），（）（）（）），（）（）（）（），（）（）（）（Matlab 程序代码如下：function tgsanci(n,s,t) %n代表元素数，s,t代表端点的一阶导。
数值分析习题(含标准答案)
]第一章 绪论姓名 学号 班级习题主要考察点：有效数字的计算、计算方法的比较选择、误差和误差限的计算。
1若误差限为5105.0-⨯,那么近似数有几位有效数字（有效数字的计算） 解：2*103400.0-⨯=x ，325*10211021---⨯=⨯≤-x x 故具有3位有效数字。
2 14159.3=π具有4位有效数字的近似值是多少（有效数字的计算） 解：10314159.0⨯= π，欲使其近似值*π具有4位有效数字，必需！41*1021-⨯≤-ππ，3*310211021--⨯+≤≤⨯-πππ，即14209.314109.3*≤≤π即取（ , ）之间的任意数，都具有4位有效数字。
3已知2031.1=a ，978.0=b 是经过四舍五入后得到的近似值，问b a +，b a ⨯有几位有效数字（有效数字的计算）解：3*1021-⨯≤-aa ，2*1021-⨯≤-b b ，而1811.2=+b a ，1766.1=⨯b a 2123****102110211021)()(---⨯≤⨯+⨯≤-+-≤+-+b b a a b a b a故b a +至少具有2位有效数字。
2123*****10210065.01022031.1102978.0)()(---⨯≤=⨯+⨯≤-+-≤-b b a a a b b a ab 故b a ⨯至少具有2位有效数字。
4设0>x ，x 的相对误差为δ，求x ln 的误差和相对误差（误差的计算）~解：已知δ=-**xx x ，则误差为 δ=-=-***ln ln xx x x x则相对误差为******ln ln 1ln ln ln xxx x xxx x δ=-=-5测得某圆柱体高度h 的值为cm h 20*=，底面半径r 的值为cm r 5*=，已知cm h h 2.0||*≤-，cm r r 1.0||*≤-，求圆柱体体积h r v2π=的绝对误差限与相对误差限。
（误差限的计算）解：*2******2),(),(h h r r r h r r h v r h v -+-≤-ππ绝对误差限为πππ252.051.02052)5,20(),(2=⨯⋅+⨯⋅⋅⋅≤-v r h v相对误差限为%420120525)5,20()5,20(),(2==⋅⋅≤-ππv v r h v 6设x 的相对误差为%a ,求nx y =的相对误差。
常州大学数值分析作业(共六章)
第二章：20.（1）用 Jacobi 迭代法解方程组 AX=b. function [x,iternum,flag]=jacobi(A,b,x0,delta,max1) %检验输入参数，初始化 if nargin<2,error('more augments are needed');end if nargin<3,x0=zeros(size(b));end if nargin<4,delta=1e-13;end if nargin<5,max1=100;end if nargin>5,error('incorrect number of input');end n=length(b);x=0*b;flag=0;iternum=0; %用Jacobi迭代法解方程组 for k=1:max1 iternum=iternum+1; for i=1:n if abs(A(i,i))<eps error('A(i,i) equal to zero,divided by zero'); end x(i)=(b(i)-A(i,[1:i-1,i+1:n])*x0([1:i-1,i+1:n]))/A(i,i); end err=norm(x-x0); relerr=err/(norm(x)+eps); x0=x; if (err<delta)||(relerr<delta) flag=1; break; end end if flag==1 disp('The Jacobi method converges.') x=x; else disp(['The Jacobi method does not converge with '... ,num2str(max1),' iterations']) end return A=[1 2 -2;1 1 1;2 2 1]; b=[1;1;1]; [x,iternum,flag]=jacobi(A,b) 输出结果为： The Jacobi method converges. ans = -3 3 1 iternum =4 flag = 1
数值分析第四版习题及答案
第四版数值分析习题第一章绪论设x>O,x 的相对误差为S ，求In x 的误差. 设x 的相对误差为2%,求x n 的相对误差. 下列各数都是经过四舍五入得到的近似数，即误差限不超过最后一位的半个单位 ，试指出它们是几位有效数字： x = 1.1021, x^ = 0.031, x^ = 385.6, x^ = 56.430, x^ = 7 1.0.利用公式(3.3)求下列各近似值的误差限：(i)x *+x ；+x 4,(ii)x *x ；x ；,(iii )x ；/x ；,其中 x ；,x ；,x 3,x ；均为第 3题所给的数.计算球体积要使相对误差限为 1%，问度量半径R 时允许的相对误差限是多少 ？设\)=28,按递推公式AY n =Y n d- _ .783100( n=1,2,…)计算到Y 00.若取7783衣27.982(五位有效数字)，试问计算^00将有多大误差？ 求方程X 2 -56X • 1 =0的两个根，使它至少具有四位有效数字 (■ 783沁27.982).\ ------ d x 当N 充分大时，怎样求N 1 x? 正方形的边长大约为 100 cm ,应怎样测量才能使其面积误差不超过 s *2设 2 假定g 是准确的，而对t 的测量有土 0.1秒的误差，证明当t 增加时s 的绝对 误差增加，而相对误差却减小. 序列｛yn｝满足递推关系y n _ 10y n _ 1(n=1,2,…)，若y0 _ X 2 1.41 (三位有效数字),计算到y 10时误差有多大？这个计算过程稳定吗？计算f = c- 2 一1)6,取' 2 : 1.4,利用下列等式计算，哪一个得到的结果最好？f (x) =1 n (x X -1)，求 f(30)的值.若开平方用六位函数表，问求对数时误差有多大改用另一等价公式ln(x_ Jx 2 T) = -ln(x +Jx 2 +1)计算,求对数时误差有多大？1. 2. 3. 4.5. 6.7.8.9.10.11.12.13.21 cm1 (、2 1)61 (32 . 2)3,99 -70、2.?若根据(2.2)定义的范德蒙行列式，令证明V n (x)是n 次多项式,它的根是X 0^L ，X nJ ,且当x= 1 , -1 , 2时，f(x)= 0 , -3,4 ,求f(x)的二次插值多项式.给出cos x,0 ° < x 90。
数值分析第六章 课后习题 常州大学
2
数值分析作业三 1、试给出下述方程的有根区间或初试近似根： 解： （1） 3 x 3 2 x 2 0 令 f x 3 x 3 2 x 2 ，
f ' x 9 x 2 2 令f ' x 0, 得：x 则，f x 在 2 或x 3 2 ， 或 ， 3 2 3 2 单调递增。 3
2 0 ，所以f x 有且仅有一个零点。 且f 3 且有f 2 0 ，f 1 0 所以，原方程的根在 - 2， - 1区间内。
2、利用二分法求上述方程的根，要求误差不超过 10-2。 解：function [c,err,yc,k]=bisect(f,a,b,epsilon) yb=f(b);ya=f(a);max1=1+round((log(b-a)-log( epsilon))/log(2)); flag=1;k=0; while flag==1 end end c=(a+b)/2; err=abs(b-a); yc=f(c); k=max1; return k=1:max1; c=(a+b)/2; yc=f(c); if yc==0 a=c; b=c; elseif yb*yc>0 b=c; yb=yc; else a=-2;b=-1; n=20;epsilon=1e-2; f=@(x)(3x^3-2*x+2); [c,err,yc,k]=bisect(f,a,b,epsilon) 输出结果： c= err = yc = k= -1.1211 0.0078 0.0150 8 end if b-a<epsilon break a=c; ya=yc;
x
3
2x
常州大学数值分析第一章习题解答
1.1解：m=3;f=@(x)digit(digit(x^4,m)- digit(x^3,m)+ digit(3*x^2,m)+ digit(x-2,m),m);g=@(x)digit(digit(digit( digit(digit(digit( (x-1)*x,m)+3,m)*x,m)+1,m)*x,m)-2,m);f(3.33)g(3.33)有ans = 121ans =121实际上，当m=2时，就可以看出这两种算法在计算的精确度上的区别：m=2;f=@(x)digit(digit(x^4,m)- digit(x^3,m)+ digit(3*x^2,m)+ digit(x-2,m),m);g=@(x)digit(digit(digit( digit(digit(digit( (x-1)*x,m)+3,m)*x,m)+1,m)*x,m)-2,m);f(3.33)g(3.33)有ans = 120ans =130，可以看出，两者在计算精度上的不同区别，数学上恒等，在数值上不一定恒等。
1.2解：（1）精确到小数点后第三位，故有4位有效数字（2）精确到小数点后第三位，故有2位有效数字（3）精确到小数点后第三位，故有0位有效数字1.3 解;记圆的面积为S，由题意有|e(S)|≤1%。
由S=πr2知：dS=2πrdr所以dS/S=(2πrdr)/(πr2)=2(dr/r)∴|e(r)|≈1／2|e(S)|≤0.5×1%=0.5%1.4 解：由题有：|e(x)|≤1/2×10^-2 ; |e(y)|≤1/2×10^-2; |e(z)||≤1/2×10^-2∴|e（S）|≈|xe(x)+ye(y)|+ |ze(z)|^2≈x|e(x)|+y|e(y)|+z^2|z(z)|^2≤4.21×0.005+1.79×1.005+2.11×2.11×0.005^2=0.03≤1/2×10^-1又S=4.21*1.79+2.11^2=11.988∴S至少具有3位有效数字。
数值分析详细答案(全)
第二章 插值法习题参考答案2.)12)(12()1)(1(4)21)(11()2)(1()3()21)(11()2)(1(0)(2+-+-⋅+------⋅-+-+-+⋅=x x x x x x x L3723652-+=x x . 3. 线性插值：取510826.0,693147.0,6.0,5.01010-=-===y y x x ，则620219.0)54.0()54.0(54.0ln 0010101-=-⋅--+=≈x x x y y y L ；二次插值：取510826.0,693147.0,916291.0,6.0,5.0,4.0210210-=-=-====y y y x x x ，则)54.0(54.0ln 2L ≈))(()54.0)(54.0())(()54.0)(54.0())(()54.0)(54.0(120210221012012010210x x x x x x y x x x x x x y x x x x x x y ----⋅+----⋅+----⋅=＝－0.616707 .6. i) 对),,1,0(,)(n k x x f k==在n x x x ,,,10 处进行n 次拉格朗日插值，则有)()(x R x P x n n k +=)())(()!1(1)(0)1(0n n ni k j j x x x x f n x x l --++=+=∑ ξ由于0)()1(=+ξn f，故有kni k j jxx x l≡∑=0)(.ii) 构造函数,)()(kt x x g -=在n x x x ,,,10 处进行n 次拉格朗日插值，有∑=-=ni j k j n x l t x x L 0)()()(.插值余项为 ∏=+-+=--nj j n n kx x n g x L t x 0)1()()!1()()()(ξ， 由于).,,2,1(,0)()1(n k g n ==+ξ故有 .)()()()(0∑=-==-ni j k j n kx l t x x L t x令,x t =即得 ∑==-ni j k jx l t x)()(.8. 截断误差].4,4[),)()((61)(2102-∈---=ξξx x x x x x e x R其中 ,,1210h x x h x x +=-= 则hx x 331+=时取得最大值321044392|))()((|max h x x x x x x x ⋅=---≤≤- .由题意， ,10)392(61|)(|6342-=⋅⋅≤h e x R所以，.006.0≤h16. ;1!7!7!7)(]2,,2,2[)7(71===ξf f .0!7)(]2,,2,2[)8(810==ξf f19. 采用牛顿插值，作均差表：i x)(i x f一阶均差 二阶均差0 1 20 1 11 0-1/2],,[))((],[)()()(210101000x x x f x x x x x x f x x x p x p --+-+=))()()((210x x x x x x Bx A ---++)2)(1()()2/1)(1(0--++--++=x x x Bx A x x x又由 ,1)1(,0)0(='='p p 得,41,43=-=B A 所以 .)3(4)(22-=x x x p第三章 函数逼近与计算习题参考答案4.设所求为()g x c =，(,)max(,),max (),min ()a x ba x bf g M c m c M f x m f x ≤≤≤≤∆=--==，由47页定理4可知()g x 在[],a b 上至少有两个正负交错的偏差点，恰好分别为()f x 的最大值和最小值处，故由1(),()2M c m c c M m -=--=+可以解得1()()2g x M m =+即为所求。
常州大学数值分析作业—第二章
20.分别用Jacobi 迭代法、Gauss-seidel 迭代法解方程组⎪⎩⎪⎨⎧=++=++=-+.122,1,122321321321x x x x x x x x x解：Jacobi 迭代法收敛⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡-=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡133321x x x ，Gauss-seidel 迭代法不收敛。
27.编写LU 分解法、改进平方根法、追赶法的Matlab 程序，并进行相关数值实验。
3.将矩阵⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡-=1100110211100201A 进行Doolittle 和Crout 分解 解：Doolittle 分解：结果如下，程序见后面。
⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡-⋅⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡-==2.100015001110020112.000010200100001LU A Crout 分解：结果如下，程序见后面。
⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡-⋅⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡-==1002.0100111002012.1100050200100001LU A 7.用改进平方根法解方程组⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡⋅⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡----000142002511013101144321x x x x 解：结果如下，程序见后面。
⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡--=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡0256.00513.00769.02821.04321x x x x 8（2）.用追赶法求解方程组⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡⋅⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡--------00001210001210001210001210001254321x x x x x解：结果如下，程序见后面。
⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡---=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡1429.02413.02857.03571.04321x x x x%检验输入参数，初始化if nargin<2,error('more augments are needed');endif nargin<3,x0=zeros(size(b));endif nargin<4,delta=1e-13;endif nargin<5,max1=100;endif nargin>5,error('incorrect number of input');endn=length(b);x=0*b;flag=0;iternum=0;%用Jacobi迭代法解方程组for k=1:max1iternum=iternum+1;for i=1:nif abs(A(i,i))<epserror('A(i,i) equal to zero,divided by zero');endx(i)=(b(i)-A(i,[1:i-1,i+1:n])*x0([1:i-1,i+1:n]))/A(i,i); x(i)=(b(i)-A(i,[1:i-1,i+1:n])*x0([1:i-1,i+1:n]))/A(i,i);enderr=norm(x-x0);relerr=err/(norm(x)+eps);x0=x;if (err<delta)||(relerr<delta)flag=1;break;endendif flag==1disp('The Jacobi method converges.')x=x;elsedisp(['The Jacobi method does not converge with '...,num2str(max1),' iterations'])endA=[1 2 -2;1 1 1;2 2 1] b=[1;1;1]A =1 2 -21 1 12 2 1b =111 jacobi(A,b)The Jacobi method converges. ans =-331%检验输入参数，初始化if nargin<2,error('more augments are needed');end if nargin<3,x0=zeros(size(b));end if nargin<4,delta=1e-13;end if nargin<5,max1=100;endif nargin>5,error('incorrect number of input');end n=length(b);flag=0;iternum=0;%用Gauss-Seidel 迭代法解方程组 for k=1:max1iternum=iternum+1; x=x0;for i=1:nif abs(A(i,i))<epserror('A(i,i) equal to zero,divided by zero'); endx0(i)=(b(i)-A(i,[1:i-1,i+1:n])*x0([1:i-1,i+1:n]))/A(i,i); enderr=norm(x-x0);relerr=err/(norm(x0)+eps);if (err<delta)||(relerr<delta) flag=1; break ; end endif flag==1disp('The Gauss-seidel method converges.') x=x0; elsedisp(['The Gauss-seidel method does not converge with '... ,num2str(max1),' iterations']) end returngseid(A,b) The Gauss-seidel method does not converge with 100 iterations ans = 1.0e+31 * -9.1905 9.2222 -0.0634 A=[1 2 -2;1 1 1;2 2 1] b=[1;1;1] A = 1 2 -2 1 1 1 2 2 1 b = 1 1 1%检验输入参数，初始化b=size(A);n=b(1);if b(1)~=b(2)error('Matrix A should be a Square matrix.'); endif n~=rank(A)error('Matrix A should be FULL RANK.');endL=zeros(n,n);U=zeros(n,n);for i=1:nU(i,i)=1;end%将矩阵A进行Crout分解for k=1:nfor i=k:ntemp_sum=0;for t=1:k-1temp_sum=temp_sum+L(i,t)*U(t,k);endL(i,k)=A(i,k)-temp_sum;endfor j=k+1:ntemp_sum=0;for t=1:k-1temp_sum=temp_sum+L(k,t)*U(t,j);endU(k,j)=(A(k,j)-temp_sum)/L(k,k);endendreturnA=[1 0 2 0;0 1 1 1;2 0 -1 1;0 0 1 1]A =1 02 00 1 1 12 0 -1 10 0 1 1 [L,U]=Crout(A)L =1.0000 0 0 00 1.0000 0 02.0000 0 -5.0000 00 0 1.0000 1.2000 U =1.0000 02.0000 00 1.0000 1.0000 1.00000 0 1.0000 -0.20000 0 0 1.0000%检验输入参数，初始化b=size(A);n=b(1);if b(1)~=b(2)error('Matrix A should be a Square matrix.'); endif n~=rank(A)error('Matrix A should be FULL RANK.');endL=zeros(n,n);U=zeros(n,n);for i=1:nU(i,i)=1;end%将矩阵A进行Doolittle分解for k=1:nfor i=k:ntemp_sum=0;for t=1:k-1temp_sum=temp_sum+L(i,t)*U(t,k);endL(i,k)=A(i,k)-temp_sum;endfor j=k+1:ntemp_sum=0;for t=1:k-1temp_sum=temp_sum+L(k,t)*U(t,j);endU(k,j)=(A(k,j)-temp_sum)/L(k,k);endendreturnA=[1 0 2 0;0 1 1 1;2 0 -1 1;0 0 1 1]A =1 02 00 1 1 12 0 -1 10 0 1 1 [L,U]=Doolittle(A)L =1.0000 0 0 00 1.0000 0 02.0000 0 1.0000 00 0 -0.2000 1.0000 U =1.0000 02.0000 00 1.0000 1.0000 1.00000 0 -5.0000 1.00000 0 0 1.2000function [x]=improvecholesky(A,b,n)%检验输入参数，初始化L=zeros(n,n);D=diag(n,0);S=L*D;for i=1:nL(i,i)=1;endfor i=1:nfor j=1:nif (eig(A)<=0)|(A(i,j)~=A(j,i))disp('Matrix A should be a Positive-definite matrix.');break;endendendD(1,1)=A(1,1);%用改进平方根法解方程组for i=2:nfor j=1:i-1S(i,j)=A(i,j)-sum(S(i,1:j-1)*L(j,1:j-1)');L(i,1:i-1)=S(i,1:i-1)/D(1:i-1,1:i-1);endD(i,i)=A(i,i)-sum(S(i,1:i-1)*L(i,1:i-1)');endy=zeros(n,1);x=zeros(n,1);for i=1:ny(i)=(b(i)-sum(L(i,1:i-1)*D(1:i-1,1:i-1)*y(1:i-1)))/D(i,i); endfor i=n:-1:1x(i)=y(i)-sum(L(i+1:n,i)'*x(i+1:n));endreturnA=[4 1 -1 0;1 3 -1 0;-1 -1 5 2;0 0 2 4] b=[1;0;0;0]A =4 1 -1 01 3 -1 0-1 -1 5 20 0 2 4b =10 n=4[x]=improvecholesky(A,b,n)n =4x =0.2821-0.07690.0513-0.0256function [L,U,x]=pursue(a,b,c,f)n=length(b);u(1)=b(1);for i=2:nif(u(i-1)~=0)l(i-1)=a(i-1)/u(i-1);u(i)=b(i)-l(i-1)*c(i-1);elsebreak;endendL=eye(n)+diag(l,-1);U=diag(u)+diag(c,1);x=zeros(n,1);y=x;y(1)=f(1);for i=2:ny(i)=f(i)-l(i-1)*y(i-1);endif(u(n)~=0)x(n)=y(n)/u(n);endfor i=n-1:-1:1x(i)=(y(i)-c(i)*x(i+1))/u(i);enda=[1 -1 -1 -1];b=[2 2 2 2 2];c=[-1 -1 -1 -1];f=[1 0 0 0 0 ];[L,U,x]=pursue(a,b,c,f)L =1.0000 0 0 0 00.5000 1.0000 0 0 00 -0.4000 1.0000 0 00 0 -0.6250 1.0000 00 0 0 -0.7273 1.0000 U =2.0000 -1.0000 0 0 00 2.5000 -1.0000 0 00 0 1.6000 -1.0000 00 0 0 1.3750 -1.00000 0 0 0 1.2727 x =0.3571-0.2857-0.2143-0.1429。
常州大学2012-2013级研究生数值分析试卷A解答及评分标准
一．（1）已知函数24()73f x x x =++，用秦九昭方法计算(2)f ；（2）秦九昭方法计算任一n 次多项式在任一点函数值至多需要多少次乘法？ （3）至少写出四种减少误差危害的常用手段。
解：（1）2422()73(31)7f x x x x x =++=++22(2)(321)2759f =⨯++=………… 5 分（2） 秦九昭方法计算任一n 次多项式在任一点函数值至多需要n 次乘法。
………… 5 分（3） A ）防止大数“吃”小数； B ）避免除数绝对值远远小于被除数绝对值的除法；C ）避免相近数相减；D ）避免使用不稳定的算法；E ）注意简化计算步骤，减少运算次数；………… 5 分二．给定方程组123311413132156x x x ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥-=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥-⎣⎦⎣⎦⎣⎦ （1）以分量形式写出解此线性方程组的Jacobi 迭代格式和Gauss -Seidel 迭代格式； （2）求1A 和A∞；（3）判断Gauss -Seidel 迭代格式的敛散性。
解：（1）Jacobi 迭代(1)()()123(1)()()213(1)()()312(4)/3(3)/3(62)/5k k k k k k k k k x x x x x x x x x +++=--=+-=-+， 0,1,2,k = Gauss-Seidel 迭代(1)()()123(1)(1)()213(1)(1)(1)312(4)/3(3)/3(62)/5k k k k k k k k k x x x x x x x x x ++++++=--=+-=-+， 0,1,2,k =………… 5 分（2）17A =，8A∞=；………… 5 分（3）因为方程组系数矩阵严格对角占优，所以Gauss -Seidel 迭代格式收敛。
………… 5 分三． 已知方程2()30x f x e x =-=，（1）证明该方程在区间[0.6,1.2]上存在唯一实根； （2）叙述牛顿法求方程()0f x =根的方法思想；（3）以初值01x =，用牛顿法求上述方程的近似解，要求误差不超过210- 。
常州大学数值分析第五章习题解答
化专012 闵建中 121029142012022302021342.(1)2()02()314==33====A A A h h A A h A A h A A h A h ⎧⎪++=⎪--=⎨⎪⎪+=⎩=解：将f(x)=1,x,x 分别代入求积公式，并令左右相等，得解得：，进一步，取f(x)x 代入上述求积公式，有左边右边0，但当f(x)x 时，上述积分公式左右不相等，因此该求积公式具有三次代数精度。
11221122.21=1=*(13=2311==*(13)=3311()((1)2(0)3())3324==,=139A A x x f x dx f f f -++-+≈-++⎰（） 解：当f(x)时，令左右相等有：2）所以有当f(x)x 时，令左右相等有：0所以有所以求积公式为：当f(x)x 时，有左边右边，左右不相等，所以其代数精度为4.解：(1) : T=0.5*(3+1)=2 S=2 C=2 err1 = 0 err2=0(2) : T=1*(1-3)=-2 S=-10/3 C=-10/3 err 1= 4/3 err2=0(3) : T=0 S=0 C=0 err1=0 err2=0(4) : T=e/2 S=1/6*(e+2exp(0.5)) C= 1 err1=-1.7183 err2= 0.00267.解：复合梯形公式T 2n 的matlab 实现：function I = trapezoid(fun,a,b,n)n = 2*n;h = (b-a)/(n-1);x = a : h :b;f = feval(fun , x);I = h * (0.5*f(1)+sum(f(2:n-1))+ 0.5*f(n));分别用复合梯形公式和复合辛普森公式计算积分，并与精确值比较，这里积分的精确值由matlab int函数求出，不手工计算。
function trapezoid_and_sinpsom clc;format longsyms xIexact = int(x*exp(x^2),x,0,1);a=0;b=1;for n=2:1:4t = trapezoid(@f,a,b,n)s = simpson(@f,a,b,n)err1 = vpa(Iexact - t,5)err2 = vpa(Iexact - s,5)endfunction y=f(x)y = x*exp(x^2);return从而得出第一小题的结果：n = 2t = 1.6697s = 0.8695err1 = -0.14144err2 = -0.0018562n = 3t = 0.923798756293777s = 0.859533825596209err1 = -0.064658err2 = -0.00039291n = 4t = 0.895892057505771s = 0.859268455239111err1 = -0.036751err2 = -0.00012754将f(x)简单修改就可以得到第2，3小题的结果8.解：function test_romberg()a=1;b=3;n=2;maxit=10;tol=0.5e-4; [R]=romberg(@f,a,b,n,tol,maxit);R = R(1,4)Rexact = log(3)err = R-Rexactreturnfunction y=f(x)y=1/x;return10.解：高斯-勒让德法求解积分的matlab实现：function I = guasslegendre(fun,a,b,n,tol)% guasslegendre: guass-legendre method to calculation integral value%Input : fun is the integrand input as a sring 'fun'% a and b are upper and lower limits of integration % n is the initial number of subintervals% tol is the tolerance%Output : I is Integral result% calculate quadrature nodesyms xp=sym2poly(diff((x^2-1)^(n+1),n+1))/(2^n*factorial(n));tk=roots(p);% calculate quadrature coefficientAk=zeros(n+1,1);for i=1:n+1xkt=tk;xkt(i)=[];pn=poly(xkt);fp=@(x)polyval(pn,x)/polyval(pn,tk(i));Ak(i)=quadl(fp,-1,1,tol);end% integral variable substitution, transformat [a, b] to [1, 1]xk=(b-a)/2*tk+(b+a)/2;% calculate the value of the integral function after variable substitutionfx=fun(xk)*(b-a)/2;%calculation integral valueI=sum(Ak.*fx);对第一小题，精确值由matlab int函数求出，不手工计算；fun = inline('x.^5+3.*x.^2+2');a = 0;b = 1;tol = 1e-5;for n =1:2;nI = guasslegendre(fun,a,b,n,tol)Iexact = int(x^5+3*x^2+2,x,0,1)err = Iexact-Iend计算结果如下：n = 1I = 3.7778Iexact = 19/6 = 3.166666666666667err = 1/72n = 2I = 3.166666666666666Iexact = 19/6 = 3.166666666666667err = 0由此可见高斯-勒让德方法可以在较少的求积节点数下得到较为精确的数值积分解！第2,3小题同理可得。
常州大学数值分析第一章习题解答
1.1解：m=3;f=@(x)digit(digit(x^4,m)- digit(x^3,m)+ digit(3*x^2,m)+ digit(x-2,m),m);g=@(x)digit(digit(digit( digit(digit(digit( (x-1)*x,m)+3,m)*x,m)+1,m)*x,m)-2,m);f(3.33)g(3.33)有ans = 121ans =121实际上，当m=2时，就可以看出这两种算法在计算的精确度上的区别：m=2;f=@(x)digit(digit(x^4,m)- digit(x^3,m)+ digit(3*x^2,m)+ digit(x-2,m),m);g=@(x)digit(digit(digit( digit(digit(digit( (x-1)*x,m)+3,m)*x,m)+1,m)*x,m)-2,m);f(3.33)g(3.33)有ans = 120ans =130，可以看出，两者在计算精度上的不同区别，数学上恒等，在数值上不一定恒等。
1.2解：（1）精确到小数点后第三位，故有4位有效数字（2）精确到小数点后第三位，故有2位有效数字（3）精确到小数点后第三位，故有0位有效数字1.3 解;记圆的面积为S，由题意有|e(S)|≤1%。
由S=πr2知：dS=2πrdr所以dS/S=(2πrdr)/(πr2)=2(dr/r)∴|e(r)|≈1／2|e(S)|≤0.5×1%=0.5%1.4 解：由题有：|e(x)|≤1/2×10^-2 ; |e(y)|≤1/2×10^-2; |e(z)||≤1/2×10^-2∴|e（S）|≈|xe(x)+ye(y)|+ |ze(z)|^2≈x|e(x)|+y|e(y)|+z^2|z(z)|^2≤4.21×0.005+1.79×1.005+2.11×2.11×0.005^2=0.03≤1/2×10^-1又S=4.21*1.79+2.11^2=11.988∴S至少具有3位有效数字。
常州大学数值分析课后习题答案第二章第三章第四章节资料
数值分析作业第二章1、用Gauss消元法求解下列方程组：2x1-x2+3x3=1,(1) 4x1+2x2+5x3=4,x1+2x2=7;(2) 解：A=[2 -1 3 1;4 2 5 4;1 2 0 7]n=size(A,1);x=zeros(n,1);flag=1;% 消元过程for k=1:n-1for i=k+1:nif abs(A(k,k))>epsA(i,k+1:n+1)=A(i,k+1:n+1)-A(k,k+1:n+1)*A(i,k)/A(k,k); elseflag=0;returnendendend% 回代过程if abs(A(n,n))>epsx(n)=A(n,n+1)/A(n,n);elseflag=0;returnendfor i=n-1:-1:1x(i)=(A(i,n+1)-A(i,i+1:n)*x(i+1:n))/A(i,i); endreturnxA = 2 -1 3 14 25 41 2 0 7x = 9-1-611x1-3x2-2x3=3,(2)-23x1+11x2+1x3=0,x1+2x2+2x3=-1;(2) 解：A=[11 -3 -2 3;-23 11 1 0;1 2 2 -1]n=size(A,1);x=zeros(n,1);flag=1;% 消元过程for k=1:n-1for i=k+1:nif abs(A(k,k))>epsA(i,k+1:n+1)=A(i,k+1:n+1)-A(k,k+1:n+1)*A(i,k)/A(k,k);elseflag=0;returnendendend% 回代过程if abs(A(n,n))>epsx(n)=A(n,n+1)/A(n,n);elseflag=0;returnendfor i=n-1:-1:1x(i)=(A(i,n+1)-A(i,i+1:n)*x(i+1:n))/A(i,i);endreturnxA = 11 -3 -2 3-23 11 1 01 2 2 -1x = 0.21240.5492-1.15544、用Cholesky分解法解方程组3 2 3 x1 52 2 0 x2 33 0 12 x3 7解：.A=[3 2 3;2 2 0;3 0 12];b=[5 3 7];lambda=eig(A);if lambda>eps&isequal(A,A')[n,n]=size(A);R=chol(A);%解R'y=by(1)=b(1)/R(1,1);if n>1for i=2:ny(i)=(b(i)-R(1:i-1,i)'*y(1:i-1)')/R(i,i);endend%解Rx=yx(n)=y(n)/R(n,n);if n>1for i=n-1:-1:1x(i)=(y(i)-R(i,i+1:n)*x(i+1:n)')/R(i,i);endendx=x';elsex=[];disp('该方法只适用于对称正定的系数矩阵!');endR= 1.7321 1.1547 1.73210 0.8165 -2.44950 0 1.7321y= 2.8868 -0.4082 0.5774x= 1.0000 0.5000 0.33335. 用列主元Doolittle分解法解方程组解：A=[3 4 5; -1 3 4; -2 3 -5;]; 3 4 5 X1 2 b=[2,-2 6]'; -1 3 4 X2 -2 [L,U,pv]=luex(A); -2 3 -5 X3 6y = L\b(pv);x = U\y结果如下：x = 11-114.已知，计算.解：A=[100 99;99 98];cond(A,inf)ans =3.9601e+04cond(A,2)ans =3.9206e+0427.编写LU分解法，改进平方根法，追赶法的Matlab程序，并进行相关数值试验。
数值分析课后习题及答案
数值分析课后习题及答案第一章绪论（12）第二章插值法（40-42）2、当时，，求的二次插值多项式。
[解]。
3、给出的数值表用线性插值及二次插值计算的近似值。
X 0.4 0.5 0.6 0.7 0.8 -0.916291 -0.693147 -0.510826 -0.357765 -0.223144 [解]若取，，则，，则，从而。
若取，，，则，，，则，从而补充题：1、令，，写出的一次插值多项式，并估计插值余项。
[解]由，可知，，余项为，故。
2、设，试利用拉格朗日插值余项定理写出以为插值节点的三次插值多项式。
[解]由插值余项定理，有，从而。
5、给定数据表：，1 2 4 6 7 4 1 0 1 1 求4次牛顿插值多项式，并写出插值余项。
[解]一阶差商二阶差商三阶差商四阶差商 1 42 1 -34 0 6 17 1 0 由差商表可得4次牛顿插值多项式为：，插值余项为。
第三章函数逼近与计算（80-82）26、用最小二乘法求一个形如的经验公式，使它与下列数据相拟合，并求均方误差。
19 25 31 38 44 19.0 32.3 49.0 73.3 97.8[解]由。
又，，，故法方程为，解得。
均方误差为。
27、观测物体的直线运动，得出以下数据：时间t（秒）0 0.9 1.9 3.0 3.9 5.0 距离s（米）0 10 30 5080 110 [解]设直线运动为二次多项式，则由。
，。
又，，，故法方程为，解得。
故直线运动为。
补充题：1、现测得通过某电阻R的电流I及其两端的电压U如下表：I ……U ……试用最小二乘原理确定电阻R的大小。
[解]电流、电阻与电压之间满足如下关系：。
应用最小二乘原理，求R使得达到最小。
对求导得到：。
令，得到电阻R为。
2、对于某个长度测量了n次，得到n个近似值，通常取平均值作为所求长度，请说明理由。
[解]令，求x使得达到最小。
对求导得到：，令，得到，这说明取平均值在最小二乘意义下误差达到最小。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
第一章：9.设2cos 1)(x xx f -=，给出计算函数值)012.0(f 的一个合适算法，并在字长m 给定的，十进制计算机上给出数值计算结果。
解：由 )2421(242)2421(1)cos(1224242x x x x x x x -=-=+--≈- 得 )2421(cos 1)(22x x x x f -≈-=10. 字长为5的十进制计算机上计算)015.0(f 和)015.0(g ，并与)015.0(f 的精确值1.0075376410479比较，说明差异存在理由，其中x e x f x 1)(-=，24621)(32x x x x g +++=。
clearf=@(x)1/2-x^2/24; f(0.012)ans =0.5000解：字长为5时的误差很大，这是因为设置的字长有限，就不可避免的使舍入误差不断积累。
把字长改为9时，误差已经大幅度减小。
这说明，加大字长可以显著减小误差。
11. 举例介绍数组矩阵常见运算。
解：举例如下clearf=@(x)digit(digit(exp(x)-1,5)/x,5);g=@(x)digit(digit(1,5)+digit(x/2,5)+digit... (digit(x^2,5)/6,5)+digit(digit(x^3,5)/24,5),5); exc=1.0075376410479; f(0.015) g(0.015)err1=f(0.015)-exc err2=g(0.015)-excans =1.0075 ans =1.0075 err1 =-3.7641e-05 err2 =-3.7641e-05A*Bans =30 30 30 30 70 70 70 70 110 110 110 110 150 150 150 150A.*B ans =1 2 3 4 10 12 14 16 27 30 33 36 52 56 60 64 A^2 ans =A.^2 ans =f=@(x)digit(digit(exp(x)-1,9)/x,9);g=@(x)digit(digit(1,9)+digit(x/2,9)+digit... (digit(x^2,9)/6,9)+digit(digit(x^3,9)/24,9),9); err1=f(0.015)-exc err2=g(0.015)-excerr1 =-1.0479e-09 err2 =-1.0479e-09clearA=[1:4;5:8;9:12;13:16]B=[1,1,1,1;2,2,2,2;3,3,3,3;4,4,4,4] A ’ A =1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 B =1 1 1 12 2 2 23 3 3 34 4 4 4 ans =%%编写m 文件使用digit 函数设置字长%% function y=digit(x,m) k=max(size(x)); y=x;for i=1:k if x(i)<0 sign=-1; elsesign=1; endx(i)=abs(x(i)); p=0;if x(i)<0.1&x(i)>eps while x(i)<0.1 x(i)=x(i)*10; p=p-1; end endif x(i)>=1while x(i)>=1 x(i)=x(i)/10; p=p+1; end endy(i)=round(x(i)*10^m)/10^m; y(i)=sign*y(i)*10^p; end return12.对任意给定的实数a 、b 、c 、试编写Matlab 程序,求方程02=++c bx ax 的根。
解：利用教材例11的方法： 当b>0时，a ac b b x 2421---=，bac b c x +--=4222。
当b<0时，bac b cx --=4221，a acb b x 2422-+-=。
13.利用1,753arctan 753<+-+-=x x x x x x 及()3/3arctan 6=π，给出一个计算π的方法，根据此方法编写程序，给出π的至少有10位有效数字的近似值。
解：根据题中所给公式，容易得到：A+B ans = 2 3 4 5 7 8 9 10 12 13 14 15 17 18 19 20 A-B ans =0 1 2 3 3 4 5 6 6 7 8 9 9 10 11 12cleara=input('输入a='); b=input('输入b='); c=input('输入c='); d=b^2-4*a*c if b>=0x1=(-b-d^0.5)/2*a x2=-2*c/(d^0.5+b) else b<0x1=-2*c/(d^0.5-b) x2=(-b+d^0.5)/2*a endclearx=3^(-0.5);y=atan(x) ; %精确值% s=0 ; %计算值% for k=1:2:100;s=s+(-1)^((k+1)/2)*(x^k)/k;err=y-s; m=abs(err) if m<=1e-11 输入a=2输入b=2 输入c=1 d = -4 x1 =-2.0000 - 2.0000i x2 =-0.5000 + 0.5000i 输入a=1 输入b=2 输入c=1 d = 0 x1 = -1 x2 = -1 输入a=1 输入b=5 输入c=2 d = 17 x1 =-4.5616 x2 =-0.4384()12)3/3(16)3/3arctan(61211--≈=-=+∑i i ni i π14.分别利用下式给出计算ln2的近似方法，编写相应程序并比较算法运行情况。
11,32)1()1ln(3211≤<-+++-=-=+∑∞=+x nx x x x n x x nn n n 11),1253(212211ln 1253112<<-+-++++=-=-+-∞=-∑x n x x x x n x x x n n n解：由运行结果可知，方法二的绝对误差比方法一的误差要小得多。
这是因为方法一给出的计算公式含有相近数相减项，损失了有效数字。
而方法二给出的计算公式避免了相近数相减，具有较好的精度。
第二章：20.（1）用Jacobi 迭代法解方程组AX=b.function [x,iternum,flag]=jacobi(A,b,x0,delta,max1) %检验输入参数，初始化if nargin<2,error('more augments are needed');end if nargin<3,x0=zeros(size(b));end if nargin<4,delta=1e-13;endif nargin<5,max1=100;endclear x=1/3; s=0;for k=1:2:100; s=s+((x)^k)/k; endy=2*s %计算值% g=log(2) %真实值% err=y-g %绝对误差% y =0.6931 g =0.6931 err =-2.2204e-16clearx=1; s=0;for k=1:100;s=s+((-1)^(k+1))*(x^k)/k; endy=s %计算值% g=log(2) %真实值% err=y-g %绝对误差% y =0.6882 g =0.6931 err =-0.005020.（1）用Gauss-Seidel 迭代法解方程组AX=b.function [x,iternum,flag]=gseid(A,b,x0,delta,max1) %检验输入参数，初始化if nargin<2,error('more augments are needed');end if nargin<3,x0=zeros(size(b));end if nargin<4,delta=1e-13;endif nargin<5,max1=100;endA=[1 2 -2;1 1 1;2 2 1]; b=[1;1;1]; [x,iternum,flag]=jacobi(A,b) 输出结果为： The Jacobi method converges. ans = -3 3 1 iternum =4 flag = 13.（1）将矩阵A 进行Crout 分解 A=[1 2 -2;1 1 1;2 2 1]; b=[1;1;1]; [x,iternum,flag]=gseid(A,b) 输出结果为： The Gauss-seidel method does not converge with 100 iterations ans = 1.0e+31 * -9.1905 9.2222 -0.0634 iternum = 4 flag = 1function [L,U]=Crout(A) %检验输入参数，初始化 b=size(A);n=b(1); if b(1)~=b(2)error('Matrix A should be a Square matrix.'); end3.（2）将矩阵A 进行Doolittle 分解。
A=[1 0 2 0;0 1 1 1;2 0 -1 1;0 0 1 1] A = 1 0 2 0 0 1 1 1 2 0 -1 1 0 0 1 1 [L,U]=Crout(A) L = 1.0000 0 0 0 0 1.0000 0 0 2.0000 0 -5.0000 0 0 0 1.0000 1.2000 U = 1.0000 0 2.0000 00 1.0000 1.0000 1.00000 0 1.0000 -0.20000 0 0 1.0000function [L,U]=Doolittle(A) %检验输入参数，初始化 b=size(A);n=b(1); if b(1)~=b(2)error('Matrix A should be a Square matrix.'); end7.用改进平方根法解方程组AX=bA=[1 0 2 0;0 1 1 1;2 0 -1 1;0 0 1 1]A =1 02 00 1 1 12 0 -1 10 0 1 1[L,U]=Doolittle(A) L = 1.0000 0 0 0 0 1.0000 0 0 2.0000 0 1.0000 0 0 0 -0.2000 1.0000 U = 1.0000 0 2.0000 00 1.0000 1.0000 1.00000 0 -5.0000 1.00000 0 0 1.2000function [x]=improvecholesky(A,b,n) %检验输入参数，初始化L=zeros(n,n);D=diag(n,0);S=L*D; for i=1:nL(i,i)=1;end8.（2）用追赶法求解方程组.function [L,U,x]=pursue(a,b,c,f) n=length(b);u(1)=b(1); for i=2:nif (u(i-1)~=0)l(i-1)=a(i-1)/u(i-1);u(i)=b(i)-l(i-1)*c(i-1); A=[4 1 -1 0;1 3 -1 0;-1 -1 5 2;0 0 2 4] b=[1;0;0;0] A = 4 1 -1 0 1 3 -1 0 -1 -1 5 2 0 0 2 4 b = 1 0 0 0n=4[x]=improvecholesky(A,b,n) n = 4 x =0.2821 -0.0769 0.0513 -0.0256第三章：1.设节点x 0=0，x 1=π/8，x 2=π/4，x 3=3π/8，x 4=π/2，试适当选取上述节点，用Lagrange 插值法分别构造cos x 在区间[0，π/2]上的一次、二次、四次差值多项式P 1（x ），P 2（x ）和P 4（x ），并分别计算P 1（π/3），P 2（π/3）和P 4（π/3）。