数值分析大作业三四五六七.docx
- 文档编号:11113520
- 上传时间:2023-02-25
- 格式:DOCX
- 页数:34
- 大小:247.31KB
数值分析大作业三四五六七.docx
《数值分析大作业三四五六七.docx》由会员分享,可在线阅读,更多相关《数值分析大作业三四五六七.docx(34页珍藏版)》请在冰豆网上搜索。
数值分析大作业三四五六七
大作业三
1.给定初值
及容许误差
,编制牛顿法解方程f(x)=0的通用程序.
解:
Matlab程序如下:
函数m文件:
fu.m
functionFu=fu(x)
Fu=x^3/3-x;
end
函数m文件:
dfu.m
functionFu=dfu(x)
Fu=x^2-1;
end
用Newton法求根的通用程序Newton.m
clear;
x0=input('请输入初值x0:
');
ep=input('请输入容许误差:
');
flag=1;
whileflag==1
x1=x0-fu(x0)/dfu(x0);
ifabs(x1-x0) flag=0; end x0=x1; end fprintf('方程的一个近似解为: %f\n',x0); 寻找最大δ值的程序: Find.m clear eps=input('请输入搜索精度: '); ep=input('请输入容许误差: '); flag=1; k=0; x0=0; whileflag==1 sigma=k*eps; x0=sigma; k=k+1; m=0; flag1=1; whileflag1==1&&m<=10^3 x1=x0-fu(x0)/dfu(x0); ifabs(x1-x0) flag1=0; end m=m+1; x0=x1; end ifflag1==1||abs(x0)>=ep flag=0; end end fprintf('最大的sigma值为: %f\n',sigma); 2.求下列方程的非零根 解: Matlab程序为: (1)主程序 clear clc formatlong x0=765; N=100; errorlim=10^(-5); x=x0-f(x0)/subs(df(),x0); n=1; whilen x=x0-f(x0)/subs(df(),x0); ifabs(x-x0)>errorlim n=n+1; else break; end x0=x; end disp(['迭代次数: n=',num2str(n)]) disp(['所求非零根: 正根x1=',num2str(x),'负根x2=',num2str(-x)]) (2)子函数非线性函数f functiony=f(x) y=log((513+0.6651*x)/(513-0.6651*x))-x/(1400*0.0918); end (3)子函数非线性函数的一阶导数df functiony=df() symsx1 y=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*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+1 s=0; forj=1: i t=1; fork=1: i ifk~=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); fori=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次插值运行结果: [bY]=func5(10) b= Columns1through4 -0.00000.00000.0027-0.0000 Columns5through8 -0.0514-0.00000.3920-0.0000 Columns9through11 -1.14330.00001.0000 Y= -(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)/+1 b为插值多项式系数向量,Y为插值多项式。 插值近似值: x1=linspace(-5,5,101); x=x1(2: 100); y=polyval(b,x) y= Columns1through12 2.70033.99944.35154.09743.49262.72371.92111.17150.52740.0154-0.3571-0.5960 Columns13through24 -0.7159-0.7368-0.6810-0.5709-0.4278-0.2704-0.11470.02700.14580.23600.29490.3227 Columns25through36 0.32170.29580.25040.19150.12550.0588-0.0027-0.0537-0.0900-0.1082-0.1062-0.0830 Columns37through48 -0.03900.02450.10520.20000.30500.41580.52800.63690.73790.82690.90020.9549 Columns49through60 0.98861.00000.98860.95490.90020.82690.73790.63690.52800.41580.30500.2000 Columns61through72 0.10520.0245-0.0390-0.0830-0.1062-0.1082-0.0900-0.0537-0.00270.05880.12550.1915 Columns73through84 0.25040.29580.32170.32270.29490.23600.14580.0270-0.1147-0.2704-0.4278-0.5709 Columns85through96 -0.6810-0.7368-0.7159-0.5960-0.35710.01540.52741.17151.92112.72373.49264.0974 Columns97through99 4.35153.99942.7003 绘制原函数和拟合多项式的图形代码: plot(x,1./(1+4.*x.^2)) holdall plot(x,y,'r') xlabel('X') ylabel('Y') title('Runge现象') gtext('原函数') gtext('十次牛顿插值多项式') 绘制结果: 误差计数并绘制误差图: holdoff ey=1./(1+4.*x.^2)-y ey= Columns1through12 -2.6900-3.9887-4.3403-4.0857-3.4804-2.7109-1.9077-1.1575-0.5128-0.00000.37330.6130 Columns13through24 0.73390.75580.70100.59210.45020.29430.14010.0000-0.1169-0.2051-0.2617-0.2870 Columns25through36 -0.2832-0.2542-0.2053-0.1424-0.0719-0.00000.06740.12540.16960.19710.20620.1962 Columns37through48 0.16790.12340.06600.0000-0.0691-0.1349-0.1902-0.2270-0.2379-0.2171-0.1649-0.0928 Columns49through60 -0.02710-0.0271-0.0928-0.1649-0.2171-0.2379-0.2270-0.1902-0.1349-0.06910.0000 Columns61through72 0.06600.12340.16790.19620.20620.19710.16960.12540.06740.0000-0.0719-0.1424 Columns73through84 -0.2053-0.2542-0.2832-0.2870-0.2617-0.2051-0.11690.00000.14010.29430.45020.5921 Columns85through96 0.70100.75580.73390.61300.37330.0000-0.5128-1.1575-1.9077-2.7109-3.4804-4.0857 Columns97through99 -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; %求解M fori=1: 1: n-1 h(i)=x(i+1)-x(i); end fori=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)); end a(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-2 A(i,i)=2; A(i,i+1)=b(i+1); A(i,i-1)=a(i+1); end C=c(2: n); C=C'; m=A\C; M (1)=m(n-1); M(2: n)=m; xx=-520: 10.4: 520; fori=1: 51 forj=1: 1: n-1 ifx(j)<=xx(i)&&xx(i) break; end end yy(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; fori=52: 101 yy(i)=-yy(102-i); end; fori=1: 50 xx(i)=-xx(i); end plot(xx,yy); holdon; fori=1: 1: n/2 x(i)=-x(i); end plot(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); fori=1: n%求函数的一阶差商 h(i)=X(i+1)-X(i); f1(i)=(Y(i+1)-Y(i))/h(i); end fori=2: n%求函数的二阶差商 f2(i)=(f1(i)-f1(i-1))/(X(i+1)-X(i-1)); d(i)=6*f2(i); end d (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); fori=1: n-1 B(i)=h(i)/(h(i)+h(i+1)); C(i)=1-B(i); end A(1,2)=1; A(n+1,n)=1; fori=1: n+1 A(i,i)=2; end fori=2: n A(i,i-1)=B(i-1); A(i,i+1)=C(i-1); end M=A\d; symsx; fori=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)); end fori=1: n disp('S(x)='); fprintf('%s(%d,%d)\n',char(Sx(i)),X(i),X(i+1)); end S=zeros(1,n); fori=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; end disp('S(i+0.5)'); disp('iX(i+0.5)S(i+0.5)'); fori=1: n fprintf('%d%.4f%.4f\n',i,X(i)+0.5,S(i)); end end 输出结果: >>X=[012345678910]; >>Y=[2.513.304.044.705.225.545.785.405.575.705.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) iX(i+0.5)S(i+0.5) 10.50002.9086 21.50003.6784 32.50004.3815 43.50004.9882 54.50005.3833 65.50005.7237 76.50005.5943 87.50005.4302 98.50005.6585 109.50005.7371 ans= [-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 y x y 2 106.42 9 110.59 3 108.26 10 110.60 5 109.58 14 110.72 6 109.50 16 110.90 7 109.86 17 110.76 9 110.00 19 111.10 10 109.93 20 111.30 选用双曲线 对使用最小二乘法数据进行拟合。 解: Matlab程序如下: functiona=nihehanshu() x0=[2356791011121416171920]; y0=[106.42108.26109.58109.50109.86110.00109.93110.59110.60110.72110.90110.76111.10111.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=A\B; 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'); holdon; 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.2446 0.9705 则拟合曲线为 拟合曲线图、散点图、误差图如下: 2、下面给出的是乌鲁木齐最近1个月早晨7: 00左右(时间)的天气预报所得到的温度,按照数据找出任意次曲线拟合方程和它的图像。 用MATLAB编程对上述数据进行最小二乘拟合。 2008年10月--11月26日 天数 1 2 3 4 5 6 7 8 9 10 温度 9 10 11 12 13 14 13 12 11 9 天数 11 12 13 14 15 16 17 18 19 20 温度 10 11 12 13 14 12 11 10 9 8 天数 21 22 23 24 25 26 27 28 29 30 温度 7 8 9 11 9 7 6 5 3 1 解: 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)%九次多项式拟合
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值 分析 作业 三四 六七