数值分析实验报告

石家庄经济学院2007/2008学年第一学期数值分析实验报告班级三班学号405417080321姓名连欣指导教师马丽１实验一牛顿插值一、实验题目已知函数表：Xi0.400.550.560.80Yi0.410.570.690.88用三次牛顿插值多项式求x=0.59时的函数值二、程序功能根据给定的数据构造三次牛顿插值多项式，并且计算出x=0.59时的函数值三、算法步1：计算dx=x-Xi,即(x-x1),(x-x2),…(x-xn)。步2：构造插商表即求D(i,j)。步3：构造一个三次牛顿插值多项式即Pn=sum(yi)。步4：计算结果。四、重要标识符说明X:函数中的自变量；Xi:插值接点；Y:在函数作用下的X的值；dx：=x-Xi；formatlong是定义，使小数点位数增多，结果更精确；prod()是连乘函数；五、程序运行实例在matlab6.5环境中，运行程序，结果如图：六、源程序%NewtonInt.m%clc;clear;２formatlong;Xi=[0.400.550.560.80];Y=[0.410.570.690.88];%Fun=In(X)x=0.59;%theinterpolativevalue;%(x-x1),(x-x2),¡­,(x-xn)n=length(Xi);yi(1)=Y(1);D(:,1)=Y(:);forj=2:nfori=j:nD(i,j)=(D(i,j-1)-D(i-1,j-1))/(Xi(i)-Xi(i-j+1));%constructadivided-differencetableendyi(j)=D(j,j)*prod(dx(1:j-1));%f(x1x2¡­xn)%(x-x1)(x-x2)¡­(x-xn)endPn=sum(yi);%Pn=f1+f(x1x2)*(x-x1)+¡­disp(Pn);七、个人实验总结第一次使用matlab环境，才初步学会了它的一些基本运用功能，比如.m文件的新建和保存，命令窗口的清空，历史窗口的调出与隐藏等。对牛顿法的计算实质也有了进一步的了解。３实验二曲线拟合的最小二乘法一、实验题目已知函数表：Xi2468Yi1.12.84.97.2使用最小二乘法计算拟合多项式y=a0+a1*x二、程序功能根据给定的数据拟合出多项式y=a0+a1x三、算法步1:根据多项式列出关于a0,a1,x,y的向量形式；步2：等式两边左乘A’；步3：利用向量算法求出a0,a1.四、重要标识符说明n:x和y的长度;A’:A的转置；reshape(a,m,n):把矩阵a变成n×n阶矩阵；polyfit(x,y,n)其中x,y为拟合数据，n为拟合多项式的阶数；ones(m,n):生成m×n阶全1矩阵；y=linspace(x1,x2,n),生成n维向量，使得y(1)=x1,y(n)=x2。五、程序运行实例在matlab6.5环境中，运行程序，结果如图：４六、源程序N=1clc;clear;x=[2468];y=[1.12.84.97.2];n=1;p=polyfit(x,y,1);xi=linspace(0,10,100);yi=polyval(p,xi);plot(x,y,'o',xi,yi);xlabel('xmm');ylabel('ymm');legend;pause;５N=4function[u]=leastfit(x,y)x=[2468];y=[1.12.84.97.2];n=4;x=x(1:n);y=y(1:n);x=reshape(x,n,1);y=reshape(y,n,1);A=[ones(n,1)x];b=y;B=A'*A;b=A'*b;u=B\b;七、个人实验总结有了第一次实验对MATLAB环境的熟悉基础，本次试验应用起来似乎更方便了些。在n不同时，利用最小二乘法计算同一多项式会得出不同的显示结果。实验三复合梯形求积一、实验题目对于函数f(x)=sinx/x,利用下表计算此函数在0-1区间的积分X01/81/43/81/25/83/47/81f(x)10.99739780.98961580.97672670.95885100.93615560.90885160.87719250.8414709６二、程序功能根据复合梯形求积公式，计算函数的积分三、算法步1：计算步长h。即h=(max(x)-min(x))/(n-1)；步2：计算各个结点对应的函数值。即f(1:n)=y(1:n)；步3：将第二步所得结果代入复合梯形公式，即I=(f(1)+2*sum(f(2:n-1))+f(n))*h/2；步4：计算结果。四、重要标识符说明x:函数中的自变量；xi:每个结点处的自变量值；y:每个自变量对应的函数值；h:步长；max()是取最大值；min()是取最小值；sum是求和函数；五、程序运行实例在matlab6.5环境中，运行程序，结果如图：六、源程序%%%%%%%%%%%%%%%%%%%%%%%%%%%ComTrapz.m%%%%%%%%%%%%%%%%%%%%%%%%%%clear;n=9;%x=[01/81/43/81/25/83/47/81];y=[10.99739780.98961580.97672670.95885100.93615560.9088516７0.87719250.8414709];h=(max(x)-min(x))/(n-1);f(1:n)=y(1:n);I=(f(1)+2*sum(f(2:n-1))+f(n))*h/2;%¸disp(I);七、个人实验总结在matlab中用了很多矩阵和向量来解决问题，更加方便，不用再定义这些矩阵，直接输入就可以。掌握节点的个数和公式，了解一些内函数的相关知识，即可理解整个程序。实验四二分法一、实验题目求方程f(x)=x3-x-1=0在区间[1.0，1.5]内的一个实根。二、程序功能利用二分法，求方程的根三、算法步1:设区间为[ab],先把(a+b)/2代入f(x)，若为0则根x=(a+b)/2，若不为0看符号。步2：若f(a)*f((a+b)/2)0则根在区间（a(a+b)/2）内，否则根在区间（(a+b)/2b）内，继续在有根区间内二分，一直到达到精度。步3：循环步骤1和3，直到步骤2成立８步4：取有根区间的中点为方程的根。四、重要标识符说明a,b：是区间的端点；emg：达到的精度；k：循环次数；x：所求近似值；Fa：将a代入方程的值；fab：将(a+b)/2代入方程的值；feval()：是求值函数；abs（）：是求绝对值函数；@：函数调用；五、程序运行实例在matlab6.5环境中，运行程序，结果如图：９六、源程序头文件functionf=fun(x)f=x^3-x-1;源程序%%%%二分法%%%%function[x,k]=demimethod(a,b,f,emg);fa=feval(f,a);fab=feval(f,(a+b)/2);k=0;whileabs(b-a)emg;iffab==0;x=(a+b)/2;return;elseiffa*fab0b=(a+b)/2;elsea=(a+b)/2;endfa=feval(f,a);fab=feval(f,(a+b)/2);k=k+1;endx=(a+b)/2;七、个人实验总结第一次用到了函数调用的方法，起初以为把它们写到一个文件里即可，后来才知道在程序中若要用到函数调用是不能在一个m文件中再定义的，需要再建立一个m文件，运行时用符号＠标识，并为变量赋值。基本上已经熟悉了MATLAB的常用功能。１０实验五高斯列主元素消去法解方程组一、实验题目解方程组：0.6428x1+0.3475x2-0.8468x3=0.41270.3475x1+1.8423x2+0.4759x3=1.7321-0.8468x1+0.4759x2+1.2147x3=-0.8621二、程序功能用高斯列主元素消去法解方程组三、算法步1：选主元即从增广矩阵第一列中选择绝对值最大的数所在的行与第一行互换，并将新矩阵的第二行和第三行的第一列数化为0，形成新的矩阵；步2：从新矩阵的第二行和第三行的第二列中选择绝对值最大的换到第二行，第三行第二列数化为0，形成新矩阵；步3：回代，先求出x3,回代求出x2,进而可得x1.四、重要标识符说明A，b：系数矩阵N：b的长度X：一个长度为n，值为0的列向量C：长度为n，值为0的横向量d1：数0A(i,i)：矩阵A中对角线的数A(j,i)：矩阵A中第j行第i列的元素zeros(n,m)：生成m*n阶0矩阵；abs：求绝对值；１１五、程序运行实例在matlab6.5环境中，运行程序，结果如图：六、源程序function[x]=Gauss(A,b)A=[0.64280.3475-0.8468;0.34751.84230.4759;-0.84680.47591.2147];b=[0.4127;1.7321;0.8621];n=length(b);X=zeros(n,1);c=zeros(1,n);d1=0;fori=1:n-1max=abs(A(i,i));m=i;forj=i+1:nifmaxabs(A(j,i));max=abs(A(j,i));m=j;endendif(m~=i)fork=i:nc(k)=A(i,k);１２A(i,k)=A(m,k);A(m,k)=c(k);endd1=b(i);b(i)=b(m);b(m)=d1;endfork=i+1:nforj=i+1:nA(k,j)=A(k,j)-A(i,j)*A(k,i)/A(i,i);endb(k)=b(k)-b(i)*A(k,i)/A(i,i);A(k,i)=0;endend%%%回代求解%%%x(n)=b(n)/A(n,n);fori=n-1:-1:1;sum=0;forj=i+1:nsum=sum+A(i,j)*x(j);endx(i)=(b(i)-sum)/A(i,i);end七、个人实验总结实验代码很长，其实到最后我也不是很明白，尤其是那些嵌套的循环语句，还好老师给了试验源代码，我们只是在此基础上稍作了些修改，最后得出了正确答案。经过五个实验的运行，已经对MATLAB环境有了进一步的了解，以后若有用到的地方也不会无从下手。