数值分析实验

1、信计、数应专业10级数值方法计算实习 要求：1、用Matlab语言或你熟悉的其他语言编写程序，使之尽可能具有通用性；2、根据上机计算实践，对所使用的数值方法的特点、性质、有效性、误差和收敛性等方面进行必要的讨论和分析；3、完成计算后写出实验报告，内容包括：课题名称、解决的问题、采用的数值方法、算法程序、数值结果、对实验结果的讨论和分析等；4、特别说明：严禁抄袭，否则一经发现，所有雷同实验报告最多评为及格。一、在区间-1,1上分别取用两组等距节点对龙格函数作多项式插值及三次样条插值，对每个值，分别画出插值函数即的图形。解答：（1） 多项式插值function y1=f1(x)y1=1./(1+2

2、5.*x.*x)function f=L(x,y,z)syms t;f=0;for(i=1:n) l=y(i); for(j=1:i-1) l=l*(t-x(j)/(x(i)-x(j); end; for(j=i+1:n) l=l*(t-x(j)/(x(i)-x(j); end; f=f+l; simplify(f); if(i=n) if(nargin=3) f=subs(f,'t',z); else f=collect(f); f=vpa(f,6); end endendMATLAB实现：>> x=linspace(-1,1,10)x =-1.0000 -0.7

3、778 -0.5556 -0.3333 -0.1111 0.1111 0.3333 0.5556 0.7778 1.0000>> y1=f1(x)y1 =0.0385 0.0620 0.1147 0.2647 0.7642 0.7642 0.2647 0.1147 0.0620 0.0385y1 =0.0385 0.0620 0.1147 0.2647 0.7642 0.7642 0.2647 0.1147 0.0620 0.0385>> L(x,y1)ans =30.7285*t4-8.26092*t2-44.9155*t6+21.6248*t8+.861538>

4、;> x=linspace(-1,1,20)x = -1.0000 -0.8947 -0.7895 -0.6842 -0.5789 -0.4737 -0.3684 -0.2632 -0.1579 -0.0526 0.0526 0.1579 0.2632 0.3684 0.4737 0.5789 0.6842 0.7895 0.8947 1.0000>> y1=f1(x)y1 =0.0385 0.0476 0.0603 0.0787 0.1066 0.1513 0.2276 0.3661 0.6160 0.9352 0.9352 0.6160 0.3661 0.2276 0.1

5、513 0.1066 0.0787 0.0603 0.0476 0.0385y1 =0.0385 0.0476 0.0603 0.0787 0.1066 0.1513 0.2276 0.3661 0.6160 0.9352 0.9352 0.6160 0.3661 0.2276 0.1513 0.1066 0.0787 0.0603 0.0476 0.0385>> L(x,y1)Ans= = 121019.*t12-146791.*t14+95604.8*t16-3055.32*t6+17172.5*t8+327.726*t4+.992681-58585.0*t10-25671.2

6、*t18-21.6235*t2>>q=polyfit(x,y1,12)>> z=polyval(q,x)>> plot(x,y1,':o',x,z,'-*') （2） 三次样条插值二、已知Wilson矩阵，且向量，则方程组有准确解。用Matlab 内部函数求，的所有特征值和；解答：>> A=10,7,8,7;7,5,6,5;8,6,10,9;7,5,9,10A = 10 7 8 7 7 5 6 5 8 6 10 9 7 5 9 10>> det(A)ans =1>> B=eig(A)B =

7、0.0102 0.8431 3.8581 30.2887>> cond=30.2887/0.0102cond =2.9695e+003令，解方程组，并求出向量和，从理论结果和实际计算结果两方面分析方程组解的相对误差与的相对误差的关系；解答：C=10,7,8.1,7.2;7.08,5.04,6,5;8,5.98,9.89,9;6.99,4.99,9,9.98C = 10.0000 7.0000 8.1000 7.2000 7.0800 5.0400 6.0000 5.0000 8.0000 5.9800 9.8900 9.0000 6.9900 4.9900 9.0000 9.980

8、0>> D=32;23;33;31D = 32 23 33 31>> rank(C)ans = 4>> rank(C,D)ans = 4>> x=CDx = -81.0000 137.0000 -34.0000 22.0000>> x=C/D? Error using => mrdivideMatrix dimensions must agree.>> e? Undefined function or variable 'e'.>> e=1;1;1;1e = 1 1 1 1>>

9、 y=x-ey = -82.0000 136.0000 -35.0000 21.0000>> norm(y)ans = 163.9695>> norm(e)ans = 2>> 163.9695/2ans = 81.9848>> a=C-Aa = 0 0 0.1000 0.2000 0.0800 0.0400 0 0 0 -0.0200 -0.1100 0 -0.0100 -0.0100 0 -0.0200>> norm(a)ans = 0.2308>> norm(A)ans = 30.2887>> 0.2308

10、/30.2887ans = 0.0076>> F=inv(A)F = 25.0000 -41.0000 10.0000 -6.0000 -41.0000 68.0000 -17.0000 10.0000 10.0000 -17.0000 5.0000 -3.0000 -6.0000 10.0000 -3.0000 2.0000>> norm(F)*norm(a)*norm(x)ans = 3.7347e+003>> E=F*aE = -3.2200 -1.7800 1.4000 5.1200 5.3400 2.9600 -2.2300 -8.4000 -1.

11、3300 -0.7500 0.4500 2.0600 0.7800 0.4400 -0.2700 -1.2400>> norm(E)ans = 12.7943>> 1-12.7943ans = -11.7943>> (3.7347e+003)/(-11.7943)ans = -316.6530所以>再改变扰动矩阵（其元素的绝对值不超过0.005），重复第2问。解答：>> c=0.001 0.005 0.001 0.003;0.001 0.002 0.002 0.001;0.002 -0.004 -0.003 0.005;-0.001 -0.

12、003 0 -0.004c = 0.0010 0.0050 0.0010 0.0030 0.0010 0.0020 0.0020 0.0010 0.0020 -0.0040 -0.0030 0.0050 -0.0010 -0.0030 0 -0.0040>> d=c+Ad = 10.0010 7.0050 8.0010 7.0030 7.0010 5.0020 6.0020 5.0010 8.0020 5.9960 9.9970 9.0050 6.9990 4.9970 9.0000 9.9960>> f=C Df = 10.0000 7.0000 8.1000 7.2

13、000 32.0000 7.0800 5.0400 6.0000 5.0000 23.0000 8.0000 5.9800 9.8900 9.0000 33.0000 6.9900 4.9900 9.0000 9.9800 31.0000>> rank(d)ans = 4>> rank(f)ans = 4>> x2=dDx2 = 0.9427 1.0904 0.9761 1.0172>> x2-eans = -0.0573 0.0904 -0.0239 0.0172>> norm(y)ans = 163.9695>> no

14、rm(e)ans = 2>> 163.9695/2ans = 81.9848>> norm(x2)ans = 2.0162>> 2.0162/2ans = 1.0081>> norm(c)ans = 0.0081>> norm(a)/norm(A)ans = 0.0076>> norm(c)/norm(A)ans = 2.6854e-004>> norm(F)*norm(c)*norm(x2)ans = 1.6157>> G=F*cG = 0.0100 0.0210 -0.0870 0.1080 -

15、0.0170 -0.0310 0.1460 -0.1800 0.0060 0.0050 -0.0390 0.0500 -0.0040 -0.0040 0.0230 -0.0310>> 1-norm(G)ans = 0.7166>> 1.6157/0.7166ans = 2.2547所以<三、解三对角线性方程组的追赶法及其应用编写解三对角线性方程组的追赶法的通用程序，并应用于方程组，检验程序的正确性；（解为）解答：function x=zgan(A,b)n = rank(A);for(i=1:n) if(A(i,i)=0) disp('Error£

16、¡'); return; endend;d = ones(n,1);a = ones(n-1,1);c = ones(n-1);for(i=1:n-1) a(i,1)=A(i+1,i); c(i,1)=A(i,i+1); d(i,1)=A(i,i);endd(n,1) = A(n,n);for(i=2:n) d(i,1)=d(i,1) - (a(i-1,1)/d(i-1,1)*c(i-1,1); b(i,1)=b(i,1) - (a(i-1,1)/d(i-1,1)*b(i-1,1);endx(n,1) = b(n,1)/d(n,1);for(i=(n-1):-1:1) x(i

17、,1) = (b(i,1)-c(i,1)*x(i+1,1)/d(i,1);endMATLAB调用：>> A=2,-1,0,0,0;-1,2,-1,0,0;0,-1,2,-1,0;0,0,-1,2,-1;0,0,0,-1,2A = 2 -1 0 0 0 -1 2 -1 0 0 0 -1 2 -1 0 0 0 -1 2 -1 0 0 0 -1 2>> b=1;0;0;0;0b = 1 0 0 0 0>> zgan(A,b)ans = 0.8333 0.6667 0.5000 0.3333 0.1667>> format>> x=5/6;2

18、/3;1/2;1/3;1/6x = 0.8333 0.6667 0.5000 0.3333 0.1667两种方法球的结果相等！得证！求微分方程边值问题的数值解（取步长），并与精确解比较（精确解为）。说明：离散化微分方程时，解答：h=1/128;u0=-2;upi=exp(pi);for x=h:h:pi-h u1=1/(1+x); u2=1/(1+x+h); s=(u2-2*u1+u0)/(h*h); t=(u2-u0)/(2*h); u=exp(x)-3*sin(x)-s+t; disp(u); d=(2*h*h*exp(x)-6*h*h*sin(x)-(2*h*h-4)*u2-(2-h)*

19、u1)/(2+h); uo=u1; u1=u2; u2=d;end MATLAB调用：>> weifen 0.9338e+004 0.9211e+004 0.9085e+004 0.8962e+004 0.8840e+004 0.3471e+004 0.3418e+004。function y=f(x)y=1./(1+x)>> x=h:h:10*hx =0.00780.01560.02340.03130.03910.04690.05470.06250.07030.0781>> y=f(x)y =0.99220.98460.97710.96970.96240.

20、95520.94810.94120.93430.9275四、公元1225年，比萨的数学家Leonardo研究了方程，得到一个根，没有人知道他用什么方法得到这个值。对于这个方程，分别用下列方法：解答：迭代法；functionx0,k,h,x=dd(g,x0,s,max)X(1)=x0;for k=2:max X(k)=feval('g',X(k-1); k,h=abs(X(k)-X(k-1)x=X(k);disp(x);if(h<s), break;endif k=max disp('out range !'); end end MATLAB输入：>&

21、gt; dd('g',1,10(-9),15)k = 2h = 0.5385x= 1.5385k = 3h = 0.2434x= 1.2950k = 4h = 0.1068x= 1.4018k = 5h = 0.0476x= 1.3542k = 6h = 0.0211x= 1.3753k = 7h = 0.0094x= 1.3659k = 8h = 0.0042x= 1.3701k = 9h = 0.0018x= 1.3682k = 10h = 8.1879e-004x= 1.3691k = 11h = 3.6341e-004x= 1.3687k = 12h = 1.6129e

22、-004x= 1.3689k = 13h = 7.1586e-005x= 1.3688k = 14h = 3.1772e-005x= 1.3688k = 15h = 1.4101e-005x= 1.3688out range !ans = 1迭代法；>> dd('h',1,10(-10),30)k = 2h = 0.5385x= 1.5385k = 3h = 0.2434x= 1.2950k = 4h = 0.1068x= 1.4018k = 5h = 0.0476x= 1.3542k = 6h = 0.0211x= 1.3753k = 7h = 0.0094x=

23、1.3659k = 8h = 0.0042x= 1.3701k = 9h = 0.0018x= 1.3682k = 10h = 8.1879e-004x= 1.3691k = 11h = 3.6341e-004x= 1.3687k = 12h = 1.6129e-004x= 1.3689k = 13h = 7.1586e-005x= 1.3688k = 14h = 3.1772e-005x= 1.3688k = 15h = 1.4101e-005x= 1.3688k = 16h = 6.2585e-006x= 1.3688k = 17h = 2.7777e-006x= 1.3688k = 18

24、h = 1.2328e-006x= 1.3688k = 19h = 5.4716e-007x= 1.3688k = 20h = 2.4285e-007x= 1.3688k = 21h = 1.0778e-007x= 1.3688k = 22h = 4.7837e-008x= 1.3688k = 23h = 2.1231e-008x= 1.3688k = 24h = 9.4231e-009x= 1.3688k = 25h = 4.1822e-009x= 1.3688k = 26h = 1.8562e-009x= 1.3688k = 27h = 8.2383e-010x= 1.3688k = 28

25、h = 3.6564e-010x= 1.3688k = 29h = 1.6228e-010x= 1.3688k = 30h = 7.2025e-011x= 1.3688ans = 1对的Steffensen加速方法；i=2;x0=1;f=inline('20/(x2+2*x+10)');y=f(x0);z=f(y);x1=x0-(y-x0).2/(z-2*y+x0);S.result=x0;x1;while abs(x1-x0)>=1e-9 x0=x1; y=f(x0); z=f(y); x1=x0-(y-x0).2./(z-2.*y+x0); i=i+1; S.resu

26、lt(i)=x1;endS.step=(0:i-1)'fprintf('diedaishu:t%dn',i-1);for j=1:i fprintf('%10dt',S.step(j);end输入：Steffensen：4 0 1 2 3 4 对的Steffensen加速方法；MATLAB输入：Steffensen：Newton法。求方程的根（可取），计算到Leonardo所得到的准确度。解答：functionp1,h,k,y=newdon(f,wei,p0,r,max)p0,feval('f',p0)for k=1:max p1=p0-

27、feval('f',p0)/feval('wei',p0); h=abs(p1-p0); p0=p1; p1,h,k,y=feval('f',p1) if(h<r)|(y=0), break,endendMATLAB输入：>> newdon('f','wei',1,10(-9),10)p0 = 1ans = -7p1 = 1.4118h = 0.4118k = 1y = 0.9176p1 = 1.3693h = 0.0424k = 2y = 0.0111p1 = 1.3688h = 5.2828e

28、-004k = 3y = 1.7045e-006p1 = 1.3688h = 8.0796e-008k = 4y = 3.9080e-014p1 = 1.3688h = 1.7764e-015k = 5y = 0ans = 1.3688五、用不同的数值方法计算积分的近似值，其中 取不同的步长，分别用复合梯形公式和复合辛普森公式计算积分，比较两个公式的计算效果，是否存在一个最小的，使得精度不能再被改善？解答：一：复合梯形公式：function s=f1(f,a,b,n)h=(b-a)/n;s=0;for k=1:(n-1) x=a+k*h; s=s+feval('f',x);en

29、ds=h*(feval('f',a)+2*h*s+feval('f',b)/2;输入：>> f1('f2',1,3,10)ans = 11.3280>> f1('f2',1,3,20)ans = 5.8960>> f1('f2',1,3,30)ans = 3.9830>> f1('f2',1,3,40)ans = 3.0070>> f1('f2',1,3,50)ans = 2.4151二：复合辛普森公式：function s=f3(f,a,b,n)h=(b-a)/n;i=0;j=0;for k=0:n-1 x=a+h/2; i=i+feval('f',x);endfor k=1:n-1 x=a+h; j=j+feval('f',x);ends=h*(feval('f',a)+feval('f