数值分析作业.docx

第二章15.方程组的性态和矩阵条件数试验。设有线性方程组ax=b，其中系数矩阵a=(aij)nxn分别为aij=xij-1(xi=1+0.1i;i,j=1,2,n)或aij=1i+j-1(i,j=1,2,n),右端向量b=（j=1na1j, j=1na2j, j=1nanj）t.利用matlab中的库函计算：（1）取n=5,10,20，分别系数矩阵的2-条件数，判别他们是否是病态矩阵？随着n的增大。矩阵形态的变化如何？（2）分别取n=5,10,20解两个线性方程组，并求的的解与精确解做比较，说明的什么？（3）取n=10，对系数矩阵中的a22和ann增加扰动10-8，求解方程组，解得变化？ 程序如下：n=input(n=);for i=1:n; for j=1:n; a(i,j)=1/(i+j-1); endendc=norm(a);d=norm(inv(a);t=c*d;b=sum(a);x=inv(a)*b;disp(t);disp(x)输入n=5时 得到t=4.7661e+005x=1.0000 1.0000 1.0000 1.0000 1.0000输入n=10时得到t=1.6025e+013x=1.0000 1.0001 1.0000 1.0000 0.9999 0.9998 0.9991 0.9999 1.0000 1.0001输入n=20时得到t= 2.3003e+018x= columns 1 through 13 0.2609 9.6667 0.9976 1.1094 0.6875 1.3750 3.3125 2.0000 16.0000 -18.0000 32.0000 -10.0000 8.0000 columns 14 through 20 -56.0000 4.0000 24.0000 39.5000 -144.0000 28.0000 -14.0000通过题意可知线性方程组的精确解为x=ones(n,1).通过对条件数t值计算可知其随n值增加而增加，n值变大矩阵变为病态矩阵，病态矩阵的结果偏差较大。（2）a22与ann发生扰动后n=input(n=);for i=1:n; for j=1:n; a(i,j)=1/(i+j-1); endenda(2,2)=a(2,2)+1*10(-8);a(n,n)=a(n,n)+1*10(-8);c=norm(a);d=norm(inv(a);t=c*d;b=sum(a);x=inv(a)*b;disp(t);disp(x)结果t= 5.2961e+011x=1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000其条件数变小，求得的解更准确。第三章14. 试分别用（1）jacobi迭代法；（2）gauss-seidel迭代法；解线形方程组（1）jacobi法 根据定义进行运算输入：a=10,1,2,3,4;1,9,-1,2,-3;2,-1,7,3,-5;3,2,3,12,-1;4,-3,-5,-1,15;b=12;-27;14;-17;12; x,k,flag=jacobi(a,b)回车得到x = 1.0000 -2.0000 3.0000 -2.0000 1.0000k = 112flag =ok其中jacobi程序为function x,k,flag = jacobi( a,b )%untitled2 summary of this function goes here% detailed explanation goes heren=length(a);k=0;x=zeros(n,1);y=zeros(n,1);flag=ok;while 1 for i=1:n y(i)=b(i); for j=1:n if j=i y(i)=y(i)-a(i,j)*x(j); end end if abs(a(i,i)1e-10|k=1000 flag=fail;return; end y(i)=y(i)/a(i,i); end if norm(y-x,inf)=1.0e-6 x0=y; y=m*x0+g;n=n+1;endfprintf(%0.6fn,y);n结果：1.000001-2.0000012.999999-2.0000000.999999n =38 第六章16. 钢包问题.炼钢厂出钢时所用的盛钢水的钢包,在使用过程中由于钢液及炉渣对 报衬耐火材料的侵蚀,使其容积不断增大。经试验，钢包的容积与相应俄使用次数的数据列表如下：使用次数x容积y使用次数x容积y23567910106.42108.26109.58109.50109.86110.00109.9311121416171920110.59110.60110.72110.90110.76111.10111.30选用双曲线对数据进行拟合，使用最小二乘法求出拟合函数，做出拟合曲线图。输入x=2,3,5,6,7,9,10,11,12,14,16,17,19,20;y=106.41,108.26,109.58,109.50,109.86,110.00,109.93,110.59,110.60,110.72,110.90,110.76,111.10,111.30;lspoly（x，y）回车得到ans = 0.0008 0.0090(即拟合方程为1y=0.009+0.008*1x)窗口继续输入：t=2:01:20;z=1./(0.009+0.0008./t);plot(x,y,ro,t,z);grid命令执行图片：其中lsploy程序如下：function c = lspoly(x,y)%untitled4 summary of this function goes here% detailed explanation goes heret=1./x;h=1./y;n=length(t);b=zeros(1:2);f=zeros(n,2);for k=(1:2) f(:,k)=t.(k-1);enda=f*f;b=f*h;c=inv(a)*b;c=flipud(c);end 第七章26.考纽螺线的形状像钟表的发条，也称回旋曲线，它在直角坐标系中的参数方程为曲线关于原点对称。取a=1，参数s的变化范围-5，5，容许误差限分别是10-3和10-7。选取适当的节点个数，利用数值积分方法计算曲线上点的坐标，并画出曲线的图形。解为：精确度为1e-3时n=201;x=zeros(n,1);y=zeros(n,1);f1=inline(cos(1/2*(t.2);f2=inline(sin(1/2*(t.2);i=1;for s= -5:0.05:5 x(i,1)=quad(f1,0,s,1e-3); y(i,1)=quad(f2,0,s,1e-3); i=i+1;endplot(x,y,m-)title(回旋曲线)xlabel(x(s)ylabel(y(s)得图精确度为1e-7时n=200;x=zeros(n,1);y=zeros(n,1);f1=inline(cos(1/2*(t.2);f2=inline(sin(1/2*(t.2);i=1;for s= -5:0.05:5 x(i,1)=quad(f1,0,s,1e-7); y(i,1)=quad(f2,0,s,1e-7); i=i+1;endplot(x,y,m-)title(回旋曲线)xlabel(x(s)ylabel(y(s)得图 第八章20.求方程在x=0.5附近的根，精确到10-8。取，用简单迭代法计算；用加快收敛的迭代格式，=0.625计算。1.(1)创建迭代函数 function p0,k,err,p = fixpt(g,p0,tol,max1)%untitled2 summary of this function goes here% detailed explanation goes herep(1)=p0;for k=2:max1 p(k)=feval(g,p(k-1); k,err=abs(p(k)-p(k-1) p=p(k); if (err fixpt(g,0.5,10(-8),100)结果为：k = 30err =1.3635e-008y = 0.5671k =31err =7.7332e-009p = columns 1 through 13 0.5000 0.6065 0.5452 0.5797 0.5601 0.5712 0.5649 0.5684 0.5664 0.5676 0.5669 0.5673 0.5671 columns 14 through 26 0.5672 0.5671 0.5672 0.5671 0.5671 0.5671 0.5671 0.5671 0.5671 0.5671 0.5671 0.5671 0.5671 columns 27 through 31 0.5671 0.5671 0.5671 0.5671 0.5671ans =0.5000结果可知循环30次，其值为0.5671。（2）fixpt函数与上题程序一致先用m文件写一个名为g.m函数fuction y=g(x);y=0.625exp(-x)+0.375*x; fixpt(g,0.5,10(-8),100)结果为：y = 0.5671k = 5err = 2.3083e-007y = 0.5671k = 6err = 4.7402e-009p = 0.5000 0.5666 0.5671 0.5671 0.5671 0.5671ans = 0.5000结果可知循环5次其值为0.5671。 第九章19. 设有常微分方程初始问题 其精确解为。选取euler和经典rk法均稳定，分别用这两种方法求解微分方程，将数值解与精确解进行比较，输出结果。rk法：f=inline(-m+2*cos(n),n,m);y=zeros(50,1);y(1)=1;h=pi/50;x=0:h:pi;for i=1:50 k1=f(x(i),y(i); k2=f(x(i)+h/2,y(i)+(h*k1)/2); k3=f(x(i)+h/2,y(i)+(h*k2)/2); k4=f(x(i)+h,y(i)+h*k3); y(i+1)=y(i)+h*(k1+2*k2+2*k3+k4)/6;end;plot(x,y,rd);hold on x1=0:pi/50:pi;y1=cos(x1)+sin(x1);plot(x