常州大学数值分析第三章
- 格式:doc
- 大小:144.39 KB
- 文档页数:2
第一章:9.设2cos 1)(xxx 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-05clearA=[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 4ans =1 5 9 132 6 10 143 7 11 154 8 12 16A*Bans = 30 30 30 30 70 70 70 70 110 110 110 110 150 150 150 150 A.*Bans =1 2 3 410 12 14 1627 30 33 3652 56 60 64A^2ans = 90 100 110 120 202 228 254 280 314 356 398 440 426 484 542 600 A.^2 ans =1 4 9 16 25 36 49 64 81 100 121 144 169 196 225 256%%编写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 returnf=@(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-0912.对任意给定的实数a 、b 、c 、试编写Matlab 程序,求方程02=++c bx ax 的根。
江苏工业学院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 =-+。
常州大学数值分析作业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。
数值分析第三章答案【篇一:常州大学数值分析作业第三章】答:matlab 程序function [a,y]=lagrange(x,y,x0) %检验输入参数if nargin 2 || nargin 3error(incorrect number of inputs); endif length(x)~=length(y)error(the length of x must be equal to it of y); endm=length(x);n=m-1;l=zeros(m,m); %计算基本插值多项式的系数for i=1:n+1 c=1;for j=1:n+1if i~=jif abs(x(i)-x(j))eps abs(x(i)-x(j))epserror(there are two two same nodes);endc=conv(c,poly(x(j)))/(x(i)-x(j));end endl(i,:)=c; end%计算lagrange插值多项式的系数 a=y*l;%计算f(x0)的近似值 if nargin==3y=polyval(a,x0);工程(专)学号:14102932enda=fliplr(a); return[a,y] = lagrange(x,y,x0); p1 = vpa(poly2sym(a),3) y[a,y] = lagrange(x,y,x0); p2=vpa(poly2sym(a),3) yp2 = x2 - 0.109x - 0.336 y =0.5174[a,y]=lagrange(x,y,x0); p4=vpa(poly2sym(a),3) yp4 =x4 + 0.00282x3 - 0.514x2 + 0.0232x + 0.0287 y =0.5001次多项式在2.8处的值。
答:matlab 程序 function[t,y0]=aitken(x,y,x0,t0) if nargin==3 t0=[]; 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; if n0==0 i0=2; elsei0=n0+1; endfor i=i0:nfor j=3:i+1t(i,j)=fun(t(j-2,1),t(i,1),t(j-2,j-1),t(i,j-1),x0); end endy0=t(n,n+1); returnfunction [y]=fun(x1,x2,y1,y2,x) y=y1+(y2-y1)*(x-x1)/(x2-x1); return%选取0、1、3、4四个节点,求三次插值多项式 x=[0,1,3,4];y=[0.5,1.25,3.5,2.75]; x0=2.8;[t,y0]=aitken(x,y,x0) t =0 0.5000 00 0 1.01.25002.6000 0 0 3.03.50003.29993.23000 4.02.75002.07502.28503.4190 y0 =3.41900000000000016、选取适当的函数y=f(x)和插值节点,编写matlab程序,分别利用lagrange插值方法,newton插值方法确定的插值多项式,并将函数y=f(x)的插值多项式和插值余项的图形画在同一坐标系中,观测节点变化对插值余项的影响。
常州大学2009~2010学年第 2 学期硕士生考试试题评分标准1. (10分)当x 充分大时, 试比较算上的差异?并叙述常见的防止误差的一些原则。
解:当x 充分大时,两个表达式在理论上恒等, 但其数值计算结果不同,前者会出现相近数相减,失去有效数位,降低计算结果精度的问题;后者避免了相近数相减的问题,尽可能地保证了计算结果的精度。
防止误差的几个基本原则主要有: 1) 防止大数“吃”小数;2) 避免除数绝对值远远小于被除数绝对值的除法; 3) 避免相近数相减;4) 避免使用不稳定的算法;注意简化计算步骤,减少运算次数;………… 5 分2. (15分)已知列表函数利用Newton 插值方法求()f x 的插值逼近多项式3()N x ,利用插值多项式近似计算(1.52)f 。
解:Newton 差商表: D =1.0000 -1.0000 -2.00003.00004.0000 3.0000 2.0000 -1.0000 -2.5000 -1.8333………… 5 分3() 1.8333^38.5000^28.6667 1.0000N x x x x =-+-+………… 5 分3(1.52)(1.52) 1.0268f N ≈=。
………… 5 分3. (10分)已知列表函数解:写出正规方程组42 5.1526 6.09a b a b +=⎧⎨+=⎩ ………… 5 分解上述正规方程组得0.9360,0.7030a b ==………… 5 分4. (15分)写出龙贝格(Romberg )方法的数值积分公式,并用龙贝格方法计算1sin 0x e dx ⎰,要求误差不超过210-。
解:龙贝格(Romberg )方法计算定积分()baf x dx ⎰的数值积分公式如下:211122221(),,22413316115156416363n n n i i i i n n nn n nn n nh b aT T f x h x x n S T T C S S R C C --=-=+=-==-=-=-∑,其中1[()()]2b aT f a f b -=+。
第三章 矩阵特征值与特征向量的计算--------学习小结一、 本章学习体会通过本章的学习,我们学到了四种矩阵特征值和特征向量的计算方法,分别是幂法、反幂法、Jacobi 方法和QR 方法。
四种方法各有其特点和适用范围。
幂法主要用于计算矩阵按模最大的特征值及其相应的特征向量;反幂法主要用于计算矩阵按模最小的特征值及其相应的特征向量;Jacobi 方法用于求实对称矩阵的全部特征值和特征向量的方法;QR 方法则适用于计算一般实矩阵的全部特征值,尤其适用于计算中小型实矩阵的全部特征值。
归结起来,这四种方法亦有其共同点,那就是都是用了迭代的方法来求矩阵的特征值和特征向量。
此外,用MA TLAB 自带的解法求解特征值和特征向量也非常快速,而且不用编辑函数建立m 文件。
其自带函数Eig 功能强大,即便得到结果是虚数也可以算出,并且结果自动正交化。
二、 本章知识梳理本章对于矩阵的特征值和特征向量的算法提出了新的思路,如幂法和反幂法、Jacobi 、QR 方法等。
本章的小结主要从方法的思想,以及一些定理展开。
2.1各种方法的运用范围1、幂法:主要用于计算矩阵按模最大的特征值和其相应的特征向量;2、反幂法:主要计算矩阵按模最小的特征值以及其相应的特征向量;3、Jacobi 方法:用于求实对称矩阵的全部特征值和特征向量的方法;4、QR 方法:适用于计算一般实矩阵的全部特征值,尤其适用于计算中小型实矩阵的全部特征值。
2.2各种方法的基本思想以及迭代公式 1.幂法幂法的基本思想:设n×n 实矩阵A 具有n 个线性无关的特征向量n x x x ,....,,21,其相应的特征值,,...,21n λλλ满足不等式n λλλλ≥≥>....321,其中iix i Ax λ=)...,3,2,1(n i =。
任取一n 维非零向量u 0,从u 0出发,按照如下的递推公式...)2,1(1===k Au u k k因n 维向量组n x x x ,....,,21线性无关,故对向量u 0必存在唯一的不全为0的数组a 1,a 2,...,a n ,使得n n x a x a x a u +++= (22110)])(...)([ (1)212211122211122110221n n n k n kn n k k n k n k k k k k k x a x a x a x a x a x a x A a x A a x A a u A u A Au u λλλλλλλλ+++=+++=+++=====--设a 1≠0,由上式可以看出,当k 充分大时有111x a u kkλ≈得迭代公式:9u A u k k = (1)从实际中来看,为了避免迭代向量u k 的模过大,(当11>λ)或过小(当11<λ),通常对u kj 进行归一化,使其范数等于1。
数值分析第3章3.1 填空题(1)当A 具有严格对角优势或具有对角线优势且矩阵不可约时,线性方程组b A x =用Jacobi 迭代法和Gauss-Seidel 迭代法均收敛;(2)当线性方程组的系数矩阵A 对称正定时,Gauss-Seidel 迭代收敛;(3)线性方程组迭代法收敛的充分必要条件是迭代矩阵的谱半径小于1;SOR 法收敛的必要条件是0<w<2;(4)用迭代法求解线性方程组,若)1(),(≥=q B q ρ时不收敛,q 接近0时收敛较快,q 接近1时收敛较慢; (5)A=⎪⎪⎭⎫⎝⎛2111;J B =⎪⎪⎭⎫ ⎝⎛--02110;s B =⎪⎪⎭⎫⎝⎛-21010;)22()(=J B ρ;)21()(=s B ρ 解:∵ A=⎪⎪⎭⎫⎝⎛2111,∴ L=⎪⎪⎭⎫⎝⎛0100, D=⎪⎪⎭⎫ ⎝⎛2001, U=⎪⎪⎭⎫⎝⎛0010 //P49D −1=1ad−bc d −b−c a∴⎪⎪⎭⎫⎝⎛--=⎪⎪⎭⎫ ⎝⎛⨯⎪⎪⎭⎫⎝⎛-=+-=-02110011021001)(1U L D B J //P49 B J 指Jacobi 迭代法的矩阵形式 又 ∵ )(U LD +=⎪⎪⎭⎫⎝⎛00211001−−−→−=-=22/12)2122)1r r r r r ⎪⎪⎭⎫ ⎝⎛-210101001=))((1U L D I -+//P50 B S 指Seidel 迭代法的矩阵形式,U L D B S 1)(-+-= ∴ ⎪⎪⎭⎫⎝⎛-=21010s B ∴ 22)(=J B ρ21)(=s B ρ //谱半径ρ A =max 1≤i≤n λi3.2 用Jacobi 迭代法和Gauss-Seidel 迭代法求解方程组,各分量第三位稳定即可。
(1) ⎝⎛012121⎪⎪⎪⎭⎫210⎪⎪⎪⎭⎫ ⎝⎛321x x x =⎪⎪⎪⎭⎫⎝⎛-453 解:利用Jacobi 迭代法,有⎪⎪⎪⎩⎪⎪⎪⎨⎧-=---=-=+++k k k k k kk x x x x x x x 21331122112122521212123,令0x =()T 0,0,0,代入有://P48k0 1 2 3 4 5 6 7 kx 1 0 1.5 2.75 3.625 4.25 4.6875 5 5.21875 k x 20 -2.5 -4.25 -5.5 -6.375 -7 -7.4375 -7.75 k x 323.254.1254.755.18755.55.71875k89101112131415k x 1 5.375 5.484375 5.5625 5.6171875 5.65625 5.68359375 5.703125 5.716796875 k x 2 -7.96875 -8.125 -8.234375 -8.3125 -8.3671875 -8.40625 -8.43359375 -8.453125 k x 35.8755.9843756.06256.11718756.156256.183593756.2031256.216796875k16 17 18 19 20 21 22 23 k x 1 5.726563 5.733398 5.738281 5.741699 5.744141 5.745850 5.747070 5.747925 k x 2 -8.466797 -8.476563 -8.483398 -8.488281 -8.491699 -8.494141 -8.495850 -8.497070 k x 36.2265636.2333986.2382816.2416996.2441416.2458506.2470706.247925由题,用Jacobi 迭代法进行第21次迭代后,前三位有效数字稳定,此时x=()T25.6,50.8,75.5-利用Gauss-Seidel 迭代法,得出⎪⎪⎪⎩⎪⎪⎪⎨⎧-=---=-=+++++1213211122112122521212123k k k k k kk x x x x x x x ,令0x =()T 0,0,0代入有//P50k0 1 2 3 4 5 6 7 k x 10 1.5 3.125 4.4375 5.09375 5.421875 5.585938 5.667969 kx 2 0 -3.25 -5.875 -7.1875 -7.84375 -8.171875 -8.335938 -8.417969 k x 33.6254.93755.593755.9218756.0859386.1679696.208984k8 9 10 11 12 13 kx 15.708984 5.729492 5.739746 5.744873 5.747437 5.748718 k x 2 -8.458984 -8.479492 -8.489746 -8.494873 -8.497437 -8.498718 k x 36.2294926.2397466.2448736.2474376.2487186.249359可见Gauss-Seidel 迭代法进行至12次迭代后稳定,x=()T 25.6,50.8,75.5-3.4 下面一些方程组的系数阵,试判断它们对Jacobi 迭代法,Gauss-Seidel 迭代法的收敛性。
一.(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- 。
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程序,并进行相关数值试验。
第三章作业
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 1
0.923879
0.707106
0.382683
(1)选择x0=0,x4=π/2的节点
y0=cosx0=1,y4=cosx4=0,可得
)
()
()()()(0101
1010
1x 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,可得
)
)(())(())(()
)(())(())(()(1202102
2101201
2010210
1x 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
,04
4∏
∑≠==--=i
j j j
i j i i x x x x y x P , 得
P3(x)=1+0.0031x-0.51542x +0.02423
x +0.02
844
x
4(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 16
1)拉格朗日差值
.选取
函数
],[),sin()cos(ππ-∈+=x x x y
x0=-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'); end
w=length(x0); n=w-1;
L=zeros(w,w); for k=1:n+1 V=1;
for j=1:n+1 if k~=j if
abs(x0(k)-x0(j))<eps
error('Divided by Zero,there are two nodes are the same'); end
V=conv(V,poly(x0(j)))/(x0(k)-x0(j));
end end L(k,:)=V;
end
C=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','The
original f(x) ','Error',0)
2)牛顿差值
.选取函数]5,5[,11
2
-∈+=
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:n
if abs(x0(k)-x0(k-j+1))<eps
error('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 end
C=D(n,n);
for k=(n-1):-1:1
C=conv(C,poly(x0(k))); m=length(C);
C(m)=C(m)+D(k,k); end
Y=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)。