当前位置:文档之家› 数值分析大作业-三、四、五、六、七

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

大作业 三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)<epflag1=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 =- (10035*x^10)/2928 + x^9/616 + (256*x^8)/93425 - x^7/76 - (693*x^6)/1312 - (3*x^5)/ + (36624*x^4)/93425 - (5*x^3)/ - (32311*x^2)/70496 + (7*x)/ + 1b为插值多项式系数向量,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现象误差图')输出结果为:大作业五1.已知直升机旋转机翼外形轮廓线上某些型值点的数据:由表中13个节点,用样条插值的方法,增加平面点为101个点,绘制出较平滑的机翼外形曲线图。

相关主题