数值分析实验报告Matlab仿真

数值分析实验报告学院：电气工程与自动化学院专业：控制理论与控制工程姓名：李亚学号：61201401622014年12月24日1实验一函数插值方法一、目的和意义1、学会常用的插值方法，求函数的近似表达式，以解决其它实际问题；2、明确插值多项式和分段插值多项式各自的优缺点；3、熟悉插值方法的程序编制；4、如果绘出插值函数的曲线，观察其光滑性。二、实验原理1、Lagrange插值公式00,()nninkkiikkixxLxyxx编写出插值多项式程序；2、给出插值多项式或分段三次插值多项式的表达式；三、实验要求对于给定的一元函数)(xfy的n+1个节点值(),0,1,,jjyfxjn。试用Lagrange公式求其插值多项式或分段二次Lagrange插值多项式。数据如下：（1）jx0.40.550.650.800.951.05jy0.410750.578150.696750.901.001.25382求五次Lagrange多项式5L()x，计算(0.596)f,(0.99)f的值。（提示：结果为(0.596)0.625732f,(0.99)1.05423f）（2）jx1234567jy0.3680.1350.0500.0180.0070.0020.001试构造Lagrange多项式6L()x，和分段三次插值多项式，计算的(1.8)f,(6.15)f值。（提示：结果为(1.8)0.164762f,(6.15)0.001266f）2四、实验过程1．进入matlab开发环境；2．根据实验内容和要求编写程序，程序如下所示，程序通过运用function函数编写，生成.m文件。调用时只需要在命令窗口调用y=Lagrange(A,input)就可以实现任意次数拉格朗日插值法求解。functiony=Lagrange(A,input)[a,b]=size(A);x=input;y=0;forj=1:aMj=1;Nj=1;fork=1:aif(k==j)continue;endMj=Mj*(x-A(k,1));Nj=Nj*(A(j,1)-A(k,1));endy=y+A(j,2)*Mj/Nj;end3．调试程序并运行程序；调用拉格朗日脚本文件对以上两个表格数据求解，表格一对应MATLAB向量A;表格二对应向量I。在命令窗口调用y=Lagrange(A,input)，求解如下面截图。3图1表一数据的解图2表二数据的解4．实验总结通过对插值法算法编程，加深了对插值方法的理解，熟悉了MATLAB编写脚本函数。通过计算机求解，能更加方便快捷求解。4实验二函数逼近与曲线拟合一、目的和意义1、掌握曲线拟合的最小二乘法；2、最小二乘法亦可用于解超定线代数方程组；3、探索拟合函数的选择与拟合精度间的关系。二、实验原理对于给定的测量数据(xi,fi)(i=1,2,…，n)，设函数分布为mjjjxaxy0)()(特别的，取)(xj为多项式jjxx)((j=0,1,…，m)则根据最小二乘法原理，可以构造泛函nimjijjimxafaaaH1010))((),,,(令0kaH(k=0,1,…，m)则可以得到法方程),(),(),(),(),(),(),(),(),(),(),(),(1010101111000100mmmmmmmmfffaaa求该解方程组，则可以得到解maaa,,,10，因此可得到数据的最小二乘解mjjjxaxf0)()(三、实验要求1、用最小二乘法进行曲线拟合；52、近似解析表达式为23123()tatatat；3、打印出拟合函数()t，并打印出()jt与()jyt的误差，1,2,,12j；4、另外选取一个近似表达式，尝试拟合效果的比较；5、绘制出曲线拟合图。四、实验步骤：1．进入matlab开发环境；2．根据实验内容和要求编写程序如下；代码一公式S(x)=a1*t+a2*t^2+a3*t^3；代码二公式S(x)=a2*t^2+a3*t^3+a4*t^4.代码一：functionerror=mintwomultiply(A)%S(x)=a1*t+a2*t^2+a3*t^3[a,b]=size(A);M=zeros(3);N=zeros(3,1);error=0;fori=1:aM(1,1)=M(1,1)+A(i,1)*A(i,1);M(1,2)=M(1,2)+A(i,1)*A(i,1)^2;M(2,1)=M(1,2);M(1,3)=M(1,3)+A(i,1)*A(i,1)^3;M(3,1)=M(1,3);M(2,2)=M(2,2)+A(i,1)^2*A(i,1)^2;M(2,3)=M(2,3)+A(i,1)^2*A(i,1)^3;M(3,2)=M(2,3);M(3,3)=M(3,3)+A(i,1)^3*A(i,1)^3;N(1,1)=N(1,1)+A(i,1)*A(i,2);N(2,1)=N(2,1)+A(i,1)^2*A(i,2);N(3,1)=N(3,1)+A(i,1)^3*A(i,2);end6%[a1,a2,a3]=solve(M,N)I=M\N;fori=1:aA(i,3)=I(1,1)*A(i,1)+I(2,1)*A(i,1)^2+I(3,1)*A(i,1)^3;error=error+(A(i,3)-A(i,2))^2;endholdon;plot(A(:,1),A(:,3),'r','LineWidth',2);plot(A(:,1),A(:,2),'b','LineWidth',2);legend('原始图像',拟合图像',2);holdoff;代码二：functionerror=mintwomultiply2(A)%S(x)=a2*t^2+a3*t^3+a4*t^4[a,b]=size(A);M=zeros(3);N=zeros(3,1);error=0;fori=1:aM(1,1)=M(1,1)+A(i,1)^2*A(i,1)^2;M(1,2)=M(1,2)+A(i,1)^2*A(i,1)^3;M(2,1)=M(1,2);M(1,3)=M(1,3)+A(i,1)^2*A(i,1)^4;M(3,1)=M(1,3);M(2,2)=M(2,2)+A(i,1)^3*A(i,1)^3;M(2,3)=M(2,3)+A(i,1)^3*A(i,1)^4;M(3,2)=M(2,3);M(3,3)=M(3,3)+A(i,1)^4*A(i,1)^4;N(1,1)=N(1,1)+A(i,1)^2*A(i,2);N(2,1)=N(2,1)+A(i,1)^3*A(i,2);7N(3,1)=N(3,1)+A(i,1)^4*A(i,2);end%[a1,a2,a3]=solve(M,N)I=M\N;fori=1:aA(i,3)=I(1,1)*A(i,1)^2+I(2,1)*A(i,1)^3+I(3,1)*A(i,1)^4error=error+(A(i,3)-A(i,2))^2;endholdon;plot(A(:,1),A(:,3),'r','LineWidth',2);plot(A(:,1),A(:,2),'b','LineWidth',2);legend('原始图像',拟合图像',2);holdoff;3．调试程序并运行程序；在某冶炼过程中，根据统计数据的含碳量与时间关系，试求含碳量与时间t的拟合曲线。t(分)05101520253035404550554(10)y01.272.162.863.443.874.154.374.514.584.024.64根据以上数据调用代码一和代码二，结果如下图，返回的结果为error为误差，拟合图像蓝的为原始图像，红的为拟合图像，具体图像见下图:8图1公式一误差图2公式一拟合图9图3公式二误差图4公式二拟合图4．实验总结从不同公式产生的误差来看，公式一明显误差小；从产生的拟合图像来看公式一更稳定，更平稳。由此可得，公式的选取对结果的影响很大。10实验三线方程组的直接解法一、目的和意义1、通过该课题的实验，体会模块化结构程序设计方法的优点；2、运用所学的计算方法，解决各类线性方程组的直接算法；3、提高分析和解决问题的能力，做到学以致用；4、学习Gauss消去法的原理；5、了解列主元的意义；6、确定什么时候系数阵要选主元。二实验数学原理：由于一般线性方程在使用Gauss消去法求解时，从求解的过程中可以看到，若)1(kkka=0，则必须进行行交换，才能使消去过程进行下去。有的时候即使)1(kkka0，但是其绝对值非常小，由于机器舍入误差的影响，消去过程也会出现不稳定得现象，导致结果不正确。因此有必要进行列主元技术，以最大可能的消除这种现象。这一技术要寻找行r，使得)1()1(max||kikkikrkaa并将第r行和第k行的元素进行交换，以使得当前的)1(kkka的数值比0要大的多。这种列主元的消去法的主要步骤如下：1.消元过程对k=1,2,…,n-1,进行如下步骤。1)选主元，记ikkirkaamax||若||rka很小，这说明方程的系数矩阵严重病态，给出警告，提示结果可能不对。2)交换增广阵A的r，k两行的元素。kjrjaa(j=k,…,n+1)3)计算消元kkkjikijijaaaaa/(i=k+1,…,n;j=k+1,……,n+1)112.回代过程对k=n,n-1,…,1,进行如下计算)/(11,nkjkkjkjnkkaxaax至此，完成了整个方程组的求解。三、实验要求1、对下述方程组分别利用Gauss顺序消去法与Gauss列主元消去法；2、应用结构程序设计编出通用程序；3、比较计算结果，分析数值解误差的原因；4、设线性方程组如下123456789104231210000865365010042213210310215131194426167332386857172635021342530116101191734212246271392012400183248631xxxxxxxxxx5123234613381921方程的解：(1,1,0,1,2,0,3,1,1,2)Tx四、实验要求1．进入matlab开发环境；2．根据实验内容和要求编写程序；代码分为两个，Gauss顺序消去法与Gauss列主元消去法；代码一：Gauss顺序消去法functionx=orderGauss(A,B)[m,n]=size(A);k=rank(A);ifm==n&&m==kfori=1:n-112forj=i+1:nB(j,:)=B(j,1)-A(j,i)*B(i,1)/A(i,i);A(j,:)=A(j,:)-A(j,i)*A(i,:)/A(i,i);endendx(n,1)=B(n,1)/A(n,n);fori=n-1:-1:1sum=0;forj=i+1:nsum=sum+A(i,j)*x(j,1);endx(i,1)=(B(i,1)-sum)/A(i,i);endend代码二：Gauss列主元消去法此代码分为两部分，即两个功能函数。[rows,columns,temp]=findmax(A)和x=Gauss(A,B)。前面函数寻找向量中的最大值，返回所在的行，列，具体值。后面再findmax基础进行列主元迭代，代码如下：第一部分：寻找向量，矩阵最大值function[rows,columns,temp]=findmax(A)[mn]=size(A);ifm~=1&&n==1rows=1;columns=