数值分析考试

一、插值法(P50):1、己知函数在下列各点的值为0.20.40.60.81.0fM0.980.920.810.640.38用4次牛顿插值法对数据进行插值。functionf=lagfun(x)a=[0.2,0.4,0.6,0.8,1.0];b=[0.98,0.92,0.81,0.64,0.38];fori=l:5L⑴=1;forj=l:5ifJ=1L⑴=L(i)*(x-a(j))/(a(i)-a(j));endendendf=0；fori=l:5f=f+L(i)*b(i);end执行文件x0=[0.2,0.4,0.6,0.8,1.0];y0=[0.98,0.92,0.81,0.64,0.38];plot(xO,yO,'o')holdongridonfplot(Jlagfun，,[0,1]);holdonx=0:0.1:1;plot(x,newton(xO,yO,x),'r');

1.牛顿插值以及三次样条插值(第一个实验题)此题要求利用给定点，及给定点的函数值进行牛顿插值以及三次样条插值。a.牛顿插值要实现牛顿插值，要用到以下代码%调用格式：yi=Lagran_(x,y,xi)%x,y数组形式的数据豪functionfi=Lagran_(x,f,xi)fi=zeros(size(xi));npl=length(f);fori=l:nplz二ones(size(xi));forj=l:nplifi~=j,z=z.*(xi-x(j))/(x(i)-x(j));endendfi=fi+z*f(i);end二、曲线拟合(P95):16、观测物体的直线运动，得出以卜.数据:时间t/s00.91.93.03.95.0距离s.zm010305080110求运动方程。解：被观测物体的运动距离与运动时间大体为线性函数关系，从而选择线性方程s=a-^bt令C)=spn〃{l,f}则"；=6，|虹=53.63,(%,0)=14.7,(%,s)=280,(°")=1078,则法方程组为(614.7><280J4.753.63；J078,从而解得。=-7.855048'b=22.25376故物体运动方程为S=22.25376J7.855048MATLAB程序：x0=[00.91.933.95];y0=[010305080110];a=polyfit(x0.y0J)x=0:0.1:5;y=polyval(a,x);%计算拟合多项式在x的值plot(x0,y0,*,x,y)运行结果：n=22.2538-7.855012017、己知实验数据如下:192531384419.032349.073.397.8用最小二乘法求形如)，=。+bx2的经验公式，并计算均方误差。解：若s=a+bx2,贝ije=spo〃{l,r}则帆卜5,|础=7277699,(%,0)=5327,(/,%)=271.4,。*)=369321.5,则法方程组为(55327、/\a_<271.4、,53277277699,0_(369321.5,从而解得a=0.9726046=0.0500351故)，=0.9726046+0.0500351F均方误差为5=[£()，(易)一刀=0.1226j=O18、在某化学反应中，由实验得分解物质浓度与时间关系如卜.:时间t/s0510152025303540455055浓度y/(xl°T)0L272.162.863.443.874.154.374.514.584.624.64用最小二乘法求y=解：观察所给数据的特点，采用方程-by=ae，,(a,b>0)两边同时取对数，则]1bIny=lna--取中=sps"l,-L=bi>\x=--则S=a+i)x腐卜11，1就=°・°62321,(%,0)=-0.603975,(%,/)=-87.674095,(^,/)=5.032489,则法方程组为,11-0.603975Y/)』一87.674095、、—0.6039750.062321人歹J一(5.032489/从而解得=-7.5587812//=7.4961692因此=5.2151048/?=//=7.49616927.4961692.・・y=5.215104阮

X/(x)011/80.99739781/40.98961583/80.97672671/20.95885105/80.93615563/40.90885167/80.877192510.841409例三、对于函数/(a-)=—,给出n=8时X的函数表，使用复合梯形公式及复合辛普森公式计算枳分并估计误差以及比较其精度。三、数值积分与数值微分(例三、对于函数/(a-)=—,给出n=8时X的函数表，使用复合梯形公式及复合辛普森公式计算枳分并估计误差以及比较其精度。1复合梯形fhtfunctionT_n=fht(a,b,n)h=(b-a)/n;fork=0:nx(k+l)=a+k*h;ifx(k+l)==0x(k+l)=10A(-10);endendT_l=h/2*(fx(x(l))+fx(x(n+l)));fori=2:nF(i)=h*fx(x(i));endT_2=sum(F)T_n=T_l+T_2文件2fx.mfunctiony=fx(x)y=sin(x)./x;endmatlab中输入的»fht(0,1,8)0.83060.9457ans=0.9457>>formatlong»fht(0,1,8)T_2=0.830598927032208T_n=0.945690863582701ans=0.945690863582701>>formatshort»fht(0,1,8)复合桑普森复合S_P_S.m文件functionS_n=S_P_S(a,bzn)h=(b-a)/n;fork=0:nx(k+l)=a+k*h;x_k(k+l)=x(k+l)+l/2*h;if(x(k+l)==0)|(x_k(k+l)==0)x(k+l)=10A(-10);x_k(k+l)=10A(-10);endendS_l=(fx(x(l))+fx(x(n+l)))*(h/6);fori=2:nF1(i)=h/3*fx(x(i));endforj=l:nF_2(j)=2*h/3*fx(x_k(j));endS_2=sum(F_l)+sum(F_2);S_n=S_l+S_2;文件名fx.xfunctiony=fx(x)y=sin(x)./x;end在matlab输入»S_l=s_p_s(0,1.8)s_l=0.9461四、用牛顿法解方程(223):用牛顿法求解方程x3-x2-l=0在1.5附近的根(要求误差《10-6)。利用Newton法求方程的根。function[x_star,index,it]=Newton(fun,x,ep,it_max)%求解非线性方程的Newton法,其中%fun(x)—需要求根的函数，%第一个分量是函数值，第二个分量是导数值%x…初始点。%ep…精度，当l(x(k)-x(k-l)l<ep时，终止计算。省缺为le・5%it_max…最大迭代次数，省缺为100%x.star-当迭代成功时，输出方程的根，%当迭代失败时，输出最后的迭代值。%index―当index=l时，表明迭代成功，%当index=0时，表明迭代失败(迭代次数>=it.max)。%it-迭代次数。ifnargin<4it_max=100;endifnargin<3ep=le-5;endindex=0;k=l;whilek<=it_maxxl=x;f=feval(fun,x);ifabs(f(2))<epbreak;endx=x-f(l)/f(2);ifabs(x-xl)<epindex=l;break;enclk=k+l;endx_star=x;it=k;在Matlab窗口输入：fun=in!ine(,[xA3-x-13*xA2-l]t);[x_star,indexJt]=Newton(funJ.5)五、雅克比迭代法与高斯一赛德尔迭代法(187):给出线性方程组，写出计算格式判断两种方法的收敛性(谱半径)解方程(精度)需要多少次达到精度5x.4-2匚+x,=-12A—J设方程组+4丛+2人\=20,2x.-3xy+10x,=3考察用雅可比迭代法，高斯-赛德尔迭代法解此方程组的收敛性；用雅可比迭代法及高斯-赛德尔迭代法解此方程组，要求当<1。7时迭代终止。‘521、[解](a)由系数矩阵-142为严格对角占优矩阵可知，使用雅可比、高2-310X/斯-赛德尔迭代法求解此方程组均收敛。[精确解为X】=_4,丛=3叫=2](b)使用雅可比迭代法：x{k+l)=D~\L+U)x(k)+D~[b‘02!5(02P5555£-102X⑴+£2001X⑴+544422-30■331\71\713八、瓦3ioJ<10>使用高斯-赛德尔迭代法:"町=（。_匕）一以工⑶+（。_匕）-%\7O10O43O->72O-23/，‘\.1OJO-11-43一40。±402一:I+\-//12000ooO111OO-11-43一40120-设方程组设方程组（a）x.+04〕+0.4x.=1&■J<0.4i]+x2+O.&q=2;0.4i]+0.8x2+心=3X+2x2一2x3=1(b)x.+x.+x=1；

▲■J2i]+2x2+X3=10.40.40.40.4、10.8可知,0.81'1[解](a)由系数矩阵A=0.40.4(°-0.4—0.4)<0-0.4一0.4、B0=D~\L+U)=-0.40-0.8—-0.40一0.8〔J[-0.4-0.801/厂0.4一0.8°;由0.42\AI-B0\=0.40.40.420.80.8A=(2-0.8)(22+0.82-0.32)可知,=(2-0.8)(2+0.4-7048)(2+0.4+V048)=0p（5o）=O.4+VO48>b从而雅可比迭代法不收敛。<100'-1fo—0.4-0.4)G=(D-L)-p—0.41000-0.8.0.40.8〔°00\z/r100Y0-0.4-0.4、fo-0.4-0.4、—-0.41000—0.8=00.16-0.64-0.4-0.8100〔°0.160.8,\7\/20.40.4,由0.642-0.16-0.162-0.8="2—0"+0.1152)=可知,\AI-G\=00=2(Z-0.48+VO.1152)(2-0.48-Jo.1152)=0P(G)=0.48+VoJBT<1,从而高斯-塞德尔迭代法收敛。4(b)由系数矩阵A=122-2、1可知,17-22、22、B°=DSU)=i-10-1—-10-1厂2-20z厂2-20>\由22|2Z-B0|=1222-21=23=0可知，p(Bo)=02从而雅可比迭代法收敛。100、■?0-22)G=(D-LylU—11000-1,221/000/‘100)—22、-22)—-11000—-1—02-3¥-2ij1。0074-2/由2\AI-G\=002-2人一23=4(矛+8)=人以一27^)(/1+27^)=0可知,-42+2p(G)=2皿>1,从而高斯-塞德尔迭代法不收敛。程序代码function[x,i,G]=kb(Azbfmax,eps)D=diag(diag(A));ID=inv(D);J=ID*(-A+D);f=ID*b;xO=zeros(rank(A)f1);x=xO;G(l,:)=xO1;fori=l:maxk=x;x=J*x+f;G(i+1,:)=x!;n=norm((x-k),inf);ifn<=epsbreak;endendG.1在niatlab中输入A=[1023;2101;3110];»b=[141120]r;»inax=100;»eps=1.0e-006;»[x,i.G]=gs(A,b,max,eps)C=001.40000.82001.49800.78660.79291.68470.73600.78431.70080.73290.78331.70180.73280.78331.70180.73280.78331.70180.73280.78331.70180.73280.78331.7018-0.2000-0.30000.73280.78331.70180.73280.78331.7018-0.20000.04000.0560-03000-0.04000.0940六、常微分方程初值问题数值解法（318）:用经典四阶K-R方法解初值问题），，=-50.v+50x-+2x,0MxVl,

步长分别取h=0.1,0.025,0.01，计算各步长R数值解在x=0.4,0.6、0.8处的误差，指出误差值随步长的变化规律(己知真解y=^e~5Qx+x2)o这时经典R-K公式为解：取步长h二0.1,程序如下：%经典的四阶R-K方法clear;F='-50*y+50*x"2+2*x';a=0;b=l;h=0.1;n=(b~a)/0.1;X=a:O.l:b;Y二zeros(1,n+1);Y(l)=l/3;fori=l:nx=X(i);y=Y(i);Kl=h*eval(F);x=x+h/2;y=y+Kl/2;K2=h*eval(F);y=Y(i)+K2/2;K3=h*eval(F);x=X(i)+h;y=Y(i)+K3;K4=h*eval(F);Y(i+1)=Y(i)+(K1+2*K2+2*K3+K4)/6;end%准确值temp=L];f=dsolve(，Dyh50*y+50*x"2+2*x','y(0)=l/3*,'x');df=zeros(1,n+1);fori=l:n+ltemp=subs(f,'x',X(i));df(i)=double(vpa(temp));enddispC步长四阶经典R-K法准确值')；disp([X',Y',df']);运行结果:步长四阶经典R-K法准确值

00.000000000010000.000000000020000.000000000030000.000000000010000.000000000050000.000000000060000.000000000070000.000000000080000.000000000090000.00000000010000%画图观察结果figure;0.000000000033330.000000000460550.000000006306250.000000086404940.000001184363000.000016235451100.000222560671340.003050935427780.041823239217100.运行结果:步长四阶经典R-K法准确值00.000000000010000.000000000020000.000000000030000.000000000010000.000000000050000.000000000060000.000000000070000.000000000080000.000000000090000.00000000010000%画图观察结果figure;0.000000000033330.000000000460550.000000006306250.000000086404940.000001184363000.000016235451100.000222560671340.003050935427780.041823239217100.573326903478097.859356300837710.000000000033330.000000000001220.000000000001000.000000000009000.000000000016000.000000000025000.000000000036000.000000000049000.000000000061000.000000000081000.000000000