数值分析作业MATLAB

1.用二分法解方程x-lnx=2在区间【2，4】内的根方法:二分法算法：f=inline('x-2-log(x)');a=2;b=4;er=b-a;ya=f(a);er0=.00001;whileerer0x0=.5*(a+b);y0=f(x0);ifya*y00b=x0;elsea=x0;ya=y0;enddisp([a,b]);er=b-a;k=k+1;end求解结果：answer1343.00003.50003.00003.25003.12503.25003.12503.18753.12503.15633.14063.15633.14063.14843.14453.14843.14453.14653.14553.14653.14603.14653.14603.14623.14613.14623.14623.14623.14623.14623.14623.14623.14623.1462最终结果为:3.14622.试编写MATLAB函数实现Newton插值，要求能输出插值多项式。对函数141)(2xxf在区间[-5，5]上实现10次多项式插值。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现象误差图')3.应用牛顿迭代法于方程f(x)-1=0,导出平方根的迭代公式，用此公式计算.算法：f=inline('1-115/x^2');f1=inline('230/x^3');x0=10;er=1;k=0;whileer0.00001x=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选用双曲线xbay11对数据进行拟合，使用最小二乘法求出拟合函数，做出拟合曲线图。【解】clear,clc;%题目条件x=[2356791011121416171920];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,'*');holdon;xx=-0.5:0.01:-0.04;yy=a(1)*xx+a(2);plot(1./xx,1./yy);使用最小二乘法拟合的曲线方程为110.0089732762624460.000841690019705yx下图为绘制出的拟合曲线，并同时将一直点用“*”表示到图中。5、利用LU分解法解方程组首先，编辑一个LU分解函数如下function[L,U]=Lu(A)%求解线性方程组的三角分解法%A为方程组的系数矩阵%L和U为分解后的下三角和上三角矩阵[n,m]=size(A);ifn~=merror('TherowsandcolumnsofmatrixAmustbeequal!');return;end%判断矩阵能否LU分解forii=1:nfori=1:iiforj=1:iiAA(i,j)=A(i,j);endendif(det(AA)==0)error('ThematrixcannotbedividedbyLU!')return;endend%开始计算，先赋初值L=eye(n);U=zeros(n,n);%计算U的第一行，L的第一列fori=1:nU(1,i)=A(1,i);L(i,1)=A(i,1)/U(1,1);end%计算U的第r行，L的第r列fori=2:nforj=i:nfork=1:i-1M(k)=L(i,k)*U(k,j);endU(i,j)=A(i,j)-sum(M);endforj=i+1:nfork=1:i-1M(k)=L(j,k)*U(k,i);endL(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);ifn~=merror('TherowsandcolumnsofmatrixAmustbeequal!');return;end%判断矩阵能否LU分解forii=1:nfori=1:iiforj=1:iiAA(i,j)=A(i,j);endendif(det(AA)==0)error('ThematrixcannotbedividedbyLU!')return;endend[L,U]=Lu(A);%直接调用自定义函数，首先将矩阵分解，A=LU%设Ly=d由于L是下三角矩阵,所以可求y(i)y(1)=d(1);fori=2:nforj=1:i-1d(i)=d