大作业 三1. 给定初值0x 及容许误差 ，编制牛顿法解方程f (x )=0的通用程序. 解：Matlab 程序如下：函数m 文件：fu.mfunction Fu=fu(x)Fu=x^3/3-x;end函数m 文件：dfu.mfunction Fu=dfu(x)Fu=x^2-1;end用Newton 法求根的通用程序Newton.mclear;x0=input('请输入初值x0:');ep=input('请输入容许误差:');flag=1;while flag==1x1=x0-fu(x0)/dfu(x0);if abs(x1-x0)<epflag=0;endx0=x1;endfprintf('方程的一个近似解为：%f\n',x0);寻找最大δ值的程序：Find.mcleareps=input('请输入搜索精度：');ep=input('请输入容许误差：');flag=1;k=0;x0=0;while flag==1sigma=k*eps;x0=sigma;k=k+1;m=0;flag1=1;while flag1==1 && m<=10^3x1=x0-fu(x0)/dfu(x0);if abs(x1-x0)<ep flag1=0;endm=m+1;x0=x1;endif flag1==1||abs(x0)>=epflag=0;endendfprintf('最大的sigma 值为：%f\n',sigma);2．求下列方程的非零根5130.6651()ln 05130.665114000.0918x x f x x +⎛⎫=-= ⎪-⨯⎝⎭解：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大作业 四试编写MATLAB 函数实现Newton 插值，要求能输出插值多项式. 对函数21()14f x x=+在区间[-5，5]上实现10次多项式插值.分析：（1）输出插值多项式。

（2）在区间[-5,5]内均匀插入99个节点，计算这些节点上函数f （x ）的近似值，并在同一张图上画出原函数和插值多项式的图形。

（3）观察龙格现象，计算插值函数在各节点处的误差，并画出误差图。

解：Matlab 程序代码如下：%此函数实现y=1/(1+4*x^2)的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+1s=0;for j=1:it=1;for k=1:iif k~=jt=(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:ns=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.0000Columns 5 through 8-0.0514 -0.0000 0.3920 -0.0000Columns 9 through 11-1.1433 0.0000 1.0000Y =b为插值多项式系数向量，Y为插值多项式。

插值近似值：x1=linspace(-5,5,101);x=x1(2:100);y=polyval(b,x)y =Columns 1 through 122.70033.99944.3515 4.0974 3.4926 2.7237 1.9211 1.1715 0.5274 0.0154 -0.3571 -0.5960Columns 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.3227Columns 25 through 360.3217 0.2958 0.2504 0.1915 0.1255 0.0588 -0.0027 -0.0537 -0.0900 -0.1082 -0.1062 -0.0830Columns 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.9549Columns 49 through 600.9886 1.0000 0.9886 0.9549 0.9002 0.8269 0.7379 0.6369 0.5280 0.4158 0.3050 0.2000Columns 61 through 720.1052 0.0245 -0.0390 -0.0830 -0.1062 -0.1082 -0.0900 -0.0537 -0.0027 0.0588 0.1255 0.1915Columns 73 through 840.2504 0.2958 0.3217 0.3227 0.2949 0.2360 0.1458 0.0270 -0.1147 -0.2704 -0.4278 -0.5709Columns 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.0974Columns 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 offey=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.6130Columns 13 through 240.7339 0.7558 0.7010 0.5921 0.4502 0.2943 0.1401 0.0000 -0.1169 -0.2051 -0.2617 -0.2870Columns 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.1962Columns 37 through 480.1679 0.1234 0.0660 0.0000 -0.0691 -0.1349 -0.1902 -0.2270 -0.2379 -0.2171 -0.1649 -0.0928Columns 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.0000Columns 61 through 720.0660 0.1234 0.1679 0.1962 0.2062 0.1971 0.1696 0.1254 0.0674 0.0000 -0.0719 -0.1424Columns 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.5921Columns 85 through 960.7010 0.7558 0.7339 0.6130 0.3733 0.0000 -0.5128 -1.1575 -1.9077 -2.7109 -3.4804 -4.0857Columns 97 through 99-4.3403 -3.9887 -2.6900plot(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-1h(i) = x(i+1)-x(i);endfor i = 2:1:n-1a(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-2A(i,i) = 2;A(i,i+1) = b(i+1);A(i,i-1) = a(i+1);endC = c(2:n);C = C';m = A\C;M(1) = m(n-1);M(2:n) = m;xx = -520:10.4:520;for i = 1:51for j = 1:1:n-1if x(j) <= xx(i) && xx(i) < x(j+1)break;endendyy(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:101yy(i) = -yy(102-i);end;for i = 1:50xx(i) = -xx(i);endplot(xx,yy);hold on;for i = 1:1:n/2x(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-1B(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+1A(i,i) = 2;endfor i = 2:nA(i,i-1) = B(i-1);A(i,i+1) = C(i-1);endM = A\d;syms x;for i = 1:nSx(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:ndisp('S(x)=');fprintf(' %s (%d,%d)\n',char(Sx(i)),X(i),X(i+1));endS = zeros(1,n);for i = 1:nx = 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:nfprintf(' %d %.4f %.4f\n',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*x^2 - 0.008514*x^3 + 2.51 (0,1)S(x)=0.8122*x - 0.01365*x^2 - 0.004458*x^3 + 2.506 (1,2)S(x)=0.8218*x - 0.01849*x^2 - 0.003652*x^3 + 2.499 (2,3)S(x)=0.317*x^2 - 0.1847*x - 0.04093*x^3 + 3.506 (3,4)S(x)=6.934*x - 1.463*x^2 + 0.1074*x^3 - 5.986 (4,5)S(x)=4.177*x^2 - 21.26*x - 0.2686*x^3 + 41.01 (5,6)S(x)=53.86*x - 8.344*x^2 + 0.427*x^3 - 109.2 (6,7)S(x)=6.282*x^2 - 48.52*x - 0.2694*x^3 + 129.6 (7,8)S(x)=14.88*x - 1.643*x^2 + 0.06076*x^3 - 39.42 (8,9) S(x)=8.966*x - 0.986*x^2 + 0.03641*x^3 - 21.67 (9,10) S(i+0.5)i X(i+0.5) S(i+0.5)1 0.5000 2.90862 1.5000 3.67843 2.5000 4.38154 3.5000 4.98825 4.5000 5.38336 5.5000 5.72377 6.5000 5.59438 7.5000 5.43029 8.5000 5.658510 9.5000 5.7371ans =[ - 0.008514*x^3 - 0.001486*x^2 + 0.8*x + 2.51, - 0.004458*x^3 - 0.01365*x^2 + 0.8122*x + 2.506, - 0.003652*x^3 - 0.01849*x^2 + 0.8218*x + 2.499, - 0.04093*x^3 + 0.317*x^2 - 0.1847*x + 3.506, 0.1074*x^3 - 1.463*x^2 + 6.934*x - 5.986, - 0.2686*x^3 + 4.177*x^2 - 21.26*x + 41.01, 0.427*x^3 - 8.344*x^2 + 53.86*x - 109.2, - 0.2694*x^3 + 6.282*x^2 - 48.52*x + 129.6, 0.06076*x^3 - 1.643*x^2 + 14.88*x - 39.42, 0.03641*x^3 - 0.986*x^2 + 8.966*x - 21.67]大作业六1、炼钢厂出钢时所用的圣刚睡的钢包，在使用过程中由于钢液及炉渣对包衬耐火材料的侵蚀，使其容积不断增大，经试验，钢包的容积与相应的使用次数的数据如下：（使用次数x，容积y）线x b a 1*y 1+=对使用最小二乘法数据进行拟合。