数值分析(第五版)计算实习题第四章作业

第四章：1、（1）：复合梯形建立m文件：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));输入：symsxf=inline('sqrt(x).*log(x);');natrapz(f,eps,1,10)输出：ans=-0.417062831779470输入：symsxf=inline('sqrt(x).*log(x);');natrapz(f,eps,1,100)输出：ans=-0.443117908008157输入：symsxf=inline('sqrt(x).*log(x);');natrapz(f,eps,1,1000)输出：ans=-0.444387538997162复合辛普森建立m文件：functiont=comsimpson(fname,a,b,n)h=(b-a)/n;fa=feval(fname,a);fb=feval(fname,b);f1=feval(fname,a+h:h:b-h+0.001*h);f2=feval(fname,a+h/2:h:b-h+0.001*h);t=h/6*(fa+fb+2*sum(f1)+4*sum(f2));输入：symsxf=inline('sqrt(x).*log(x);');formatlong；comsimpson(f,eps,1,10)输出：ans=-0.435297890074689输入：symsxf=inline('sqrt(x).*log(x);');comsimpson(f,eps,1,100)输出：ans=-0.444161178415673输入：symsxf=inline('sqrt(x).*log(x);');comsimpson(f,eps,1,1000)输出：ans=-0.444434117614180（2）龙贝格建立m文件：function[RT,R,wugu,h]=Romberg(fun,a,b,wucha,m)%RT是龙贝格积分表%R是数值积分值%wugu是误差估计%h是最小步长%fun是被积函数%ab是积分下、上限%m是龙贝格积分表中行最大数目%wucha是两次相邻迭代值的绝对误差限n=1;h=b-a;wugu=1;x=a;k=0;RT=zeros(4,4);RT(1,1)=h*(feval(fun,a)+feval(fun,b))/2;while((wuguwucha)&(km)|(k4))k=k+1;h=h/2;s=0;forj=1:nx=a+h*(2*j-1);s=s+feval(fun,x);endRT(k+1,1)=RT(k,1)/2+h*s;n=2*n;fori=1:kRT(k+1,i+1)=((4^i)*RT(k+1,i)-RT(k,i))/(4^i-1);endwugu=abs(RT(k+1,k)-RT(k+1,k+1));endR=RT(k+1,k+1);输入：fun=inline('sqrt(x).*log(x)');[RT,R,wugu,h]=Romberg(fun,eps,1,1e-5,13)输出：RT=1至5列-0.0000002685461450000-0.245064670140209-0.326752804004897000-0.358104125949240-0.395783944552250-0.40038602058874100-0.408090073087781-0.424752055467295-0.426683262861631-0.4271006794056450-0.429474601629505-0.436602777810080-0.437392825966266-0.437562819031419-0.437603847029951-0.438389494461832-0.441361125405941-0.441678348578999-0.441746372747455-0.4417627788404596列00000-0.441766844267449R=-0.441766844267449wugu=4.065426989774412e-06h=0.031250000000000（3）自适应辛普森输入：f=inline('sqrt(x).*log(x)');q=quad(f,0,1,1e-4)输出：q=-0.4439755729517282.（1）复合辛普森建立m文件functionq=combinesimpson2(F,x0,a,b,n)%复合Simpson多元求积公式%F—被积函数%x0—被积函数自变量%[a,b]积分区间%n—区间份数x=linspace(a,b,n+1);q=0;fork=1:nq=q+subs(F,x0,x(k))+4*subs(F,x0,(x(k)+x(k+1))/2)+subs(F,x0,x(k+1));endq=q*(b-a)/n/6;输入：clearsymsxy;F=exp(-x.*y);s=combinesimpson2(combinesimpson2(F,'x',0,1,4),'y',0,1,4)输出：s=exp(-1)/576+exp(-1/2)/144+exp(-1/4)/72+exp(-3/4)/144+exp(-1/8)/36+exp(-3/8)/36+exp(-5/8)/72+exp(-7/8)/72+(5*exp(-1/16))/144+exp(-3/16)/24+exp(-5/16)/36+exp(-7/16)/36+exp(-9/16)/144+exp(-1/32)/36+exp(-3/32)/18+exp(-5/32)/36+exp(-7/32)/36+exp(-9/32)/36+exp(-15/32)/36+exp(-21/32)/36+exp(-1/64)/36+exp(-3/64)/18+exp(-5/64)/18+exp(-7/64)/18+exp(-9/64)/36+exp(-15/64)/18+exp(-21/64)/18+exp(-25/64)/36+exp(-35/64)/18+exp(-49/64)/36+47/576double(s)ans=0.796599967946203高斯求积公式functionq=gaussquad(F,x0,a,b,n)%Gauss求积公式%F—被积函数%x0—被积函数自变量%[a,b]积分区间%n—节点个数symst;F=subs(F,x0,(b-a)/2*t+(a+b)/2);[x,A]=gausspoints(n);q=(b-a)/2*sum(A.*subs(F,t,x));输入：clearsymsxy;F=exp(-x.*y);s=gaussquad(gaussquad(F,x,0,1,4),y,0,1,4)输出：s=0.7966（2）复合辛普森输入：symsxy;f=exp(-x.*y);s=combinesimpson2(combinesimpson2(f,y,0,sqrt(1-x^2),4),x,0,1,4)输出：s=(3^(1/2)*(exp(-3^(1/2)/4)+2*exp(-3^(1/2)/8)+2*exp(-3^(1/2)/16)+2*exp(-(3*3^(1/2))/16)+4*exp(-3^(1/2)/32)+4*exp(-(3*3^(1/2))/32)+4*exp(-(5*3^(1/2))/32)+4*exp(-(7*3^(1/2))/32)+1))/576+(7^(1/2)*(exp(-(3*7^(1/2))/16)+2*exp(-(3*7^(1/2))/32)+2*exp(-(3*7^(1/2))/64)+2*exp(-(9*7^(1/2))/64)+4*exp(-(3*7^(1/2))/128)+4*exp(-(9*7^(1/2))/128)+4*exp(-(15*7^(1/2))/128)+4*exp(-(21*7^(1/2))/128)+1))/1152+(15^(1/2)*(exp(-15^(1/2)/16)+2*exp(-15^(1/2)/32)+2*exp(-15^(1/2)/64)+2*exp(-(3*15^(1/2))/64)+4*exp(-15^(1/2)/128)+4*exp(-(3*15^(1/2))/128)+4*exp(-(5*15^(1/2))/128)+4*exp(-(7*15^(1/2))/128)+1))/1152+(15^(1/2)*(exp(-(7*15^(1/2))/64)+2*exp(-(7*15^(1/2))/128)+2*exp(-(7*15^(1/2))/256)+2*exp(-(21*15^(1/2))/256)+4*exp(-(7*15^(1/2))/512)+4*exp(-(21*15^(1/2))/512)+4*exp(-(35*15^(1/2))/512)+4*exp(-(49*15^(1/2))/512)+1))/1152+(39^(1/2)*(exp(-(5*39^(1/2))/64)+2*exp(-(5*39^(1/2))/128)+2*exp(-(5*39^(1/2))/256)+2*exp(-(15*39^(1/2))/256)+4*exp(-(5*39^(1/2))/512)+4*exp(-(15*39^(1/2))/512)+4*exp(-(25*39^(1/2))/512)+4*exp(-(35*39^(1/2))/512)+1))/1152+(55^(1/2)*(exp(-(3*55^(1/2))/64)+2*exp(-(3*55^(1/2))/128)+2*exp(-(3*55^(1/2))/256)+2*exp(-(9*55^(1/2))/256)+4*exp(-(3*55^(1/2))/512)+4*exp(-(9*55^(1/2))/512)+4*exp(-(15*55^(1/2))/512)+4*exp(-(21*55^(1/2))/512)+1))/1152+(63^(1/2)*(exp(-63^(1/2)/64)+2*exp(-63^(1/2)/128)+2*exp(-63^(1/2)/256)+2*exp(-(3*63^(1/2))/256)+4*exp(-63^(1/2)/512)+4*exp(-(3*63^(1/2))/512)+4*exp(-(5*63^(1/2))/512)+4*exp(-(7*63^(1/2))/512)+1))/1152+1/24double(s)ans=0.670113633359095