数值分析大作业 三、四、五、六、七(35页).doc

-大作业 三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 = - (7319042784910035*x10)/147573952589676412928 + x9/18446744073709551616 + (256*x8)/93425 - x7/1152921504606846976 - (28947735013693*x6)/562949953421312 - (3*x5)/72057594037927936 + (36624*x4)/93425 - (5*x3)/36028797018963968 - (5148893614132311*x2)/4503599627370496 + (7*x)/36028797018963968 + 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.001486*x2 - 0.008514*x3 + 2.51 (0,1)S(x)= 0.8122*x - 0.01365*x2 - 0.004458*x3 + 2.506 (1,2)S(x)= 0.8218*x - 0.01849*x2 - 0.003652*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.008514*x3 - 0.001486*x2 + 0.8*x + 2.51, - 0.004458*x3 - 0.01365*x2 + 0.8122*x + 2.506, - 0.003652*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.008973276262446 0.000841690019705则拟合曲线为 拟合曲线图、散点图、误差图如下：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.000291429309783 -0.009415524223645 Columns 3 through 4 0.292004142810601 10.208538587848912a2 = Columns 1 through 2 -0.000000003138372 0.000000520938217 Columns 3 through 4 -0.000035908922354 0.001336820044746 Columns 5 through 6 -0.029221194837007 0.381184877642325 Columns 7 through 8 -2.877296607421855 11.503362068342057 Columns 9 through 10 -19.950300302531350 20.434416122867351a3 = Columns 1 through 2 0.000000000000019 -0.000000000004224 Columns 3 through 4 0.000000000421568 -0.000000025043991 Columns 5 through 6 0.000000985890612 -0.000027096197490 Columns 7 through 8 0.000533732138195 -0.007616987378772 Columns 9 through 10 0.078769002923289 -0.585435616904112 Columns 11 through 12 3.079422065442977 -11.201853510022081 Columns 13 through 14 27.180413929381590 -41.322246162561328 Columns 15 through 16 35.849674463847762 -4.072038248460558b1 = Columns 1 through 2 10.490835777126085 10.752553342097270 Columns 3 through 4 10.991942706903769 11.207255295686885 Columns 5 through 6 11.396742532587918 11.558655841748170 Columns 7 through 8 11.691246647308944 11.792766373411540 Columns 9 through 10 11.861466444197262 11.895598283807411 Columns 11 through 12 11.893413316383288 11.853162966066195 Columns 13 through 14 11.773098656997437 11.651471813318311 Columns 15 through 16 11.486533859170121 11.276536218694170 Columns 17 through 18 11.019730316031760 10.714367575324189 Columns 19 through 20 10.358699420712764 9.950977276338783 Columns 21 through 22 9.489452566343548 8.972376714868364 Columns 23 through 24 8.398001146054529 7.764577284043350 Columns 25 through 26 7.070356552976123 6.313590376994154 Columns 27 through 28 5.492530180238742 4.605427386851188 Columns 29 through 30 3.650533420972803 2.626099706744875b2 = Columns 1 through 2 9.463446392983759 8.773862633430774 Columns 3 through 4 11.101355276104371 13.121450915718160 Columns 5 through 6 13.808995551263358 13.310996466695141 Columns 7 through 8 12.218048948080927 11.136022854345800 Columns 9 through 10 10.474994026394281 10.386576634782475 Columns 11 through 12 10.793844829246524 11.469925464287773 Columns 13 through 14 12.132098233639173 12.527855252812390 Columns 15 through 16 12.497848983074590 12.007995392020991 Columns 17 through 18 11.149198395600301 10.108221923752364 Columns 19 through 20 9.117157396954028 8.391716993829025 Columns 21 through 22 8.070226830600895 8.165699061425361 Columns 23 through 24 8.542727944878635 8.929182105649307 Columns 25 through 26 8.969753551969951 8.324373488882909 Columns 27 through 28 6.809315594398981 4.572479200277375 Columns 29 through 30 2.286877742084474 1.337751913240449b3 = Columns 1 through 2 8.999596533473458 10.005535778432630 Columns 3 through 4 10.975000499091013 12.035470425513552 Columns 5 through 6 13.067199404249777 13.693579592992023 Columns 7 through 8 13.363208660348477 12.052126723871538