数值分析实验报告Matlab仿真(共20页).doc
精选优质文档-倾情为你奉上数值分析实验报告学院：电气工程与自动化学院专业：控制理论与控制工程姓名：李亚学号：2014 年 12 月24日实验一 函数插值方法 一、目的和意义 1、 学会常用的插值方法，求函数的近似表达式，以解决其它实际问题； 2、 明确插值多项式和分段插值多项式各自的优缺点； 3、 熟悉插值方法的程序编制； 4、 如果绘出插值函数的曲线，观察其光滑性。二、实验原理 1、 Lagrange插值公式 编写出插值多项式程序； 2、 给出插值多项式或分段三次插值多项式的表达式； 三、实验要求对于给定的一元函数的n+1个节点值。试用Lagrange公式求其插值多项式或分段二次Lagrange插值多项式。数据如下： （1） 0.4 0.55 0.65 0.80 0.95 1.05 0.41075 0.578150.696750.90 1.00 1.25382 求五次Lagrange多项式，计算,的值。（提示：结果为, ） （2） 1 2 3 4 5 6 7 0.368 0.135 0.050 0.018 0.007 0.002 0.001 试构造Lagrange多项式，和分段三次插值多项式，计算的,值。（提示：结果为, ）四、实验过程1进入matlab开发环境；2根据实验内容和要求编写程序，程序如下所示，程序通过运用function函数编写，生成.m文件。调用时只需要在命令窗口调用y=Lagrange(A,input)就可以实现任意次数拉格朗日插值法求解。function y=Lagrange(A,input)a,b=size(A);x=input;y=0;for j=1:a Mj=1; Nj=1; for k=1:a if(k=j) continue; end Mj=Mj*(x-A(k,1); Nj=Nj*(A(j,1)-A(k,1); end y=y+A(j,2)*Mj/Nj;end3调试程序并运行程序；调用拉格朗日脚本文件对以上两个表格数据求解，表格一对应MATLAB向量A;表格二对应向量I。在命令窗口调用y=Lagrange(A,input)，求解如下面截图。图1 表一数据的解图2 表二数据的解4实验总结 通过对插值法算法编程，加深了对插值方法的理解，熟悉了MATLAB编写脚本函数。通过计算机求解，能更加方便快捷求解。 实验二 函数逼近与曲线拟合 一、目的和意义 1、掌握曲线拟合的最小二乘法； 2、最小二乘法亦可用于解超定线代数方程组； 3、探索拟合函数的选择与拟合精度间的关系。二、实验原理 对于给定的测量数据(xi,fi)(i=1,2,，n)，设函数分布为特别的，取为多项式 (j=0, 1,，m)则根据最小二乘法原理，可以构造泛函令 (k=0, 1,，m)则可以得到法方程求该解方程组，则可以得到解，因此可得到数据的最小二乘解三、实验要求 1、用最小二乘法进行曲线拟合； 2、近似解析表达式为；3、打印出拟合函数，并打印出与的误差，； 4、另外选取一个近似表达式，尝试拟合效果的比较； 5、绘制出曲线拟合图。 四、实验步骤： 1进入matlab开发环境；2根据实验内容和要求编写程序如下；代码一公式S(x)=a1*t+a2*t2+a3*t3；代码二公式S(x)=a2*t2+a3*t3+a4*t4.代码一：function error=mintwomultiply(A)%S(x)=a1*t+a2*t2+a3*t3a,b=size(A);M=zeros(3);N=zeros(3,1);error=0;for i=1:a M(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);end%a1,a2,a3=solve(M,N)I=MN;for i=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;endhold on;plot(A(:,1),A(:,3),'r','LineWidth',2);plot(A(:,1),A(:,2),'b','LineWidth',2);legend('原始图像',拟合图像',2);hold off;代码二：function error=mintwomultiply2(A)%S(x)=a2*t2+a3*t3+a4*t4a,b=size(A);M=zeros(3);N=zeros(3,1);error=0;for i=1:a M(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); N(3,1)=N(3,1)+A(i,1)4*A(i,2);end%a1,a2,a3=solve(M,N)I=MN;for i=1:a A(i,3)=I(1,1)*A(i,1)2+I(2,1)*A(i,1)3+I(3,1)*A(i,1)4 error=error+(A(i,3)-A(i,2)2;endhold on;plot(A(:,1),A(:,3),'r','LineWidth',2);plot(A(:,1),A(:,2),'b','LineWidth',2);legend('原始图像',拟合图像',2);hold off;3调试程序并运行程序；在某冶炼过程中，根据统计数据的含碳量与时间关系，试求含碳量与时间t的拟合曲线。t(分)0 5 10 15 20 25 30 35 40 45 50 55 0 1.27 2.16 2.86 3.44 3.87 4.15 4.37 4.51 4.58 4.02 4.64 根据以上数据调用代码一和代码二，结果如下图，返回的结果为error为误差，拟合图像蓝的为原始图像，红的为拟合图像，具体图像见下图:图1 公式一误差图2 公式一拟合图图3 公式二误差图4 公式二拟合图4实验总结从不同公式产生的误差来看，公式一明显误差小；从产生的拟合图像来看公式一更稳定，更平稳。由此可得，公式的选取对结果的影响很大。实验三 线方程组的直接解法一、目的和意义 1、通过该课题的实验，体会模块化结构程序设计方法的优点； 2、运用所学的计算方法，解决各类线性方程组的直接算法； 3、提高分析和解决问题的能力，做到学以致用；4、学习Gauss消去法的原理；5、了解列主元的意义；6、确定什么时候系数阵要选主元。二 实验数学原理：由于一般线性方程在使用Gauss消去法求解时，从求解的过程中可以看到，若=0，则必须进行行交换，才能使消去过程进行下去。有的时候即使0，但是其绝对值非常小，由于机器舍入误差的影响，消去过程也会出现不稳定得现象，导致结果不正确。因此有必要进行列主元技术，以最大可能的消除这种现象。这一技术要寻找行r，使得并将第r行和第k行的元素进行交换，以使得当前的的数值比0要大的多。这种列主元的消去法的主要步骤如下：1. 消元过程对k=1,2,n-1,进行如下步骤。1) 选主元，记若很小，这说明方程的系数矩阵严重病态，给出警告，提示结果可能不对。2) 交换增广阵A的r，k两行的元素。 (j=k,n+1)3) 计算消元 (i=k+1,n; j=k+1,n+1)2. 回代过程对k= n, n-1,1,进行如下计算至此，完成了整个方程组的求解。三、实验要求 1、对下述方程组分别利用Gauss顺序消去法与Gauss列主元消去法；2、应用结构程序设计编出通用程序； 3、比较计算结果，分析数值解误差的原因； 4、设线性方程组如下 方程的解：四、实验要求 1进入matlab开发环境；2根据实验内容和要求编写程序；代码分为两个，Gauss顺序消去法与Gauss列主元消去法；代码一：Gauss顺序消去法function x=orderGauss(A,B)m,n=size(A);k=rank(A);if m=n&&m=k for i=1:n-1 for j=i+1:n B(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); end end x(n,1)=B(n,1)/A(n,n); for i=n-1:-1:1 sum=0; for j=i+1:n sum=sum+A(i,j)*x(j,1); end x(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)m n=size(A);if m=1&&n=1 rows=1; columns=1; temp=A(1,1); for i=2:m if abs(A(i,1)>abs(temp) temp=A(i,1); rows=i; end endelseif m=1&&n=1 rows=1; columns=1; temp=A(1,1); for i=2:n if abs(A(1,i)>abs(temp) temp=A(1,i); columns=i; end end else rows=1; columns=1; temp=A(1,1); for i=1:m for j=1:n if abs(A(i,j)>abs(temp) temp=A(i,j); rows=i; columns=j; end end endend第二部分：高斯消元法求解function x=Gauss(A,B)m,n=size(A);k=rank(A);if m=n&&m=k for i=1:n-1 a,c=findmax(A(i:m,i); a=a+i-1; if c=0 continue; end if a=i temp=A(i,:); A(i,:)=A(a,:); A(a,:)=temp; temp1=B(i,1); B(i,1)=B(a,1); B(a,1)=temp1; end for j=i+1:n B(j,:)=B(j,1)-A(j,i)*B(i,1)/c; A(j,:)=A(j,:)-A(j,i)*A(i,:)/c; end end x(n,1)=B(n,1)/A(n,n); for i=n-1:-1:1 sum=0; for j=i+1:n sum=sum+A(i,j)*x(j,1); end x(i,1)=(B(i,1)-sum)/A(i,i); endend3调试程序并运行程序；用Gauss顺序消去法与Gauss列主元消去法求解方程组的解，结果如下图： 图1 Gauss顺序消去法的解 图2 Gauss列主元消去法的解4实验总结对比试验结果，两种方法的解和方程的解都一直。这是由于计算机有足够的精度，如果没有足够的精度，误差就容易被放大。产生不正确的结果。列主元消去法可以克服这种缺陷，求出相对精确的解。实验四 解线性方程组的迭代法一、目的和意义 1、通过上机计算体会迭代法求解线性方程组的特点，并能和消去法比较； 2、运用所学的迭代法算法，解决各类线性方程组，编出算法程序； 3、体会上机计算时，终止步骤或k >（给予的迭代次数），对迭代法敛散性的意义； 4、体会初始解，松弛因子的选取，对计算结果的影响。二、实验原理设线性方程组 (1)的系数矩阵A可逆且主对角元素均不为零,令 并将A分解成 (2)从而(1)可写成 令 其中. (3)以为迭代矩阵的迭代法(公式) (4)称为雅可比(Jacobi)迭代法(公式),用向量的分量来表示,(4)为 (5)其中为初始向量.三、要求 1、体会迭代法求解线性方程组，并能与消去法做以比较；四、实验过程 1进入matlab开发环境2根据实验内容和要求编写程序此代码分为两部分，即两个功能函数。y=Jacobicondition(A)和X=Jacobi(A,B,X,input)。前面函数为判断条件，判断矩阵A是否可以进行雅克比迭代。后面X=Jacobi(A,B,X,input)进行迭代求解，具体代码如下：代码一：function y=Jacobicondition(A)a=size(A);k=rank(A);D=A;L=A;U=A;for i=1:a for j=1:k if i=j D(i,j)=0; end if i>j L(i,j)=-L(i,j); else L(i,j)=0; end if i<j U(i,j)=-U(i,j); else U(i,j)=0; end endendJ=D(L+U);m n=eig(J);for i=1:a if n(i,i)>=1; y=1; break; else y=0; endend代码二：function X=Jacobi(A,B,X,input)y=Jacobicondition(A);a b=size(A);k=rank(A);temp=zeros(a,input+1);temp(:,1)=X(:,1)+temp(:,1);if a=b&&a=k&&y=0 for i=1:input for j=1:a sum=0; for m=1:a if j=m continue; end sum=sum+A(j,m)*temp(m,i); end temp(j,i+1)=(B(j,1)-sum)/A(j,j); end end X=temp(:,input+1);else disp('不满足雅克比条件');end3调试程序并运行程序图1 方程组的系数矩阵 图2 迭代次数为5,10,30,50的解专心-专注-专业
