数值分析(王能超)matlab上机作业

引论习题题目：用二分法求方程x^3-3x-1=0在[1,2]内的近似根，要求误差不超过10-3程序：function[root,n]=dichotomy(a,b,eps)%dichotomy：二分法求根函数%f:带求解方程；[a,b]求解区间；n为二分次数ifnargin2disp('输入变量不够');return;elseifnargin3disp('输入变量过多');return;elsenargin==2eps=1.0e-3;endendifabdisp('区间输入错误')returnendc=(a+b)/2;iff(c)==0y=cn=1returnendn=0;whileabs(b-a)epsiff(a)*f(c)0a=c;elseb=c;endc=(a+b)/2;n=n+1;endroot=c;end运算结果：定义函数f如下：function[yy]=f(x)%定义了测试函数%x为输入变量yy=x^3-3*x-1;end在commondwindow中运行结果如下：[root,n]=dichotomy(1,2,1.0e-3)root=1.8794n=10习题一题目一给出概率积分202xxyedx的数据表i0123xi0.460.470.480.49yi0.48465550.49374520.50274980.5116683用二次插值计算，试问：（1）当x=0.472时该积分值等于多少？（2）当x为何值时积分值等于0.5？程序：function[y,t]=interpolation(X,Y,x)%拉格朗日插值函数%X为自变量数组，Y为函数值数组，x为插值点，t为插值基函数n=length(X);y=0;t=zeros(1,n);fori=1:nt(i)=1;forj=1:nifi==jcontinueendt(i)=t(i)*(x-X(j))/(X(i)-X(j));endy=y+t(i)*Y(i);endend运算结果：在commondwindow中运行结果如下：（1）当x=0.472时该积分值等于多少？clearX=[0.460.470.480.49];Y=[0.48465550.49374520.50274980.5116683];x=0.472;[y,t]=interpolation(X,Y,x)y=0.4956t=-0.04800.86400.2160-0.0320（2）当x为何值时积分值等于0.5？clearX=[0.460.470.480.49];Y=[0.48465550.49374520.50274980.5116683];y=0.5;[x,t]=interpolation(Y,X,y)x=0.4769t=-0.04520.33560.7707-0.0611题目二构造适合下列数据表的三次插值样条函数S（x）:x-1013y-1135y'61程序：function[y,m]=spline(X,Y,x,m1,mn)%样条插值函数%X，Y为输入的自变量、函数值数组；%x为插值点(数组)，y为插值结果(数组)，m为各插值节点的一阶导数值%m1为X(1)的一阶导数;mn为X(n)的一阶导数n=length(X);alpha=zeros(n-2,1);beta=zeros(n-2,1);fori=1:(n-2)alpha(i)=(X(i+1)-X(i))/(X(i+1)-X(i)+X(i+2)-X(i+1));beta(i)=3*((1-alpha(i))*(Y(i+1)-Y(i))/(X(i+1)-X(i))+alpha(i)*(Y(i+2)-Y(i+1))/(X(i+2)-X(i+1)));endm=zeros(n,1);m(1)=m1;m(n)=mn;A=zeros((n-2),(n-2));B=zeros((n-2),1);forj=1:(n-2)A(j,j)=2;endfork=2:(n-2)A(k,(k-1))=1-alpha(k);endforl=1:(n-3)A(l,(l+1))=alpha(l);endB(1)=beta(1)-(1-alpha(1))*m1;B(n-2)=beta(n-2)-alpha(n-2)*mn;forp=2:(n-3)B(p)=beta(p);endm(2:(n-1))=A\B;s=length(x);y=zeros(1,s);fort=1:sr=1;forq=1:nifx(t)=X(q)&&x(t)=X(q+1)break;endr=r+1;endxx=(x(t)-X(r))/(X(r+1)-X(r));y(t)=(xx-1)^2*(2*xx+1)*Y(r)+xx^2*(-2*xx+3)*Y(r+1)+(X(r+1)-X(r))*xx*(xx-1)^2*m(r)+(X(r+1)-X(r))*xx^2*(xx-1)*m(r+1);endend运算结果：在commondwindow中运行结果如下：clearX=[-1013];Y=[-1135];m1=6;mn=1;x=-1:0.1:3;[y,m]=spline(X,Y,x,m1,mn);plot(x,y)求得结果如下（整理之后）：x-1-0.9-0.8-0.7-0.6-0.5-0.4-0.3-0.2-0.1y-1-0.47-0.0560.2510.4720.6250.7280.7990.8560.917x00.10.20.30.40.50.60.70.80.9y11.1191.2721.4531.6561.8752.1042.3372.5682.791x11.11.21.31.41.51.61.71.81.9y33.193.3623.5173.6563.7813.8943.9964.0884.172x22.12.22.32.42.52.62.72.82.9y4.254.3234.3924.4594.5264.5944.6644.7384.8184.905x3y5其中红色数据表示插值节点将以上结果绘制成图，如下所示:习题二题目编写一通用型复化辛普生公式，能够对任意长度的等间距离散数据进行积分运算。要求：命令格式：y=simpson(x,h);y：积分结果x：输入数组h：间距程序：function[y]=simpson(x,h)%一通用型复化辛普森求积公式%y:积分结果；x:输入数组；h:间距ifnargin2disp('输入参数不够，请准确输入积分数据和步长');return;endn=length(x);%X=zeros(n+1);ifn==1y=x;return;elseifmod(n,2)==1%X=zeros(n+1);A=zeros((n+1)/2);B=zeros((n-1)/2);A=x(1:2:(end-2));B=x(2:2:(end-1));y=2*h/6*(x(end)-x(1)+4*sum(B)+2*sum(A));%X(1:n)=x;elsemod(n,2)==1X=zeros(n-1);X=x(1:(end-1));m=length(X);A=zeros((m+1)/2);B=zeros((m-1)/2);A=X(1:2:(end-2));B=X(2:2:(end-1));y=2*h/6*(X(end)-X(1)+4*sum(B)+2*sum(A))+(x(n)+x(n-1))/2*h;endendend运算结果：以计算10sin()xIdxx为例在commondwindow中运行结果如下：cleart=0:(1/8):1;x=sin(t)./t;h=1/8;[y]=simpson(x,h)y=0.9461以上积分结果与教材P65所给结果完全一致习题三题目分别用显式和隐式的二阶亚当姆斯方法求解初值问题y'=1-y，y(0)=0，令y(0.2)=0.181，取h=0.2计算y(1.0)。程序：1、显式亚当姆斯方法：function[X,Y]=Adams_shown(x0,y0,y1,h,N)%显式二阶亚当姆斯求解微分方程%x0，y0,y1为输入初值，h为步长，N为计算点数（包括除x0,y0在内的输入初值）%X,Y分别为输出数组,f为y'=f(x,y)X=zeros(N+1,1);Y=zeros(N+1,1);%为X,Y预置空间X=linspace(x0,x0+N*h,N+1);X=X(:);Y(1)=y0;Y(2)=y1;forn=2:Nyy0=f(x0,y0);yy1=f(x0+h,y1);y2=y1+h/2*(3*yy1-yy0);yy2=f(x0+2*h,y2);Y(n+1)=y2;x0=x0+h;y0=y1;y1=y2;endend2、隐式亚当姆斯方法：function[X,Y]=Adams_hidden(x0,y0,h,N)%隐式二阶亚当姆斯求解微分方程%x0，y0为输入初值，h为步长，N为计算点数（不包括输入初值）%X,Y分别为输出数组,f为y'=f(x,y)X=zeros(N+1,1);Y=zeros(N+1,1);%为X,Y预置空间X=linspace(x0,x0+N*h,N+1);X=X(:);Y(1)=y0;forn=1:Nyy0=f(x0,y0);%y1=solve('y0+h/2*(f(x0+h,y1)+yy0)-y1','y1');y1=(y0+h/2*(1+yy0))/(1+h/2);yy1=f(x0+h,y1);Y(n+1)=y1;x0=x0+h;y0=y1;endend3、二阶亚当姆斯预报—校正方法：function[X,Y]=Adams_shown_hidden(x0,y0,y1,h,N)%二阶亚当姆斯预报-校正求解微分方程%x0，y0,y1为输入初值，h为步长，N为计算点数（包括除x0,y0在内的输入初值）%X,Y分别为输出数组,f为y'=f(x,y)X=zeros(N+1,1);Y=zeros(N+1,1);%为X,Y预置空间X=linspace(x0,x0+N*h,N+1);X=X(:);Y(1)=y0;Y(2)=y1;forn=2:Nyy0=f(x0,y0);yy1=f(x0+h,y1);y2=y1+h/2*(3*yy1-yy0);yy2=f(x0+2*h,y2);y2=y1+h/2*(yy2+yy1);yy2=f(x0+2*h,y2);Y(n+1)=y2;x0=x0+h;y0=y1;y1=y2;endend运算结果：定义函数y'=1-y如下：function[yy]=f(x,y)%定义了测试函数%x，y为输入变量yy=1-y;end在commondwindow中运行结果如下：clearx0=0;y0=0;y1=0.181;h=0.2;N=5;[X1,Y1]=Adams_shown(x0,y0,y1,h,N);[X2,Y2]=Adams_hidden(x0,y0,h,N);[X3,Y3]=Adams_shown_hidden(x0,y0,y1,h,N);plot(X1,Y1,'r',X2,Y2,'b',X3,Y3,'g');legend('显式亚当姆斯','隐式亚当姆斯',,'二阶亚当姆斯校正-预报方法')最后所得结果如下表所示：X00.20.40.60.81显式00.1810.32670.44680.54540.6265隐式00.18180.33060.45230.55190.6334校正-预报00.1810.33020.45230.55210.6337根据上述数据绘制成图，如下所示：习题四题目用牛顿法求下列方程的根，要求计算结果小数点后有4位有效数字（1）30310,2xxx（2）20320,1xxxex程序：functionroot=NewtonRoot(f,x0,eps)%牛顿法求方程的根%x0为迭代初值%f为方程表达式：f(x)=0%eps为迭代精度%root为满足精度要求