数值分析大作业-三、四、五、六、七(共38页).docx

精选优质文档-倾情为你奉上大作业 三1. 给定初值及容许误差，编制牛顿法解方程f(x)=0的通用程序.解：Matlab程序如下：函数m文件：fu.mfunction Fu=fu(x)Fu=x3/3-x;end函数m文件：dfu.mfunction Fu=dfu(x)Fu=x2-1;end用Newton法求根的通用程序Newton.mclear;x0=input('请输入初值x0:');ep=input('请输入容许误差:');flag=1;while flag=1 x1=x0-fu(x0)/dfu(x0); if abs(x1-x0)<epflag=0; endx0=x1;endfprintf('方程的一个近似解为：%fn',x0); 寻找最大值的程序：Find.mcleareps=input('请输入搜索精度：');ep=input('请输入容许误差：');flag=1;k=0;x0=0;while flag=1 sigma=k*eps;x0=sigma; k=k+1; m=0; flag1=1; while flag1=1 && m<=103 x1=x0-fu(x0)/dfu(x0); if abs(x1-x0)<ep flag1=0; end m=m+1; x0=x1; end if flag1=1|abs(x0)>=ep flag=0; endendfprintf('最大的sigma值为：%fn',sigma);2求下列方程的非零根解：Matlab程序为：（1） 主程序clearclcformat longx0=765;N=100;errorlim=10(-5);x=x0-f(x0)/subs(df(),x0);n=1;while n<Nx=x0-f(x0)/subs(df(),x0);if abs(x-x0)>errorlimn=n+1;elsebreak;endx0=x;enddisp('迭代次数: n=',num2str(n)disp('所求非零根: 正根x1=',num2str(x),' 负根x2=',num2str(-x)（2）子函数 非线性函数ffunction y=f(x)y=log(513+0.6651*x)/(513-0.6651*x)-x/(1400*0.0918);end（3） 子函数 非线性函数的一阶导数dffunction y=df()syms x1y=log(513+0.6651*x1)/(513-0.6651*x1)-x1/(1400*0.0918);y=diff(y);end运行结果如下：迭代次数: n=5所求非零根: 正根x1=767.3861 负根x2=-767.3861大作业 四分析：（1）输出插值多项式。 （2）在区间-5,5内均匀插入99个节点，计算这些节点上函数f（x）的近似值，并在同一张图上画出原函数和插值多项式的图形。 （3）观察龙格现象，计算插值函数在各节点处的误差，并画出误差图。解：Matlab程序代码如下：%此函数实现y=1/(1+4*x2)的n次Newton插值，n由调用函数时指定%函数输出为插值结果的系数向量（行向量）和插值多项式function t y=func5(n)x0=linspace(-5,5,n+1)'y0=1./(1.+4.*x0.2);b=zeros(1,n+1);for i=1:n+1 s=0; for j=1:i t=1; for k=1:i if k=j t=(x0(j)-x0(k)*t; end; end; s=s+y0(j)/t; end; b(i)=s;end; t=linspace(0,0,n+1);for i=1:n s=linspace(0,0,n+1); s(n+1-i:n+1)=b(i+1).*poly(x0(1:i); t=t+s;end;t(n+1)=t(n+1)+b(1);y=poly2sym(t);10次插值运行结果： b Y=func5(10)b = Columns 1 through 4 -0.0000 0.0000 0.0027 -0.0000 Columns 5 through 8 -0.0514 -0.0000 0.3920 -0.0000 Columns 9 through 11 -1.1433 0.0000 1.0000Y = - (10035*x10)/ + x9/ + (256*x8)/93425 - x7/ - (693*x6)/1312 - (3*x5)/ + (36624*x4)/93425 - (5*x3)/ - (32311*x2)/70496 + (7*x)/ + 1b为插值多项式系数向量，Y为插值多项式。插值近似值： x1=linspace(-5,5,101); x=x1(2:100); y=polyval(b,x)y = Columns 1 through 12 2.7003 3.9994 4.3515 4.0974 3.4926 2.7237 1.9211 1.1715 0.5274 0.0154 -0.3571 -0.5960 Columns 13 through 24 -0.7159 -0.7368 -0.6810 -0.5709 -0.4278 -0.2704 -0.1147 0.0270 0.1458 0.2360 0.2949 0.3227 Columns 25 through 36 0.3217 0.2958 0.2504 0.1915 0.1255 0.0588 -0.0027 -0.0537 -0.0900 -0.1082 -0.1062 -0.0830 Columns 37 through 48 -0.0390 0.0245 0.1052 0.2000 0.3050 0.4158 0.5280 0.6369 0.7379 0.8269 0.9002 0.9549 Columns 49 through 60 0.9886 1.0000 0.9886 0.9549 0.9002 0.8269 0.7379 0.6369 0.5280 0.4158 0.3050 0.2000 Columns 61 through 72 0.1052 0.0245 -0.0390 -0.0830 -0.1062 -0.1082 -0.0900 -0.0537 -0.0027 0.0588 0.1255 0.1915 Columns 73 through 84 0.2504 0.2958 0.3217 0.3227 0.2949 0.2360 0.1458 0.0270 -0.1147 -0.2704 -0.4278 -0.5709 Columns 85 through 96 -0.6810 -0.7368 -0.7159 -0.5960 -0.3571 0.0154 0.5274 1.1715 1.9211 2.7237 3.4926 4.0974 Columns 97 through 994.3515 3.9994 2.7003绘制原函数和拟合多项式的图形代码：plot(x,1./(1+4.*x.2)hold allplot(x,y,'r')xlabel('X')ylabel('Y')title('Runge现象')gtext('原函数')gtext('十次牛顿插值多项式')绘制结果：误差计数并绘制误差图： hold off ey=1./(1+4.*x.2)-yey = Columns 1 through 12 -2.6900 -3.9887 -4.3403 -4.0857 -3.4804 -2.7109 -1.9077 -1.1575 -0.5128 -0.0000 0.3733 0.6130 Columns 13 through 24 0.7339 0.7558 0.7010 0.5921 0.4502 0.2943 0.1401 0.0000 -0.1169 -0.2051 -0.2617 -0.2870 Columns 25 through 36 -0.2832 -0.2542 -0.2053 -0.1424 -0.0719 -0.0000 0.0674 0.1254 0.1696 0.1971 0.2062 0.1962 Columns 37 through 48 0.1679 0.1234 0.0660 0.0000 -0.0691 -0.1349 -0.1902 -0.2270 -0.2379 -0.2171 -0.1649 -0.0928 Columns 49 through 60 -0.0271 0 -0.0271 -0.0928 -0.1649 -0.2171 -0.2379 -0.2270 -0.1902 -0.1349 -0.0691 0.0000 Columns 61 through 72 0.0660 0.1234 0.1679 0.1962 0.2062 0.1971 0.1696 0.1254 0.0674 0.0000 -0.0719 -0.1424 Columns 73 through 84 -0.2053 -0.2542 -0.2832 -0.2870 -0.2617 -0.2051 -0.1169 0.0000 0.1401 0.2943 0.4502 0.5921 Columns 85 through 96 0.7010 0.7558 0.7339 0.6130 0.3733 0.0000 -0.5128 -1.1575 -1.9077 -2.7109 -3.4804 -4.0857 Columns 97 through 99 -4.3403 -3.9887 -2.6900 plot(x,ey) xlabel('X') ylabel('ey') title('Runge现象误差图')输出结果为：大作业 五解：Matlab程序为：x = -520,-280,-156.6,-78,-39.62,-3.1,0,3.1,39.62,78,156.6,280,520'y = 0,-30,-36,-35,-28.44,-9.4,0,9.4,28.44,35,36,30,0'n = 13;%求解Mfor i = 1:1:n-1 h(i) = x(i+1)-x(i);endfor i = 2:1:n-1 a(i) = h(i-1)/(h(i-1)+h(i); b(i) = 1-a(i); c(i) = 6*(y(i+1)-y(i)/h(i)-(y(i)-y(i-1)/h(i-1)/(h(i-1)+h(i);enda(n) = h(n-1)/(h(1)+h(n-1);b(n) = h(1)/(h(1)+h(n-1);c(n) = 6/(h(1)+h(n-1)*(y(2)-y(1)/h(1)-(y(n)-y(n-1)/h(n-1);A(1,1) = 2;A(1,2) = b(2);A(1,n-1) = a(2);A(n-1,n-2) = a(n);A(n-1,n-1) = 2;A(n-1,1) = b(n);for i = 2:1:n-2 A(i,i) = 2; A(i,i+1) = b(i+1); A(i,i-1) = a(i+1);endC = c(2:n);C = C'm = AC;M(1) = m(n-1);M(2:n) = m;xx = -520:10.4:520;for i = 1:51 for j = 1:1:n-1 if x(j) <= xx(i) && xx(i) < x(j+1) break; end endyy(i) = M(j+1)*(xx(i)-x(j)3/(6*h(j)-M(j)*(xx(i)-x(j+1)3/(6*h(j)+(y(j+1)-M(j+1)*h(j)2/6)*(xx(i)-x(j)/h(j)-(y(j)-M(j)*h(j)2/6)*(xx(i)-x(j+1)/h(j);end;for i = 52:101 yy(i) = -yy(102-i);end;for i = 1:50 xx(i) = -xx(i);endplot(xx,yy);hold on;for i = 1:1:n/2 x(i) = -x(i);endplot(x,y,'bd');title('机翼外形曲线');输出结果：运行jywx.m文件，得到2. （1）编制求第一型3 次样条插值函数的通用程序； （2）已知汽车门曲线型值点的数据如下：解：（1）Matlab编制求第一型3 次样条插值函数的通用程序：function Sx = Threch(X,Y,dy0,dyn)%X为输入变量x的数值%Y为函数值y的数值%dy0为左端一阶导数值%dyn为右端一阶导数值%Sx为输出的函数表达式n = length(X)-1;d = zeros(n+1,1);h = zeros(1,n-1);f1 = zeros(1,n-1);f2 = zeros(1,n-2);for i = 1:n %求函数的一阶差商 h(i) = X(i+1)-X(i); f1(i) = (Y(i+1)-Y(i)/h(i);endfor i = 2:n %求函数的二阶差商 f2(i) = (f1(i)-f1(i-1)/(X(i+1)-X(i-1); d(i) = 6*f2(i);endd(1) = 6*(f1(1)-dy0)/h(1);d(n+1) = 6*(dyn-f1(n-1)/h(n-1);%赋初值A = zeros(n+1,n+1);B = zeros(1,n-1);C = zeros(1,n-1);for i = 1:n-1 B(i) = h(i)/(h(i)+h(i+1); C(i) = 1-B(i);endA(1,2) = 1;A(n+1,n) = 1;for i = 1:n+1 A(i,i) = 2;endfor i = 2:n A(i,i-1) = B(i-1); A(i,i+1) = C(i-1);endM = Ad;syms x;for i = 1:n Sx(i) = collect(Y(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i)*(x-X(i)+M(i)/2*(x-X(i)2+(M(i+1)-M(i)/(6*h(i)*(x-X(i)3); digits(4); Sx(i) = vpa(Sx(i);endfor i = 1:n disp('S(x)='); fprintf(' %s (%d,%d)n',char(Sx(i),X(i),X(i+1);endS = zeros(1,n);for i = 1:n x = X(i)+0.5; S(i) = Y(i)+(f1(i)-(M(i)/3+M(i+1)/6)*h(i)*(x-X(i)+M(i)/2*(x-X(i)2+(M(i+1)-M(i)/(6*h(i)*(x-X(i)3;enddisp('S(i+0.5)');disp(' i X(i+0.5) S(i+0.5) ');for i = 1:n fprintf(' %d %.4f %.4fn',i,X(i)+0.5,S(i);endend输出结果：>> X = 0 1 2 3 4 5 6 7 8 9 10;>> Y = 2.51 3.30 4.04 4.70 5.22 5.54 5.78 5.40 5.57 5.70 5.80;>> Threch(X,Y,0.8,0.2)S(x)= 0.8*x - 0.*x2 - 0.*x3 + 2.51 (0,1)S(x)= 0.8122*x - 0.01365*x2 - 0.*x3 + 2.506 (1,2)S(x)= 0.8218*x - 0.01849*x2 - 0.*x3 + 2.499 (2,3)S(x)= 0.317*x2 - 0.1847*x - 0.04093*x3 + 3.506 (3,4)S(x)= 6.934*x - 1.463*x2 + 0.1074*x3 - 5.986 (4,5)S(x)= 4.177*x2 - 21.26*x - 0.2686*x3 + 41.01 (5,6)S(x)= 53.86*x - 8.344*x2 + 0.427*x3 - 109.2 (6,7)S(x)= 6.282*x2 - 48.52*x - 0.2694*x3 + 129.6 (7,8)S(x)= 14.88*x - 1.643*x2 + 0.06076*x3 - 39.42 (8,9)S(x)= 8.966*x - 0.986*x2 + 0.03641*x3 - 21.67 (9,10)S(i+0.5) i X(i+0.5) S(i+0.5) 1 0.5000 2.9086 2 1.5000 3.6784 3 2.5000 4.3815 4 3.5000 4.9882 5 4.5000 5.3833 6 5.5000 5.7237 7 6.5000 5.5943 8 7.5000 5.4302 9 8.5000 5.6585 10 9.5000 5.7371 ans = - 0.*x3 - 0.*x2 + 0.8*x + 2.51, - 0.*x3 - 0.01365*x2 + 0.8122*x + 2.506, - 0.*x3 - 0.01849*x2 + 0.8218*x + 2.499, - 0.04093*x3 + 0.317*x2 - 0.1847*x + 3.506, 0.1074*x3 - 1.463*x2 + 6.934*x - 5.986, - 0.2686*x3 + 4.177*x2 - 21.26*x + 41.01, 0.427*x3 - 8.344*x2 + 53.86*x - 109.2, - 0.2694*x3 + 6.282*x2 - 48.52*x + 129.6, 0.06076*x3 - 1.643*x2 + 14.88*x - 39.42, 0.03641*x3 - 0.986*x2 + 8.966*x - 21.67大作业 六1、 炼钢厂出钢时所用的圣刚睡的钢包，在使用过程中由于钢液及炉渣对包衬耐火材料的侵蚀，使其容积不断增大，经试验，钢包的容积与相应的使用次数的数据如下：（使用次数x，容积y）xyxy2106.429110.593108.2610110.605109.5814110.726109.5016110.907109.8617110.769110.0019111.1010109.9320111.30 选用双曲线对使用最小二乘法数据进行拟合。解：Matlab程序如下：function a=nihehanshu()x0=2 3 5 6 7 9 10 11 12 14 16 17 19 20;y0=106.42 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;A=zeros(2,2);B=zeros(2,1);a=zeros(2,1); x=1./x0; y=1./y0;A(1,1)=14; A(1,2)=sum(x);A(2,1)=A(1,2);A(2,2)=sum(x.2);B(1)=sum(y); B(2)=sum(x.*y);a=AB; y=1./(a(1)+a(2)*1./x0);subplot(1,2,2);plot(x0,y0-y,'bd-'); title('拟合曲线误差');subplot(1,2,1);plot(x0,y0,'go'); hold on;x=2:0.5:20;y=1./(a(1)+a(2)*1./x);plot(x,y,'r*-'); legend('散点' ,'拟合曲线图 1/y=a(1)+a(2)*1/x');title('最小二乘法拟合曲线');试验所求的系数为：nihehanshu ans = 0.2446 0.9705则拟合曲线为 拟合曲线图、散点图、误差图如下：2、 下面给出的是乌鲁木齐最近1个月早晨7:00左右（新疆时间）的天气预报所得到的温度，按照数据找出任意次曲线拟合方程和它的图像。用MATLAB编程对上述数据进行最小二乘拟合。2008年10月-11月26日天数12345678910温度910111213141312119天数11121314151617181920温度101112131412111098天数21222324252627282930温度78911976531解：Matlab的程序如下：x=1:1:30;y=9,10,11,12,13,14,13,12,11,9,10,11,12,13,14,12,11,10,9,8,7,8,9,11,9,7,6,5,3,1;a1=polyfit(x,y,3) %三次多项式拟合%a2= polyfit(x,y,9) %九次多项式拟合%a3= polyfit(x,y,15) %十五次多项式拟合%b1= polyval(a1,x)b2= polyval(a2,x)b3= polyval(a3,x)r1= sum(y-b1).2) %三次多项式误差平方和%r2= sum(y-b2).2) %九次次多项式误差平方和%r3= sum(y-b3).2) %十五次多项式误差平方和%plot(x,y,'*') %用*画出x,y图像%hold onplot(x,b1,'r') %用红色线画出x,b1图像%hold onplot(x,b2,'g') %用绿色线画出x,b2图像%hold onplot(x,b3,'b:o') %用蓝色o线画出x,b3图像%试验结果为：a1 = Columns 1 through 2 -0.9783 -0.3645 Columns 3 through 4 0.0601 10.8912a2 = Columns 1 through 2 -0.8372 0.8217 Columns 3 through 4 -0.2354 0.4746 Columns 5 through 6 -0.7007 0.2325 Columns 7 through 8 -2.1855 11.2057 Columns 9 through 10 -19.1350 20.7351a3 = Columns 1 through 2 0.0019 -0.4224 Columns 3 through 4 0.1568 -0.3991 Columns 5 through 6 0.0612 -0.7490 Columns 7 through 8 0.8195 -0.8772 Columns 9 through 10 0.3289 -0.4112 Columns 11 through 12 3.2977 -11.2081 Columns 13 through 14 27.1590 -41.1328 Columns 15 through 16 35.7762 -4.0558b1 = Columns 1 through 2 10.6085 10.7270 Columns 3 through 4 10.3769 11.6885 Columns 5 through 6 11.7918 11.8170 Columns 7 through 8 11.8944 11.1540 Columns 9 through 10 11.7262 11.7411 Columns 11 through 12 11.3288 11.6195 Columns 13 through 14 11.7437 11.8311 Columns 15 through 16 11.0121 11.4170 Columns 17 through 18 11.1760 10.4189 Columns 19 through 20 10.2764 9.8783 Columns 21 through 22 9.3548 8.8364 Columns 23 through 24 8.4529 7.3350 Columns 25 through 26 7.6123 6.4154 Columns 27 through 28 5.8742 4.1188 Columns 29 through 30 3.2803 2.4875b2 = Columns 1 through 2 9.3759 8.0774 Columns 3 through 4 11.4371 13.8160 Columns 5 through 6 13.3358 13.5141 Columns 7 through 8 12.0927 11.5800 Columns 9 through 10 10.4281 10.2475 Columns 11 through 12 10.6524 11.7773 Columns 13 through 14 12.9173 12.2390 Columns 15 through 16 12.4590 12.0991 Columns 17 through 18 11.0301 10.2364 Columns 19 through 20 9.4028 8.9025 Columns 21 through 22 8.0895 8.5361 Columns 23 through 24 8.8635 8.9307 Columns 25 through 26 8.9951 8.2909 Columns 27 through 28 6.8981 4.7375 Columns 29 through 30 2.4474 1.0449b3 = Columns 1 through 2 8.3458 10.2630 Columns 3 through 4 10.1013 12.3552 Columns 5 through 6 13.9777 13.2023 Columns 7 through 8 13.8477 12.1538 Columns 9 through 10 10.0319 9.6735 Columns 11 through 12 9.9620 10.3149 Columns 13 through 14 12.1375 13.0373 Columns 15 through 16 13.0648 12.5279 Columns 17 through 18 11.0928 9.6272 Columns 19 through 20 8.7726 7.2293 Columns 21 through 22 7.0918 8.2312 Columns 23 through 24 9.0858 10.5636 Columns 25 through 26 9.4581 7.8568 Columns 27 through 28 5.6671 5.3263 Columns 29 through 30 2.9322 1.3906r1 = 67.7722r2 = 20.0223r3 = 3.3674其中的最后图像为：大作业 七用Romberg求积法计算积分的近似值，要求误差不超过。解：Matlab程序为：%被积函数m文件：fx.mfunction Fx=fx(x)Fx=1/(1+100*x*x);end%Romberg求积法计算积分的通用程序function Romberg()clear;a=input('请输入积分下限：a=');b=input('请输入积分上限：b=');eps=input('请输入允许精度：eps=');%=计算Tn=%function Tn=T(n)Tn=0;h=(b-a)/n;x=zeros(1,n+1);for k=1:n+1 x(k)=a+(k-1)*h;endfor j=1:n Tn=Tn+h*(fx(x(j)+fx(x(j+1)/2;endend%=计算Sn=%function Sn=S(n)Sn=4/3*T(2*n)-1/3*T(n);end%=计算Cn=%function Cn=C(n)Cn=16/15*S(2*n)-1/15*S(n);end%=计算Rn=%function Rn=R(n)Rn=64/63*C(2*n)-1/63*C(n);end%=计算满足允许精度的Rn，并打印输出=%i=1;flag=1;while flag=1if abs(R(2i)-R(2(i-1)/255<eps flag=0;endi=i+1;endfprintf('该积分的值为：%fn', R(2(i-1); End运行结