数值分析各类问题题解展示及程序.docx

第二章 非线性方程求根题目：分别用二分法、牛顿迭代法、割线法、斯蒂芬森迭代法求方程的根x=1,观察不同初始值下的收敛性，并给出你的结论。解：二分法函数：function c,err,yc=bisct(f,a,b,delta)%f是所要求解函数，a、b为有根区间左右端点%delta是允许的误差界，c为所求近似解%yc为f在c处的值，err为c的误差估计if nargin0 disp( (a,b)不是有根区间 ); return,endmax1=1+round(log(b-a)-log(delta)/log(2);for k=1:max1 c=(a+b)/2; yc=feval(f,c); if yc=0 a=c; b=c; break, elseif yb*yc0 b=c; yb=yc; else a=c; ya=yc; end if(b-a) bisct(f,0.5,1.5)(a,b)不是有根区间 bisct(f,0,1)c = 1ans = 1 fun=(x)(x2+1)*(x-1)6; fplot(fun,0,2)因函数恒大于等于0，所以二分法不适用，只有在端点取值正好是根时得到结果。牛顿法函数：function p,err,k,y=newton(f,def,p0,delta,max1)%f为非线性函数，df为f微商，p0为初始值，delta为误差限，max1为迭代最大次数%p为牛顿法求得近似值，err为p0误差估计，k为迭代次数，y=f（p）p0for k=1:max1 p=p0-feval(f,p0)/feval(def,p0); err=abs(p-p0); p0=p; p,err,k,y=feval(f,p) if(err newton(f,df,0.1,10(-6),200)p = 1.0000err = 9.4200e-007k = 66y = 2.1835e-032以2为初始值： newton(f,df,2,10(-6),20)p = 1.0000err = 9.8840e-007k = 68y = 2.9137e-032割线法函数：function p,err,k,y=secant(f,p0,p,delta,max1)p0,p,feval(f,p0),feval(f,p),for k=1:max1 p2=p-feval(f,p)*(p-p0)/(feval(f,p)-feval(f,p0); err=abs(p2-p); p0=p; p=p2; p,err,k,y=feval(f,p) if(err f=inline(x2+1)*(x-1)6); secant(f,1.5,2,10(-6),100)p = 1.0000err = 9.1034e-007k = 91y = 1.9036e-031 secant(f,0,0.5,10(-6),100)p = 1.0000err = 9.7610e-007k = 88y = 2.8928e-031 secant(f,0,2,10(-6),100)p = 1.0000err = 8.9194e-007k = 101y = 1.6840e-031斯蒂芬森迭代法程序：function r,n=stevenfun(f,x0,tol)%史蒂芬森迭代法求解非线性方程的函数，f为给定迭代函数%x0为给定初始值，tol为给定精度，r为求得的根，n为迭代步数if(nargin=2) tol=1e-4;endeps=1;r=x0;n=0;while(epstol) n=n+1; r1=r; y=subs(sym(f),findsym(sym(f),r1)+r1; z=subs(sym(f),findsym(sym(f),y)+y; r=r1-(y-r1)2/(z-2*y+r1); eps=abs(r-r1);end求解：分别以0.5和1.5为初始点，进行迭代： x0=1.5; fun=(x)(x2+1)*(x-1)6; r,n=stevenfun(fun,x0)r = NaNn = 19 x0=0.5; r,n=stevenfun(fun,x0)r = NaNn =18结论：(x2+1)*(x-1)6在（0,2）上恒大于等于0（仅在1处取0），所以二分法和斯蒂芬森加速法都无法求解。牛顿法和割线法可以求解，其中牛顿法的收敛速度更快。第三章插值与拟合1. 区间-1,1做等距划分：xk= -1+kh k=0,1,n, h= n2以xk为节点对函数 fx= 15+ x2 进行插值逼近。（1）、分别取n=1,5,10,20,25，用牛顿插值对fx进行逼近，并在同一坐标系下做出函数的图形，进行比较，写出插值函数对fx的逼近程度与结点个数的关系，并分析原因；（2）、试用三次样条插值对fx进行逼近，在同一坐标下画出图形，观察样条插值函数对fx的逼近程度与结点个数的关系；（3）、整体差值有何局限性，如何避免？解：（1）牛顿插值函数：function c,d=newploy(x,y)%x为n个节点横坐标组成向量，y为对应纵坐标向量%c为所求牛顿插值多项式的系数构成的向量n=length(x);d=zeros(n,n);d(:,1)=y;for j=2:n for k=j:n d(k,j)=(d(k,j-1)-d(k-1,j-1)/(x(k)-x(k-j+1); endendc=d(n,n);for k=(n-1):-1:1 c=conv(c,poly(x(k); m=length(c); c(m)=c(m)+d(k,k);end求解及作图命令： n=1 5 10 20 25;for k=1:5h=2/n(k);x=-1:h:1;y=1./(5+x.2);newploy(x,y)plot(x,y);hold on;endtitle(不同节点牛顿插值函数图形比较);结果：不同节点插值多项式系数：（n=1） 0 0.1667（n=5） 0 0.0062 -0.0000 -0.0395 0 0.2000（n=10） -0.0000 0 0.0003 0.0000 -0.0016 0.0000 0.0080 0.0000 -0.0400 0 0.2000（n=20） 0.0000 -0.0000 -0.0000 -0.0000 0.0000 -0.0000 -0.0000 0.0000 0.0000 0.0000 -0.0001 0.0000 0.0003 0 -0.0016 0 0.0080 0 -0.0400 0 0.2000（n=25） 0 -0.0000 0 -0.0000 -0.0000 0.0000 -0.0000 -0.0000 -0.0000 0.0000 0.0000 -0.0000 -0.0000 0.0000 -0.0000 -0.0001 -0.0000 0.0003 0.0000 -0.0016 -0.0000 0.0080 0 -0.0400 -0.0000 0.2000讨论：由图形看出插值节点越多，插值函数对f（x）的逼近越好。本题中插值函数在节点增多的情况下，暂未发现龙格现象，但注意到所得插值函数的高阶项系数几乎都为0。（2）三次样条插值n=1 5 10 20 25;for k=1:5 x=-1:(1/n(k):1;y=1./(5+x.2);xi=linspace(-1,1);yi=spline(x,y,xi);plot(x,y,rp,xi,yi);hold on;end讨论：由图像看出插值节点增多时，三次样条插值对f(x)的逼近更好。整体上三次样条插值效果比牛顿插值要好。（3）整体插值时，对任意的插值节点，当n趋于无穷，高次插值函数不一定收敛于f（x），有时会出现龙格现象。所以节点多时用分段低次插值更好。2、已知一组数据如下，求其拟合曲线。i01234567891023478101114161819106.42108.2109.5110109.93110.49110.59110.6110.76111111.2（1） 求以上数据形如的拟合曲线及其平方误差;（2） 求以上数据形如的拟合曲线及其平方误差;（3） 通过观察（1）（2）的结果，写出你对数据拟合的认识。解：（1）多项式拟合function f=niheduo(x,xi)n=length(xi);for i=1:n f(i)=x(1)+x(2)*xi(i)+x(3)*xi(i)2;end输入命令： xi=2 3 4 7 8 10 11 14 16 18 19;yi=106.42 108.2 109.5 110 109.93 110.49 110.59 110.6 110.76 111 111.2 ;x0=1 1 1;x,resnorm,residual,exitflag,output,lambda,jacobian=lsqcurvefit(niheduo,x0,xi,yi)结果：x = 106.2927 0.6264 -0.0205resnorm = 2.7796residual = Columns 1 through 10 1.0436 -0.2123 -1.0292 -0.3251 0.0645 0.0208 0.1176 0.4526 0.3181 -0.0601 Column 11 -0.3905exitflag = 1output = firstorderopt: 3.6806e-012 iterations: 5 funcCount: 24 cgiterations: 0 algorithm: large-scale: trust-region reflective Newton message: 1x427 charlambda = lower: 3x1 doubleupper: 3x1 double作图： X=2:0.1:19;Y=x(1)+x(2)*X+x(3)*X.2;plot(xi,yi,rp);grid on; hold on;plot(X,Y);legend(观测数据点,拟合曲线);得：结论：拟合曲线为,平方误差为2.7796。（2）按所给函数拟合function f=nihe(x,xi)n=length(xi);for i=1:n f(i)=x(1)*exp(-x(2)/xi(i);end输入命令：x0=1 1;x,resnorm,residual,exitflag,output,lambda,jacobian=lsqcurvefit(nihe,x0,xi,yi)结果：x = 111.4922 0.0902resnorm =0.4719注：其他结果略去作图：X=2:0.1:19;Y=x(1)*exp(-x(2)./X);plot(xi,yi,rp);grid on; hold on;plot(X,Y);legend(观测数据点,拟合曲线);得：结论：拟合曲线为，平方误差为0.4719。（3）由前两问结果发现（2）拟合函数的平方误差比（1）的小了不少，说明选择合适的拟合曲线类型可以明显减小误差，先分析要拟合的函数类型是很必要的。第四章数值微分和数值积分题目：设计区间分半求积算法、龙贝格求积算法和自适应辛普森求积算法的程序，观察n=1，10，100，500，积分的结果，并给出相应的评价。解：（1）、区间分半求积算法：function x=fenban(f,a,b,delta)if nargindelta t0=t1, h=h/2; t1=t0/2; for k=1:2n t1=t1+h*feval(f,a+(2*k-1)*h); end n=n+1;endx=t1 ; %结果n %迭代次数求解： k=1 10 100 500; f=(x)(x2)*cos(k*x) ; x=fenban(f,-1,1)x =0.4783 -0.1402 -0.0092 0.0014n = 8 10 12 13(2)、龙贝格求积算法函数：function r,quad,err,h=romber(f,a,b,n,delta)%f为被积函数，a、b为积分上下限，n+1为T数表列数，delta为允许误差%r是T数表，quad是所求积分值m=1;h=b-a; err=1; j=0;r=zeros(4,4);r(1,1)=h*(feval(f,a)+feval(f,b)/2;while(errdelta)&(jn)|(jfun=(x)x2*cos(n*x) romber(fun,-1,1,7,1e-5)quad = 0.4783 -0.1402 0.6112 -0.0190迭代次数：5 8 5 6（3）、自适应辛普森求积算法函数：function j,num,tol = changesps(f,a,b,tol)%变步长simpson积分函数，a、b为积分区间，tol为积分精度if(nargin=3) tol=1e-5;endn=1; h=(b-a)/2;t1=0; x=0;t2=feval(f,a)+feval(f,b)/h;while abs(t2-t1)tol n=n+1; h=(b-a)/n; t1=t2; t2=0; for i=1:n-1 x=a+h*i; x1=x+h; t2=t2+(h/6)*(feval(f,x)+4*feval(f,(x+x1)/2)+feval(f,x1); endendj=t2; num=n;求解： f=(x)x2*cos(1*x); s,num,tol=changesps(f,-1,1,1e-5)s = 0.4750num = 329tol = 1.0000e-005对不同的n，f=(x)x2*cos(1*x)中1替换为10、100、500，得:s = 0.4750 -0.1361 -0.0146 -0.0024num = 329 415 169 119（4）、由各算法迭代次数可知，在精度误差一致的情况下，龙贝格求积算法的收敛速度最快，区间分半算法次之，自适应辛普森算法最慢。第五章常微分方程初值问题的数值解法1、分别取步长h=0.5，0.1，0.05，0.01，。，用显式欧拉法求解，y(1)=1,计算到t=2，并与精确解比较，观察欧拉法的收敛性。解：显式欧拉法函数：function E=Eule(fun,x0,y0,xN,N)%fun为一阶微分方程的函数，x0、y0为初始条件，xN为取值范围一个端点，h为区间步长%N为区间个数，x为Xn构成的向量，y为Yn构成的向量x=zeros(1,N+1); y=zeros(1,N+1);x(1)=x0;y(1)=y0;h=(xN-x0)/N;for n=1:N x(n+1)=x(n)+h; y(n+1)=y(n)+h*feval(fun,x(n),y(n);endT=y(N+1)求解： f=(x,y)x*y(1/3);h=0.5;s=Eule(f,1,1,2,1/h);xi=1:h:2;plot(xi,s,p); hold on;h=0.1;s=Eule(f,1,1,2,1/h);xi=1:h:2;plot(xi,s,*); hold on;h=0.05;s=Eule(f,1,1,2,1/h);xi=1:h:2;plot(xi,s,-); hold on;h=0.01;s=Eule(f,1,1,2,1/h);xi=1:h:2;plot(xi,s,:); hold on;h=0.005;s=Eule(f,1,1,2,1/h);xi=1:h:2;plot(xi,s,-.); hold on; X=1:0.1:2;Y=(X.2+2)/3).(3/2);plot(X,Y,rp);结果：各步长t=2处求解的值：T = 2.3585 2.7239 2.7754 2.8177 2.8231 由题目所给精确解计算t=2点的值为: y=(22+2)/3)(3/2)y = 2.8284可见步长h取0.005时，t=2处的值与精确解最为接近。再结合图像，不难发现步长越短，欧拉法求得的解与精确解越接近，此显式欧拉法收敛。2.分别取步长 h= 0.1 , 0.05 , 0.01,0.0005 ，用显示欧拉法和隐式欧拉法求解y= -5y ,y(0)=1由结果分析算法的稳定性。（注：原题所给初值过大，为更好展示，作了修改）。 解：显式欧拉法函数如下：function x,y=Euler(fun,x0,y0,xN,N)%fun为一阶微分方程的函数，x0、y0为初始条件，xN为取值范围一个端点，h为区间步长%N为区间个数，x为Xn构成的向量，y为Yn构成的向量x=zeros(1,N+1); y=zeros(1,N+1);x(1)=x0;y(1)=y0;h=(xN-x0)/N;for n=1:N x(n+1)=x(n)+h; y(n+1)=y(n)+h*feval(fun,x(n),y(n);end用显式欧拉法作出图象： f=(x,y)-5*y;h=0.1;x,y=Euler(f,0,1,1,1/h);plot(x,y,.); hold on;h=0.05;x,y=Euler(f,0,1,1,1/h);plot(x,y,*); hold on;h=0.01;x,y=Euler(f,0,1,1,1/h);plot(x,y,-); hold on;h=0.005;x,y=Euler(f,0,1,1,1/h);plot(x,y,:); hold on;隐式欧拉法函数如下：function x,y=euler2(fun,x0,y0,xN,N)%欧拉向后公式，fun为一阶微分方程函数，x0、y0为初始条件%xN为取值范围一端点，h为步长，n为区间个数x=zeros(1,N+1); y=zeros(1,N+1);x(1)=x0;y(1)=y0;h=(xN-x0)/N;for n=1:N x(n+1)=x(n)+h; z0=y(n)+h*feval(fun,x(n),y(n); for k=1:3 z1=y(n)+h*feval(fun,x(n+1),z0); if abs(z1-z0)k temp=Aug(k,:); Aug(k,:)=Aug(r,:); Aug(r,:)=temp; end if Aug(k,k)=0, error(对角元出现0), end % 把增广矩阵消元成为上三角 for p = k+1:n Aug(p,:)=Aug(p,:)-Aug(k,:)*Aug(p,k)/Aug(k,k); end end % 解上三角方程组 A = Aug(:,1:n); b = Aug(:,n+1); x(n) = b(n)/A(n,n); for k = n-1:-1:1 x(k)=b(k); for p=n:-1:k+1 x(k) = x(k)-A(k,p)*x(p); end x(k)=x(k)/A(k,k); end求解希尔伯特矩阵，假定右端项元素全为1，则： A=hilb(n);b=ones(n,1);x=gaosi(A,b)n值分别取5、10、20、50、100，得:x = 1.0e+003 * 0.0050 -0.1200 0.6300 -1.1200 0.6300x = 1.0e+006 * -0.0000 0.0010 -0.0238 0.2402 -1.2610 3.7831 -6.7257 7.0002 -3.9377 0.9237x = 1.0e+008 * -0.0000 0.0000 -0.0001 -0.0008 0.0239 -0.2175 1.0227 -2.7057 3.8790 -2.3342 -0.0956 -0.4873 0.0819 6.1175 -9.6550 5.3460 -3.1653 5.4801 -4.5446 1.2551x = 1.0e+010 * Columns 1 through 10 0.0000 -0.0000 0.0000 -0.0005 0.0048 -0.0276 0.0944 -0.1877 0.1831 -0.0126。 Columns 41 through 50 -1.7255 0.0547 0.3567 0.1077 0.8827 -0.4095 0.3212 -0.4466 -0.3768 0.3457x = 1.0e+010 * Columns 1 through 10 -0.0000 0.0000 -0.0001 0.0015 -0.0143 0.0770 -0.2408 0.4117 -0.2839 -0.0555。 Columns 91 through 100 2.0122 -1.6636 -1.4601 -0.7345 -2.1395 1.7803 -0.8700 0.5395 0.3762 0.0392带入验证： A*xans = 1.0000 1.0000 1.0000 1.0000 。 正确第八章 矩阵特征值问题的解法题目：给定矩阵（1）.用幂法计算主特征值和特征向量；（2）.选择不同初值，观察所需迭代次数和结果，说明幂法对初值选择的要求；（3）.对矩阵B=A-pI(p分别取1、2、3、4)，幂法求B的主特征值和主特征向量，相同精度下，比较收敛速度，与（1）比较；（4）.用反幂法计算该矩阵按模最小的特征值和相应特征向量。解:幂法程序：function m,u,flag=pows(A,u0,delta,max1)%A为所给矩阵，delta为精度要求，max1为最大迭代次数%m为绝对值最大特征值，u为m的特征向量，flag标识成功或失败if nargin4 max1=100; endif nargin3 delta=1e-5; endn=length(A);flag=失败; k=0; m1=0; u=u0;while k=max1 v=A*u; vmax,i=max(abs(v); m=v(i);u=v/m; if abs(m-m1)delta flag=成功;break; end m1=m; k=k+1;endk反幂法程序：f