数值分析作业MATLAB.doc

1.用二分法解方程 x-lnx=2 在区间【2 ，4】内的根 方法: 二分法 算法：f=inline(x-2-log(x);a=2;b=4;er=b-a; ya=f(a);er0=.00001;while erer0 x0=.5*(a+b); y0=f(x0); if ya*y0 answer1 3 4 3.0000 3.5000 3.0000 3.2500 3.1250 3.2500 3.1250 3.1875 3.1250 3.1563 3.1406 3.1563 3.1406 3.1484 3.1445 3.1484 3.1445 3.1465 3.1455 3.1465 3.1460 3.1465 3.1460 3.1462 3.1461 3.1462 3.1462 3.1462 3.1462 3.1462 3.1462 3.1462 3.1462 3.1462最终结果为: 3.14622.试编写MATLAB函数实现Newton插值，要求能输出插值多项式。对函数 在区间-5，5上实现10次多项式插值。Matlab程序代码如下：%此函数实现y=1/(1+4*x2)的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=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);for i=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次插值运行结果： b Y=func5(10)b = Columns 1 through 4 -0.0000 0.0000 0.0027 -0.0000 Columns 5 through 8 -0.0514 -0.0000 0.3920 -0.0000 Columns 9 through 11 -1.1433 0.0000 1.0000Y = - (7319042784910035*x10)/147573952589676412928 + x9/18446744073709551616 + (256*x8)/93425 - x7/1152921504606846976 - (28947735013693*x6)/562949953421312 - (3*x5)/72057594037927936 + (36624*x4)/93425 - (5*x3)/36028797018963968 - (5148893614132311*x2)/4503599627370496 + (7*x)/36028797018963968 + 1b为插值多项式系数向量，Y为插值多项式。插值近似值： x1=linspace(-5,5,101); x=x1(2:100); y=polyval(b,x)y = Columns 1 through 12 2.7003 3.9994 4.3515 4.0974 3.4926 2.7237 1.9211 1.1715 0.5274 0.0154 -0.3571 -0.5960 Columns 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.3227 Columns 25 through 36 0.3217 0.2958 0.2504 0.1915 0.1255 0.0588 -0.0027 -0.0537 -0.0900 -0.1082 -0.1062 -0.0830 Columns 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.9549 Columns 49 through 60 0.9886 1.0000 0.9886 0.9549 0.9002 0.8269 0.7379 0.6369 0.5280 0.4158 0.3050 0.2000 Columns 61 through 72 0.1052 0.0245 -0.0390 -0.0830 -0.1062 -0.1082 -0.0900 -0.0537 -0.0027 0.0588 0.1255 0.1915 Columns 73 through 84 0.2504 0.2958 0.3217 0.3227 0.2949 0.2360 0.1458 0.0270 -0.1147 -0.2704 -0.4278 -0.5709 Columns 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.0974 Columns 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 off ey=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.6130 Columns 13 through 24 0.7339 0.7558 0.7010 0.5921 0.4502 0.2943 0.1401 0.0000 -0.1169 -0.2051 -0.2617 -0.2870 Columns 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.1962 Columns 37 through 48 0.1679 0.1234 0.0660 0.0000 -0.0691 -0.1349 -0.1902 -0.2270 -0.2379 -0.2171 -0.1649 -0.0928 Columns 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.0000 Columns 61 through 72 0.0660 0.1234 0.1679 0.1962 0.2062 0.1971 0.1696 0.1254 0.0674 0.0000 -0.0719 -0.1424 Columns 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.5921 Columns 85 through 96 0.7010 0.7558 0.7339 0.6130 0.3733 0.0000 -0.5128 -1.1575 -1.9077 -2.7109 -3.4804 -4.0857 Columns 97 through 99 -4.3403 -3.9887 -2.6900 plot(x,ey) xlabel(X) ylabel(ey) title(Runge现象误差图)3.应用牛顿迭代法于方程f(x)-1=0,导出平方根的迭代公式，用此公式计算. 算法：f=inline(1-115/x2);f1=inline(230/x3);x0=10;er=1;k=0;while er0.00001 x=x0-f(x0)/f1(x0); er=abs(x-x0) x0=x; disp(x);end求解结果： answer9er =0.652210.6522er =0.070910.7231er =7.1604e-0410.7238er =7.1729e-0810.7238最终结果：10.72384.实验数据使用次数x容积y使用次数x容积y2106.4211110.593108.2612110.605109.5814110.726109.5016110.907109.8617110.769110.0019111.1010109.9320111.30选用双曲线对数据进行拟合，使用最小二乘法求出拟合函数，做出拟合曲线图。【解】clear,clc;%题目条件x=2 3 5 6 7 9 10 11 12 14 16 17 19 20;y=106.42,108.26,109.58,109.50,109.86,110.00,109.93. 110.59,110.60,110.72,110.90,110.76,111.10,111.30;%使用最小二乘法求出1次多项式拟合系数a=polyfit(1./x,1./y,1);%绘制拟合图像xx=0.04:0.01:0.5;yy=a(1)*xx + a(2);plot(1./xx,1./yy,x,y,*);hold on;xx=-0.5:0.01:-0.04;yy=a(1)*xx + a(2);plot(1./xx,1./yy);使用最小二乘法拟合的曲线方程为下图为绘制出的拟合曲线，并同时将一直点用“*”表示到图中。5、利用LU分解法解方程组首先，编辑一个LU分解函数如下functionL,U=Lu(A)% 求解线性方程组的三角分解法% A为方程组的系数矩阵%L和U为分解后的下三角和上三角矩阵n,m=size(A);if n=m error(The rows and columns of matrix A must be equal!); return;end%判断矩阵能否LU分解for ii=1:n for i=1:ii for j=1:ii AA(i,j)=A(i,j); end endif (det(AA)=0) error(The matrix can not be divided by LU!) return;endend%开始计算，先赋初值L=eye(n);U=zeros(n,n);%计算U的第一行，L的第一列for i=1:n U(1,i)=A(1,i); L(i,1)=A(i,1)/U(1,1);end%计算U的第r行，L的第r列for i=2:n for j=i:n for k=1:i-1 M(k)=L(i,k)*U(k,j); end U(i,j)=A(i,j)-sum(M); end for j=i+1:n for k=1:i-1 M(k)=L(j,k)*U(k,i); end L(j,i)=(A(j,i)-sum(M)/U(i,i); endend然后，编辑一个通过LU分解法解线性方程组的函数如下function L,U,x=Lu_x(A,d)%三角分解法求解线性方程组，LU法解线性方程组Ax=LUx=d%A为方程组的系数矩阵%d为方程组的右端项%L和U为分解后的下三角和上三角矩阵%x为线性方程组的解n,m=size(A);if n=m error(The rows and columns of matrix A must be equal!); return;end%判断矩阵能否LU分解for ii=1:n for i=1:ii for j=1:ii AA(i,j)=A(i,j); end endif (det(AA)=0) error(The matrix can not be divided by LU!) return;endendL,U=Lu(A); %直接调用自定义函数，首先将矩阵分解，A=LU%设Ly=d由于L是下三角矩阵,所以可求y(i)y(1)=d(1);for i=2:n for j=1:i-1 d(i)=d(i)-L(i,j)*y(j); end y(i)=d(i);end%设Ux=y,由于U是上三角矩阵，所以可求x(i)x(n)=y(n)/U(n,n);for i=(n-1):-1:1 for j=n:-1:i+1