数值分析大作业三四五六七HEN system office room 【HEN16H-HENS2AHENS8Q8-HENH1688】大作业 三1.给定初值0x 及容许误差，编制牛顿法解方程f (x )=0的通用程序.解：Matlab 程序如下：函数m 文件：function Fu=fu(x) Fu=x^3/3-x; end函数m 文件：function Fu=dfu(x) Fu=x^2-1; end用Newton 法求根的通用程序 clear;x0=input('请输入初值x0:'); ep=input('请输入容许误差:'); flag=1;while flag==1x1=x0-fu(x0)/dfu(x0); if abs(x1-x0)<ep flag=0; end x0=x1; endfprintf('方程的一个近似解为：%f\n',x0); 寻找最大δ值的程序： cleareps=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<=10^3 x1=x0-fu(x0)/dfu(x0);if abs(x1-x0)<epend m=m+1; x0=x1; endif flag1==1||abs(x0)>=ep flag=0; end endfprintf('最大的sigma 值为：%f\n',sigma);2．求下列方程的非零根5130.6651()ln 05130.665114000.0918x x f x x +⎛⎫=-= ⎪-⨯⎝⎭解：Matlab 程序为：（1）主程序 clear clc format long x0=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)>errorlim n=n+1; else break ; end x0=x; enddisp(['迭代次数: n=',num2str(n)])disp(['所求非零根: 正根x1=',num2str(x),' 负根x2=',num2str(-x)]) （2）子函数 非线性函数f function y=f(x)y=log((513+*x)/*x))-x/(1400*; end（3）子函数 非线性函数的一阶导数df function y=df() syms x1y=log((513+*x1)/*x1))-x1/(1400*; y=diff(y);运行结果如下：迭代次数: n=5所求非零根: 正根x1= 负根x2=大作业 四试编写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+1 s=0; for j=1:i t=1; for k=1:i if 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)Columns 1 through 4Columns 5 through 8Columns 9 through 11Y =b为插值多项式系数向量，Y为插值多项式。

插值近似值：x1=linspace(-5,5,101);x=x1(2:100);y=polyval(b,x)y =Columns 1 through 12Columns 13 through 24Columns 25 through 36Columns 37 through 48Columns 49 through 60Columns 61 through 72Columns 73 through 84Columns 85 through 96Columns 97 through 99绘制原函数和拟合多项式的图形代码：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 13 through 24Columns 25 through 36Columns 37 through 48Columns 49 through 60Columns 61 through 72Columns 73 through 84Columns 85 through 96Columns 97 through 99plot(x,ey)xlabel('X')ylabel('ey')title('Runge现象误差图')输出结果为：大作业五解：Matlab程序为：x = [-520,-280,,-78,,,0,,,78,,280,520]';y = [0,-30,-36,-35,,,0,,,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;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::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('机翼外形曲线');输出结果：运行文件，得到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;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)+;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+');disp(' i X(i+ S(i+ ');for i = 1:nfprintf(' %d %.4f %.4f\n',i,X(i)+,S(i));endend输出结果：>> X = [0 1 2 3 4 5 6 7 8 9 10];>> Y = [ ];>> Threch(X,Y,,S(x)=*x - *x^2 - *x^3 + (0,1)S(x)=*x - *x^2 - *x^3 + (1,2)S(x)=*x - *x^2 - *x^3 + (2,3)S(x)=*x^2 - *x - *x^3 + (3,4)S(x)=*x - *x^2 + *x^3 - (4,5)S(x)=*x^2 - *x - *x^3 + (5,6)S(x)=*x - *x^2 + *x^3 - (6,7)S(x)=*x^2 - *x - *x^3 + (7,8)S(x)=*x - *x^2 + *x^3 - (8,9)S(x)=*x - *x^2 + *x^3 - (9,10)S(i+i X(i+ S(i+12345678910[ - *x^3 - *x^2 + *x + , - *x^3 - *x^2 + *x + , - *x^3 - *x^2 + *x + , - *x^3 + *x^2 - *x + , *x^3 - *x^2 + *x - , - *x^3 + *x^2 - *x + , *x^3 - *x^2 + *x - , - *x^3 + *x^2 - *x + , *x^3 - *x^2 + *x - , *x^3 - *x^2 + *x - ]大作业六1、炼钢厂出钢时所用的圣刚睡的钢包，在使用过程中由于钢液及炉渣对包衬耐火材料的侵蚀，使其容积不断增大，经试验，钢包的容积与相应的使用次数的数据如下：（使用次数x，容积y）选用线1=ya对使用最小二乘法数据进行拟合。