数值分析大作业-三、四、五、六、七

大作业三1.给定初值0x及容许误差，编制牛顿法解方程f(x)=0的通用程序.解：Matlab程序如下：函数m文件：fu.mfunctionFu=fu(x)Fu=x^3/3-x;end函数m文件：dfu.mfunctionFu=dfu(x)Fu=x^2-1;end用Newton法求根的通用程序Newton.mclear;x0=input('请输入初值x0:');ep=input('请输入容许误差:');flag=1;whileflag==1x1=x0-fu(x0)/dfu(x0);ifabs(x1-x0)epflag=0;endx0=x1;endfprintf('方程的一个近似解为：%f\n',x0);寻找最大δ值的程序：Find.mcleareps=input('请输入搜索精度：');ep=input('请输入容许误差：');flag=1;k=0;x0=0;whileflag==1sigma=k*eps;x0=sigma;k=k+1;m=0;flag1=1;whileflag1==1&&m=10^3x1=x0-fu(x0)/dfu(x0);ifabs(x1-x0)epflag1=0;endm=m+1;x0=x1;endifflag1==1||abs(x0)=epflag=0;endendfprintf('最大的sigma值为：%f\n',sigma);2．求下列方程的非零根5130.6651()ln05130.665114000.0918xxfxx解：Matlab程序为：（1）主程序clearclcformatlongx0=765;N=100;errorlim=10^(-5);x=x0-f(x0)/subs(df(),x0);n=1;whilenNx=x0-f(x0)/subs(df(),x0);ifabs(x-x0)errorlimn=n+1;elsebreak;endx0=x;enddisp(['迭代次数:n=',num2str(n)])disp(['所求非零根:正根x1=',num2str(x),'负根x2=',num2str(-x)])（2）子函数非线性函数ffunctiony=f(x)y=log((513+0.6651*x)/(513-0.6651*x))-x/(1400*0.0918);end（3）子函数非线性函数的一阶导数dffunctiony=df()symsx1y=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()14fxx在区间[-5，5]上实现10次多项式插值.分析：（1）输出插值多项式。（2）在区间[-5,5]内均匀插入99个节点，计算这些节点上函数f（x）的近似值，并在同一张图上画出原函数和插值多项式的图形。（3）观察龙格现象，计算插值函数在各节点处的误差，并画出误差图。解：Matlab程序代码如下：%此函数实现y=1/(1+4*x^2)的n次Newton插值，n由调用函数时指定%函数输出为插值结果的系数向量（行向量）和插值多项式function[ty]=func5(n)x0=linspace(-5,5,n+1)';y0=1./(1.+4.*x0.^2);b=zeros(1,n+1);fori=1:n+1s=0;forj=1:it=1;fork=1:iifk~=jt=(x0(j)-x0(k))*t;end;end;s=s+y0(j)/t;end;b(i)=s;end;t=linspace(0,0,n+1);fori=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次插值运行结果：[bY]=func5(10)b=Columns1through4-0.00000.00000.0027-0.0000Columns5through8-0.0514-0.00000.3920-0.0000Columns9through11-1.14330.00001.0000Y=-(7319042784910035*x^10)/147573952589676412928+x^9/18446744073709551616+(256*x^8)/93425-x^7/1152921504606846976-(28947735013693*x^6)/562949953421312-(3*x^5)/72057594037927936+(36624*x^4)/93425-(5*x^3)/36028797018963968-(5148893614132311*x^2)/4503599627370496+(7*x)/36028797018963968+1b为插值多项式系数向量，Y为插值多项式。插值近似值：x1=linspace(-5,5,101);x=x1(2:100);y=polyval(b,x)y=Columns1through122.70033.99944.35154.09743.49262.72371.92111.17150.52740.0154-0.3571-0.5960Columns13through24-0.7159-0.7368-0.6810-0.5709-0.4278-0.2704-0.11470.02700.14580.23600.29490.3227Columns25through360.32170.29580.25040.19150.12550.0588-0.0027-0.0537-0.0900-0.1082-0.1062-0.0830Columns37through48-0.03900.02450.10520.20000.30500.41580.52800.63690.73790.82690.90020.9549Columns49through600.98861.00000.98860.95490.90020.82690.73790.63690.52800.41580.30500.2000Columns61through720.10520.0245-0.0390-0.0830-0.1062-0.1082-0.0900-0.0537-0.00270.05880.12550.1915Columns73through840.25040.29580.32170.32270.29490.23600.14580.0270-0.1147-0.2704-0.4278-0.5709Columns85through96-0.6810-0.7368-0.7159-0.5960-0.35710.01540.52741.17151.92112.72373.49264.0974Columns97through994.35153.99942.7003绘制原函数和拟合多项式的图形代码：plot(x,1./(1+4.*x.^2))holdallplot(x,y,'r')xlabel('X')ylabel('Y')title('Runge现象')gtext('原函数')gtext('十次牛顿插值多项式')绘制结果：误差计数并绘制误差图：holdoffey=1./(1+4.*x.^2)-yey=Columns1through12-2.6900-3.9887-4.3403-4.0857-3.4804-2.7109-1.9077-1.1575-0.5128-0.00000.37330.6130Columns13through240.73390.75580.70100.59210.45020.29430.14010.0000-0.1169-0.2051-0.2617-0.2870Columns25through36-0.2832-0.2542-0.2053-0.1424-0.0719-0.00000.06740.12540.16960.19710.20620.1962Columns37through480.16790.12340.06600.0000-0.0691-0.1349-0.1902-0.2270-0.2379-0.2171-0.1649-0.0928Columns49through60-0.02710-0.0271-0.0928-0.1649-0.2171-0.2379-0.2270-0.1902-0.1349-0.06910.0000Columns61through720.06600.12340.16790.19620.20620.19710.16960.12540.06740.0000-0.0719-0.1424Columns73through84-0.2053-0.2542-0.2832-0.2870-0.2617-0.2051-0.11690.00000.14010.29430.45020.5921Columns85through960.70100.75580.73390.61300.37330.0000-0.5128-1.1575-1.9077-2.7109-3.4804-4.0857Columns97through99-4.3403-3.9887-2.6900plot(x,ey)xlabel('X')ylabel('ey')title('Runge现象误差图')输出结果为：大作业五1.已知直升机旋转机翼外形轮廓线上某些型值点的数据：kkxkykkxky152002280-303156.6-36478-35539.62-28.4463.1-9.470083.19.4939.6228.4410783511156.6361228030135200由表中13个节点，用样条插值的方法，增加平面点为101个点，绘制出较平滑的机翼外形曲线图。解：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;%求解Mfori=1:1:n-1h(i)=x(i+1)-x(i);endfori=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);fori=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;fori=1:51forj=1:1:n-1ifx(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)