数值试验部分分析

1、第二章数值代数2P45-46实验题2(1)(3),5(A为12阶改为5阶)，8,10(n=300改为n=10);第三章迭代法25P71-72实验题2,3;习题2,12第四章数据建模34P106实验题1,2,4,8;习题3第五章数值微积分43P135实验题1(1)-(4)4,5,8;第六章数值分析及其MATLAE数验53P166实验题1(1),5,6第二章数值代数实验2求下列矩阵的行列式、逆、特征值、特征向量、各种范数、和条件数：(1)41-132-6J"31(1)>>A=41-1;32-6;1-53;a=det(A),B=inv(A),V,D=eig(A),t=eig(A

2、)%矩阵的行列式a、逆B、特征值t、特征向量Vnorm(A),norm(A,1),norm(A,inf)%分别为矩阵A的2,1,00-范数cond(A),cond(A,1),cond(A,inf)%分别为矩阵A的2,1,00-条件数a=+-94B=0.2553-0.02130.04260.1596-0.1383-0.22340.1809-0.2234-0.05320.0185-0.9009-0.3066-0.7693-0.1240-0.7248-0.6386-0.41580.6170D=-3.05270003.67600008.3766t=-3.05273.67608.3766ans=8.60

3、8910.000011.0000ans=3.74945.95745.7340(2)57651710876810957910>>A=5765;71087;68109;57910;a=det(A),B=inv(A),V,D=eig(A),t=eig(A)%矩阵的行列式a、逆B、特征值t、特征向量Vnorm(A),norm(A,1),norm(A,inf)%分别为矩阵A的2,1,00-范数cond(A),cond(A,1),cond(A,inf)%分别为矩阵A的2,1,00-条件数a=1.0000B=68.0000-41.0000-17.000010.0000-41.000025.000

4、010.0000-6.0000-17.000010.00005.0000-3.000010.0000-6.0000-3.00002.0000V=0.83040.09330.39630.3803-0.5016-0.30170.61490.5286-0.20860.7603-0.27160.55200.1237-0.5676-0.62540.5209D=0.010200.843100003.8581000030.2887t=0.01020.84313.858130.2887ans=30.288733.000033.0000ans=1.0e+03*2.98414.48804.4880实验5aj=1/

5、(i+j-1)。1,并分析结果(5)。Hilbert矩阵是著名的病态矩阵，n阶Hilbert矩阵定义为A=(aj),其中设A为5阶Hilbert矩阵，计算cond(A),A-1,norm(A-1A-E)及|A|A-1|的精度。再比较MATLA球解HIlbert矩阵及其逆函数hilb(5),invhilb>>A=hilb(5)A=1.00000.50000.33330.25000.20000.50000.33330.25000.20000.16670.33330.25000.20000.16670.14290.25000.20000.16670.14290.12500.20000.1

6、6670.14290.12500.1111>>cond(A)ans=4.7661e+005>>inv(A)ans=1.0e+005*0.0002-0.00300.0105-0.01400.0063-0.00300.0480-0.18900.2688-0.12600.0105-0.18900.7938-1.17600.5670-0.01400.2688-1.17601.7920-0.88200.0063-0.12600.5670-0.88200.4410>>norm(inv(A)*A-eye(5)ans=4.8038e-012>>det(A)*de

7、t(inv(A)-1ans=2.4158e-013>>det(hilb(5)ans=3.7493e-012>>det(inv(hilb(5)ans=2.6672e+011实验8分别用顺序Gauss消去法（程序2.1）和选列主元Gauss消去法（程序2.2）求解下列方程组，验证计算结果，并分析误差产生的原因：0.310-1559.145.291-6.1311.2912顺序Gauss消去法计算3r-xj一59.17一-12X246.7852X3-111-X4一-2,>>A=0.3*10A(-15)59.1431;5.291-6.13-12;11.2952;121

8、1;b=59.17;46.78;1;2;>>x=nagauss(A,b)a=1.0e+018*0.00000.00000.00000.00000.00000-1.0430-0.0529-0.0176-1.04360-2.2079-0.1120-0.0373-2.20900-0.1971-0.0100-0.0033-0.1972a=1.0e+018*0.00000.00000.00000.00000.00000-1.0430-0.0529-0.0176-1.043600-0.0000-0.000000000.00000a=1.0e+018*0.00000.00000.00000.00

9、000.00000-1.0430-0.0529-0.0176-1.043600-0.0000-0.000000000.00000x=23.68481.000500选列主元Gauss消去法计算>>A=0.3*10A(-15)59.1431;5.291-6.13-12;11,2952;1211;b=59.17;46.78;1;2;>>x=nagauss2(A,b)a=11.20009.00005.00002.00001.00000-10.3817-3.36211.055246.3076059.14003.00001.000059.170001.19640.55360.821

10、41.9107a=11.20009.00005.00002.00001.0000059.14003.00001.000059.170000-2.83541.230756.6946000.49290.80120.7137a=11.20009.00005.00002.00001.0000059.14003.00001.000059.170000-2.83541.230756.69460001.015110.5689x=3.84571.6095-15.476110.4113误差产生的原因由于第一列的主元素0.3X10-15的绝对值过于小，当它在消元过程作分母时把中间过程数值放大3.3X1015倍，使

11、中间结果“吃掉了”原始数据，从而造成数值不稳定。Gauss消去法针对以上问题，考虑选用绝对值大的数作为主元素。即选用选列主元实验10（追赶法的速度）考虑n阶三对角方程组（1）用选列主元Gauss消去法（程序2.2）求解;编写追赶法程序并求解；用矩阵除法求解；比较三种方法的计算时间和精度.用选列主元Gauss消去法（程序2.2）求解;>>A=2*ones（1,10）;%对角线元素B=1*ones(1,9);%对角线上方的元素，个数比A少一个C=1*ones(1,9);%对角线下方的元素，个数比A少一个D=diag(A)+diag(B,1)+diag(C,-1);b=3;4*ones(

12、8,1);3;x=nagauss2(D,b)第三章迭代法实验2用二分法(程序3.1)和Newton迭代法(程序3.2)求下列方程的正跟:xln(一x2-1x)-yx2-1-0.5x=0二分法functionx=nabisect(fname,a,b,e)ifnargin<4,e=1e-4;end;fa=feval(fname,a);fb=feval(fname,b);iffa*fb>0,error;endx=(a+b)/2while(b-a)>(2*e)fx=feval(fname,x);iffa*fx<0,b=x;fb=fx;elsea=x;fa=fx;endx=(a+

13、b)/2end>>fun=inline('x*log(xA2-1)A(1/2)+x)-(xA2-1)A(1/2)-0.5*x');x=nabisect(fun,2,3,0.05)x=2.5000x=2.2500x=2.1250x=2.0625x=2.0938x=2.0938Newton迭代法functionx=nanewton(fname,dfname,x0,e,N)ifnargin<5,N=500;endifnargin<4,e=1e-4;endx=x0;x0=x+2*e;k=0;whileabs(x0-x)>e&k<N,k=k+1

14、;x0=x;x=x0-feval(fname,x0)/feval(dfname,x0);disp(x)endifk=N,warning;end> >symsxydyy=x*log(xA2-1)A(1/2)+x)-(xA2-1)A(1/2)-0.5*x;dy=diff(y,x)dy=log(x+(xA2-1)A(1/2)-x/(xA2-1)A(1/2)+(x*(x/(xA2-1)A(1/2)+1)/(x+仅人2-1)A(1/2)> 1/2>>fun=inline('x*log(xA2-1)A(1/2)+x)-(xA2-1)A(1/2)-0.5*x')

15、;dfun=inline('log(x+国人2-1)A(1/2)-x/(xA2-1)A(1/2)+(x*(x/(xA2-1)A(1/2)+1)/(x+仅人2> 1)A(1/2)-1/2');nanewton(fun,dfun,2,0.5e-3)2.12012.11552.1155ans=2.1155验3已知函数f(x)=x4-2x在(-2,2)内有两个根。作图求取初值，用MATLAB命令求这两个根。> >fun=inline('xA4-2Ax','x')fun=Inlinefunction:fun(x)=*人4-2人*>

16、>fzero(fun,-22)Errorusingfzero(line274)Thefunctionvaluesattheintervalendpointsmustdifferinsign.> >fplot(fun,-22);gridon;> >fzero(fun,-20),fzero(fun,0,2)ans=-0.8613ans=1.2396> >fzero(fun,-1),fzero(fun,1)%参数x0也可以用一个点ans=-0.8613ans=1.2396> >x,f,h=fsolve(fun,-1),x,f,h=fsolve(f

17、un,1)%也可以用fsolve求解Equationsolved.fsolvecompletedbecausethevectoroffunctionvaluesisnearzeroasmeasuredbythedefaultvalueofthefunctiontolerance,andtheproblemappearsregularasmeasuredbythegradient.stoppingcriteriadetails>x=-0.8613f=3.7825e-12h=1Equationsolved.fsolvecompletedbecausethevectoroffunctionva

18、luesisnearzeroasmeasuredbythedefaultvalueofthefunctiontolerance,andtheproblemappearsregularasmeasuredbythegradient.<stoppingcriteriadetails>x=1.2396f=4.0459e-10h=1习题2.讨论下列求x32x-5=0在2,3上根的三种迭代格式的合理性与收敛性:(1)xk=3x015(2)xk=2(3)xk=(2xk,*5)x。-2并取x°=2.5,由一种收敛的迭代格式计算使结果有4位有效数字。33x-5.x-5斛：(1)易知f(x

19、)=0ux=,所以迭代格式合理。记g(x)=,由于当x=2,3,22、32一一g(x)=-x>6,定理3.1的条件不能满足。事实上，取Xo=2.5,计算得2x1=5.3125,x2=72.4664,x3=190272.0118,可知该迭代是发散的5.5(2)由于f(x)=0ux=力，所以迭代格式合理。记g(x)=力，由于当x=2,3,x-2x-2g(x)=-10x(x2-2)2g'(2)=5：1,定理3.1的条件不能满足11(3)由于f(x)=0=x=(2x+5)3,所以迭代格式合理。记g(x)=(2x+5)3,由于当xw2,3,g(x)=17-(2x+5)33111=3=L1,

20、g'(x)0,得g(x)单调上升，而g(2)=93和g(3)=1131口2,3,所以g:2,3T2,3。根据定理3.1,Vxow2,3,Xk=(2x+5)3,收敛于方程在(2,3】上的唯一解x”,且有误差估计k.Lx-x1-LXk-Xk4=-2|xk-Xy._lLi-1-LXi-X0-23k'：Xi-Xo1先验估计式Xk-x-1.wW,可用于计算前估算需迭代的次数。后验估计式一xkxkq<5,用于在2计算过程中判断数值解得精度。本例需要使结果有四位有效数字，则由X*12,3,知W=-x10'3o当2取X0=2.5时，得Xi=2.0801,从而由先验估计式得k之5。

21、进一步计算得x2=2.0924,x3=2.0942人=2.0945,由后验估计式知x*定x4=2.0945已经有4位有效数字。习题12取原点为初值，分别由Gauss-Seidel迭代格式与切=0.9时的逐次超松弛迭代格式解线性方程组(精度0.5x10-2)5xi+2x2+x3=-12-xi+3x2+2x3=172xi-3x2+4x3=-9雅可比迭代A=521;-132;2-34;b=-12,17,-9'x0=ones(3,1);x,n=jacobi(A,b,x0)x=-4.00003.00002.0000n=67Gauss-Seidel迭代格式程序functionx=nags(A,b,

22、x0,e,N)n=length(b);ifnargin<5,N=500;endifnargin<4,e=1e-4;endifnargin<3,x0=zeros(n,1);endx=x0;x0=x+2*e;k=0;Al=tril(A);iAl=inv(Al);whilenorm(x0-x,inf)>e&k<N,k=k+1;x0=x;x=-iAl*(A-Al)*x0+iAl*b;disp(x')endifk=N,warning;end>>A=521;-132;2-34;b=-1217-9'formatlong;x=nags(A,b,

23、0,0,0',0.5e-2)-2.4000000000000004.866666666666666-4.8666666666666672.311111111111111-3.7077777777777783.152962962962963-4.0549074074074073.0026234567901232.6000000000000001.9166666666666661.9686111111111112.029421296296296-4.006933641975309-3.988634387860083-4.005705454818244-3.998187232974680-4

24、.000296720122337x=-4.0002967201223373.0000708063995982.0002014648608672.9780745884773663.0124400291495202.9956667045324653.0008690851361263.0000708063995981.9870227623456792.0036472157921811.9996027558084701.9997454303394342.000201464860867逐次超松弛迭代格式程序functionx=nasor(A,b,omega,x0,e,N)n=length(b);ifna

25、rgin<6,N=500;endifnargin<5,e=1e-4;endifnargin<4,x0=zeros(n,1);endifnargin<3,omega=1.5;endx=x0;x0=x+2*e;k=0;L=tril(A,-1);U=triu(A,1);whilenorm(x0-x,inf)>e&k<N,k=k+1,x0=x;fori=1:nx1(i)=(b(i)-L(i,1:i-1)*x(1:i-1,1)-U(i,i+1:n)*x0(i+1:n,1)/A(i,i);x(i)=(1-omega)*x0(i)+omega*x1(i);endd

26、isp(x')endifk=N,warning;end>>A=521;-132;2-34;b=-1217-9'nasor(A,b,0.9,0,0,0',0.5e-2),formatshortk=1 -2.16004.45201.9521k=2 -4.33013.07492.1943k=3 -4.09502.86241.9693k=4-3.95443.01831.9888-4.00003.00852.0047k=6-4.00392.99692.0001k=7-3.99932.99981.9996ans=-3.99932.99981.9996第四章数据建模实验1

27、应用MATLAB由值和拟合命令解习题3、11、12、13、14。习题3不用计算器，利用函数f(x)=JX在100,121,144的值，分别用线性插值和抛物插值求JH5得近似值，并估计计算结果的误差。functionyy=nalagr(x,y,xx)m=length(x);n=length(y);ifm=n,error;ends=0;fori=1:nt=ones(1,length(xx);forj=1:nifj=i,t=t.*(xx-x(j)/(x(i)-x。)；endends=s+t*y(i);endyy=s;>>x1=100121;y1=1011;xx=115;yy1=nalag

28、r(x1,y1,xx)x2=100121144;y2=101112;yy2=nalagr(x2,y2,xx)fplot('xA(1/2)',100144);holdon;plot(x,y,'o',xx,yy1,1gA',xx,yy2,'rv');holdoff;yy1=10.714285714285715yy2=10.72275550536420112106100105110115120125130135140115=10.723805294763608115-yy1=10.723805294763608-10.7142857142857

29、15=0.009519580477892115-yy2=10.723805294763608-10.722755505364201=0.001049789399406习题11求三阶样条S(x)分别满足下列插值条件:x-1012f(x)-1010f(x)0-1x-101f(x)-101f'(»0-1>>clear;>>x=-1,0,1,2;>>y=-1,0,1,0;>>pp=csape(x,y,'complete',0,-1)>>pp.coefs>>fnplt(pp)pp=form:'

30、;pp'breaks:-1012coefs:3x4doublepieces:3order:4dim:1ans=-0.46671.46670-1.0000-0.60000.06671.533300.8667-1.7333-0.13331.0000s(x)=«-0.4667(x1)31.4667(x1)2-1,-1-x-032-0.6x0.0667x1.533x,0<x<1320.8667(x-1)-1.7333(x-1)一0.1333(x1)1,1三x三2(2)>>clear;x=-1,0,1;y=-1,0,1;pp=csape(x,y,'sec

31、ond',0,-1)pp.coefsfnplt(pp)pp=form:'pp'breaks:-101coefs:2x4doublepieces:2order:4dim:1ans=0.04170.0000-0.20830.12500.9583-1.00001.083303_s(x)="0.0417(x+1)3+0.9583(x+1)-1,-1<x<0-0.2083x3+0.1250x2+1.0833x,0<x<13习题12已知cos0=1,cos-=,cos-=0.5，由取小一乘法求cosx的拟合曲线，并比623较它们的拟合误差:(1)中

32、(x)=a+bx；(2)中(x)=a+bx2(3)(x)=aebx(1)>>fun=inline('c(1)+c(2)*x','c','x');x=0,pi/6,pi/3;y=1,3A(1/2)/2,0.5;c=lsqcurvefit(fun,0,0,x,y)norm(feval(fun,c,x)-y)A2Localminimumfound.Optimizationcompletedbecausethesizeofthegradientislessthanthedefaultvalueofthefunctiontolerance.&l

33、t;stoppingcriteriadetails>c=1.0387-0.4775ans=0.0090>>fun=inline('c(1)+c*x.AZl'cl'x');x=0,pi/6,pi/3;y=1,3A(1/2)/2,0.5;c=lsqcurvefit(fun,0,0,x,y)norm(feval(fun,c,x)-y)A2Localminimumfound.Optimizationcompletedbecausethesizeofthegradientislessthanthedefaultvalueofthefunctiontole

34、rance.<stoppingcriteriadetails>c=0.9959-0.4534ans=4.9565e-05(3)>>fun=inline('c(1)*exp(c(2)*x),c,x,);x=0,pi/6,pi/3;y=1,3A(1/2)/2,0.5;c=lsqcurvefit(fun,0,0,x,y)norm(feval(fun,c,x)-y)A2Localminimumpossible.lsqcurvefitstoppedbecausethefinalchangeinthesumofsquaresrelativetotsinitialvaluei

35、slessthanthedefaultvalueofthefunctiontolerance.<stoppingcriteriadetails>c=1.0367-0.5761ans=0.0157(1)的误差为0.0090(2)的误差为4.9565e-05(3)的误差为0.0157<(1)<(3)实验2利用Lagrange®值（程序4.1）作出Runge现象（例4.4）的图像。一、1例4.4设f(x)=分别讨论将-5,55等分与10等分后Lagrange插值的效果。1xfunctionyy=nalagr(x,y,xx)m=length(x);n=length(y

36、);ifm=n,error;ends=0;fori=1:nt=ones(1,length(xx);forj=1:nifj=i,t=t.*(xx-x(j)/(x(i)-x(j);endends=s+t*y(i);endyy=s;>>x=linspace(-5,5,6);y=1./(1+x.A2);xx=-5:0.001:5;yy=nalagr(x,y,xx);fplot('1./(1+x.A2)',-5,5);holdon;plot(xx,yy,'-');x=linspace(-5,5,11);y=1./(1+x.A2);xx=-5:0.001:5;y

37、y=nalagr(x,y,xx);plot(xx,yy,'-');2_Q§LJLJILJL.-5-4-3-2-1012345实验4用样条插值(程序4.2)解习题11(1)> >naspline(-1012,-1010,0,-1)ans=01.5333-0.1333-1.0000实验8解：> >t=0:1:24;> >T=15141414141516182022232528313231292725242220181716;> >p=polyfit(t,T,2)P=-0.09362.59438.4157> >ti

38、=0:0.1:24;> >Ti=polyval(p,ti);subplot(2,2,1);> >plot(t,T,'o',ti,Ti,'k');> >title('polyfit');> >p=polyfit(t,T,3);> >Ti=polyval(p,ti);subplot(2,2,2);> >plot(t,T,'o',ti,Ti,'k');polyfit4040第五章数值微积分实验1用MATLAB命令计算下列积分：12dx0x2:c.20

39、exp2xsinxdx222_11X!,(1)(-expdx(2)乂2冗工2!1(x'dx(4)1)>>clear;x=0:0.1:1;y=1/(2*pi)A(1/2)*exp(-x.A2/2);trapz(x,y)ans=0.3411(2)>>clear;x=eps:0.1:1;y=sin(x)./x;trapz(x,y)ans=0.9458(3)>>clear;x=0:0.1:1;y=x'(-x);trapz(x,y)ans=1.2868(4)>>clear;x=0:0.1:2*pi;y=exp(2*x).*sin(x).A2

40、;trapz(x,y)ans=3.5764e+04实验(4)(Simpson)积分法编制一个定步长Simpson法数值积分程序，并取n=10计算实验1(1),比较n=20的定步长梯形法(程序5.1)的精度。定步长Simpson%复化simpson积分公式fun=inline('1/(2*pi)A(1/2)*exp(-x.A2/2)');淅入函数a=0;b=1;%输入区间m=10;%输入等分段数h=(b-a)/m/2;%计算步长%计算以下积分值S=feval(fun,a)+feval(fun,b);fork=2:mS=S+2*feval(fun,a+2*(k-1)*h);endf

41、ork=1:mS=S+4*feval(fun,a+(2*k-1)*h);endS=S*h/3> >quad_4Simpson10S=0.3413> >formatlong> >SS=0.341344762887082定步长梯形法functiont=natrapz(fname,a,b,n)h=(b-a)/n;fa=feval(fname,a);fb=feval(fname,b);f=feval(fname,a+h:h:b-h+0.001*h);t=h*(0.5*(fa+fb)+sum(f);> >formatlong;natrapz(x)1/(2*

42、pi)A(1/2)*exp(-x.A2/2),eps,1,10),formatshort;ans=0.341143036535970> >formatlong;natrapz(inline('1/(2*pi)A(1/2)*exp(-x.A2/2)'),eps,1,20),formatshort;ans=0.341294331299304定步长Simpson0.341344762887082定步长梯形法0.341294331299304实验5编制一个变步长梯形法程序，解例5.14比较与Romberg积分法(程序5.3)的计算量变步长梯形法functionT2,k,n=

43、mvtrap(f,a,b,eps)h=b-a;T1=h*(feval(f,a)+feval(f,b)/2;T2=T1/2+(h/2)*feval(f,a+h/2);k=1;n=2;while(abs(T2-T1)>eps)h=h/2;T1=T2;T2=T1/2+(h/2)*sum(feval(f,a+h/2:h:b-h/2);k=k+1;n=2*n;end>>formatlongf=inline('sin(x)./x');T,k,n=mvtrap(f,eps,1,0.5e-6)T=0.946082974628235k=9n=512Romberg积分法funct

44、iont=naromberg(fname,a,b,e)ifnargin<4,e=1e-4;end;i=1;j=1;h=b-a;T(i,1)=h/2*(feval(fname,a)+feval(fname,b);T(i+1,1)=T(i,1)/2+sum(feval(fname,a+h/2:h:b-h/2+0.001*h)*h/2;T(i+1,j+1)=4Aj*T(i+1,j)/(4Aj-1)-T(i,j)/(4Aj-1);whileabs(T(i+1,i+1)-T(i,i)>ei=i+1;h=h/2;T(i+1,1)=T(i,1)/2+sum(feval(fname,a+h/2:h

45、:b-h/2+0.001*h)*h/2;forj=1:iT(i+1,j+1)=4Aj*T(i+1,j)/(4Aj-1)-T(i,j)/(4Aj-1);endendTt=T(i+1,j+1);>>formatlong;naromberg(inline('sin(x)./x'),eps,1,0.5e-6),formatshort;T=0.9207354924039480000.9397932848061770.946145882273587000.9445135216653890.9460869339517940.94608300406367400.9456908635

46、827010.9460833108884720.9460830693509170.946083070387222ans=0.946083070387222实验8分别用变步长梯形法，Romberg积分法(程序5.3),自适应Simpson积分法(程序5.4)计算下列积分，比较计算量(精度1e-8):42013xxexp-1.5xdx变步长梯形法functionT2,k,n=mvtrap(f,a,b,eps)h=b-a;T1=h*(feval(f,a)+feval(f,b)/2;T2=T1/2+(h/2)*feval(f,a+h/2);k=1;n=2;while(abs(T2-T1)>eps

47、)h=h/2;T1=T2;T2=T1/2+(h/2)*sum(feval(f,a+h/2:h:b-h/2);k=k+1;n=2*n;end>>formatlongf=inline('13.*(x-x.A2).*exp(-1.5*x)');T,k,n=mvtrap(f,0,4,1e-8)T=-1.548788373509369k=17131072Romberg积分法>>formatlong;naromberg(inline('13.*(x-x.A2).*exp(-1.5*x)'),0,4,1e-8),formatshort;Columns1

48、through6>>formatlong;naromberg(inline('13.*(x-x.A2).*exp(-1.5*x)'),0,4,1e-8),formatshort;T=-0.773370679119904-2.975612894688877-3.70969363321186800-2.354308177327338-2.147206604873492-2.0430408029842670-1.795077324392272-1.608667040080584-1.572764402427723-1.565299697656984-1.613549786

49、323728-1.553040606967546-1.549332178093344-1.548960238024544-1.565184666489712-1.549062959878374-1.548797783405762-1.5487893009504040000000000001.5488961617122600001.548788630608937-1.54878852549544400-1.552900423435264-1.548805675750448-1.548788523475253-1.548788376492229-1.548788372866903-1.548788

50、372614956-1.548788372577623-1.549817198012654-1.548789456205117-1.548788374902095-1.548788372543791-1.548788372528307-1.548788372527976-1.548788372527955-1.548788372527952ans=-1.548788372527952自适应Simpson法>>q,srm,err=naadapt(inline('13.*(x-x.A2).*exp(-1.5*x),),0,4,1e-8)q=-1.5488srm=err=3.25

51、24e-09> >formatlong> >qq=-1.548788372179998第六章实验1用oode45、ode23和ode113解下列微分方程:(Dy'=x+y,y(0)=1,0cx<3(要求输出x=1,2,3点的y值);> >odefun=inline('x+y','x','y');> >ode45(odefun,0,3,1;ode45(odefun,0,3,1;IError:Unbalancedorunexpectedparenthesisorbracket.Didyou

52、mean:> >ode45(odefun,0,3,1);> >t,y=ode45(odefun,0:3,1);t,yans=01.00001.00003.43662.000011.77813.000036.171140>>t,y=ode23(odefun,0:3,1);t,yans=01.00001.00003.43412.000011.76253.000036.1040ans=01.00001.00003.43642.000011.77713.000036.1683实验5分别用定步长的4阶经典Runge-Kutta格式(程序6.4,步长h=0.01)和变步

53、长的阶经典Runge-Kutta格式(程序6.5,要求精度10-5)解习题3习题3解初值y'=x+y,0MxW0.4,:y(0)=1,并与精确解y(x)=-x-1+2ex比较变步长的4阶经典Runge-Kutta格式>>clear;dyfun=inline('x+y');formatlong;x,y=nark4V(dyfun,0,0.4,1,1e-5,0.01);x,y,formatshort;ans=01.0000000000000000.3200000000000001.4342513150256010.4000000000000001.58364482

54、5768755定步长的4阶经典Runge-Kutta格式>>clear;dyfun=inline('x+y');x,y=nark4(dyfun,0,0.4,1,0.01);x,yans=实验6考虑Euler格式(程序6.1),隐式Euler格式(程序6.2),中点Euler格式，改进Euler格式(程序6.3),4阶经典Runge-Kutta格式(程序6.4)求解微分方程y二-50y取步长h=0.1,0.05,0.025,0.001,试验其稳定性。Euler格式(程序6.1),h=0.1>>clear;dyfun=inline('0*x-50*y');x,y=naeuler(dyfun,0,0.2,1,0.1);x,yans=01.00000.1000-4.00000.200016.0000h=0.05>>x,y=naeuler(dyfun,0,0.2,