2016春计算方法实验指导-作业.doc

计算方法实验指导材料目 录第1次 解线性方程组的直接解法- 2 -1.1 例题- 2 -1.2 Matlab解线性方程组常用命令介绍- 6 -1.3 实验作业- 7 -第2次 解线性方程组的迭代解法- 8 -2.1 例题- 8 -2.2 Matlab迭代解法常用函数介绍- 12 -2.3 实验作业- 12 -第3次 非线性方程求根- 13 -3.1 例题- 13 -3.2 Matlab非线性方程求根的命令- 22 -3.3 实验作业- 22 -第4次 插值法234.1 例题234.2 Matlab插值函数介绍294.3 实验作业30第5次 利用最小2乘法进行曲线拟合315.1 例题315.2 Matlab数据拟合命令介绍345.3 实验作业35第6次 数值积分与数值微分366.1 例题366.2 Matlab数值积分函数介绍426.3 实验作业43第1次 解线性方程组的直接解法Gauss消元法算法如下：Gauss列主元消元法算法：在以上的算法中，增加调换行的步骤即可。在每步消元之前要按列选出主元。1.1 例题例1.1 用Gauss消元法解方程组：解：直接建立求解该方程组的M文件Gauss.m如下% 求解例题1.1% 高斯法求解线性方程组Ax=b% A为输入矩阵系数,b为方程组右端系数% 方程组的解保存在x变量中% 先输入方程系数A=1 2 3;2 7 5;1 4 9;b=1 6 -3;m,n=size(A);%检查系数正确性if m=n error(矩阵A的行数和列数必须相同); return;endif m=size(b) error(b的大小必须和A的行数或A的列数相同); return;end%再检查方程是否存在唯一解if rank(A)=rank(A,b) error(A矩阵的秩和增广矩阵的秩不相同,方程不存在唯一解); return;end%这里采用增广矩阵行变换的方式求解c=n+1;A(:,c)=b; %消元过程for k=1:n-1A(k+1:n, k:c)=A(k+1:n, k:c)-(A(k+1:n,k)/ A(k,k)*A(k, k:c); End%回代结果x=zeros(length(b),1); x(n)=A(n,c)/A(n,n);for k=n-1:-1:1x(k)=(A(k,c)-A(k,k+1:n)*x(k+1:n)/A(k,k);end%显示计算结果disp(x=);disp(x);直接运行上面的M文件或在Matlab命令窗口中直接输入Gauss即可得出结果.在Matlab命令窗口中输入Gauss得出结果如下: Gaussx= 2.0000 1.0000 -1.0000扩展：Matlab求解线性方程的几种命令如下（方程组的一般形式可用矩阵和向量表示成，但运用下列方法的前提必须保证所求解的方程为恰定方程，即方程组存在唯一解）运用求逆思想：: 或 ；左除法：原理上是运用高斯消元法求解，但Matlab在实际执行过程中是通过分解法进行的（即先将矩阵A作分解,再回代计算）：；符号矩阵法：这种计算方法最接近精确值，但计算速度最慢：；将矩阵施行初等行变换化成行简化阶梯形的办法：可以这样实现之: ; 上面四种常用的办法示例如下： A=1 2 3;2 7 5;1 4 9 % 上面示例方程组系数A = 1 2 3 2 7 5 1 4 9 b=1 6 -3 % 方程组右端的系数b = 1 6-3 x1_1=inv(A)*b,x1_2=A(-1)*b % 方法一，求逆思想x1_1 = 2.0000 1.0000 -1.0000x1_2 = 2.0000 1.0000 -1.0000 x2=Ab % 方法二，左除思想x2 = 2 1-1 x3=sym(A)sym(b) % 方法三，符号法 x3 = 2 1 -1 C=A,b,rref(C) % 方法四，行简化阶梯形思想,最后输出结果的一列为解C = 1 2 3 1 2 7 5 6 1 4 9 -3ans = 1 0 0 2 0 1 0 1 0 0 1 -1例1.2 ，求解：输入矩阵：x=1 0.5 -0.3 x = 1.0000 0.5000 -0.3000 计算其1-范数norm_1=norm(x,1) norm_1 = 1.8000 计算其2-范数norm_2=norm(x) norm_2 = 1.1576 计算其无穷大范数norm_inf=norm(x,inf) norm_inf = 1 还可以计算其无穷小范数(即各分量绝对值中的最小值)norm_minusInf=norm(x,-inf) norm_minusInf = 0.3000 例1.3 设矩阵，求解：先输入矩阵:A=-2 1 0;1 -2 1;0 1 -2 A = -2 1 0 1 -2 1 0 1 -2 求A的1-范数(列和范数):norm_1=norm(A,1) norm_1 = 4 求解A的无穷大范数(行和范数):norm_inf=norm(A,inf) norm_inf = 4 求A的2-范数(最大特征值):norm_2=norm(A,2) norm_2 = 3.4142 还可以求解A的F-范数:norm_F=norm(A,fro) norm_F = 4 谱半径可以通过按其概念进行计算:对其特征值的绝对值取最大值即可:R_A=max(abs(eig(A) R_A = 3.4142 1.2 Matlab解线性方程组常用命令介绍1. 求秩命令rank在解线性方程组时，为了判断是否存在解，我们应先判断矩阵的秩调用格式为:c=rank(A)2. 矩阵零空间向量null当线性方程组的秩时，方程组有无穷多个解，使用x=null(A)可得到满足的一个解向量，可为符号矩阵或数值矩阵：x=null(A) 或 x=null(sym(A)3. 方阵的LU分解命令lu调用格式为：L,U=lu(A)L为准下三角矩阵,U为上三角矩阵,满足A=LUL,U,P=lu(A)L为下三角矩阵,U为上三角矩阵,P为变换方阵元素位置的换位阵,满足PA=PL其它调用格式请用help lu获得更多信息4. Cholsky分解cholL=chol(A)其中L为一个下三角矩阵,满足. 必须为正定矩阵5. 条件数condc=cond(A,p)A可为向量或矩阵,P取值为下列之一:1向量或矩阵的返回1条件数.2返回向量或矩阵的2条件数.inf返回向量或矩阵的条件数.fro返回向量或矩阵的F条件数.6. 奇异值分解svdU,S,V=svd(A)将A分解为正交矩阵U,对角矩阵S和正交矩阵V的乘积,使得A=USVT.7. 线性方程组的符号解linsolveX=linsolve(A,b)它等价于X=sym(A)sym(b).返回方程组的符号解1.3 实验作业实验名称：解线性方程组的列主元素高斯消去法实验目的：通过数值实验，从中体会解线性方程组选主元的必要性，以及方程组系数矩阵和右端向量的微小变化对解向量的影响实验内容：解线性方程组：实验要求：（1） 用Matlab语言编写程序，用列主元高斯消去法求解上述方程组，输出Ax=b中矩阵A、b、解向量x及detA（2） 将方程组中系数3.01改为3.00，0.987改为0.990，用列主元高斯消去法求解系数变换后的方程组，输出列主元行交换次序、解向量x及detA，并与（1）中结果比较（3）用MATLAB的内部函数inv求出系数矩阵的逆矩阵，再输入命令x=inv(A)*b,即可求出上述各个方程组的解，并与列主元高斯消去法求出的解进行比较，体会选主元的方法具有良好的数值稳定性用MATLAB的内部函数det求出系数行列式的值，并与（1）、（2）、中输出的系数行列式的值进行比较第2次 解线性方程组的迭代解法2.1 例题Jacobi迭代法算法：流程图如下。高斯-塞德尔迭代算法：计算步骤与流程图与雅可比迭代法大致相同，只是一旦求出变元 的某个新值 后, 就改用新值 替代老值 进行这一步剩下的计算。例2.1 用Jacobi迭代法解以下方程：解：编制迭代计算的M文件程序如下：% Jacobi迭代法求解例3.1% A为方程组的增广矩阵clc;A=10 -2 -1 3;-2 10 -1 15;-1 -2 5 10MAXTIME=50;%最多进行50次迭代eps=1e-5;%迭代误差n,m=size(A); x=zeros(n,1);%迭代初值 y=zeros(n,1); k=0; % 进入迭代计算 disp(迭代过程X的值情况如下:) disp(X=); while 1 disp(x); for i=1:1:n s=0.0; for j=1:1:n if j=i s=s+A(i,j)*x(j); end y(i)=(A(i,n+1)-s)/A(i,i); end end for i=1:1:n maxeps=max(0,abs(x(i)-y(i); % 检查是否满足迭代精度要求 end if maxepsMAXTIME % 超过最大迭代次数退出 error(超过最大迭代次数,退出); return; end end运行该程序结果如下：A = 10 -2 -1 3 -2 10 -1 15 -1 -2 5 10迭代过程X的值情况如下：X= 0 0 0 0.3000 1.5000 2.0000 0.8000 1.7600 2.6600 0.9180 1.9260 2.8640 0.9716 1.9700 2.9540 0.9894 1.9897 2.9823 0.9962 1.9961 2.9938 0.9986 1.9986 2.9977 0.9995 1.9995 2.9992 0.9998 1.9998 2.9997 0.9999 1.9999 2.9999 1.0000 2.0000 3.0000 1.0000 2.0000 3.0000容易看出迭代计算最后结果为：例2.2 用Gauss-Seidel迭代法计算例2.1并作比较.解：编制求解程序Gauss_Seidel.m如下：% Gauss_Seidel.m% Gauss_Seidel迭代法求解例2.2% A为方程组的增广矩阵clc;format long;A=10 -2 -1 3;-2 10 -1 15;-1 -2 5 10n,m=size(A);% 最多进行50次迭代Maxtime=50;%控制误差Eps=10E-5;% 初始迭代值x=zeros(1,n);disp(x=);% 迭代次数小于最大迭代次数,进入迭代for k=1:Maxtime disp(x); for i=1:n s=0.0; for j=1:n if i=j s=s+A(i,j)*x(j);%计算和 end end x(i)=(A(i,n+1)-s)/A(i,i);%求出此时迭代的值 end% 因为方程的精确解为整数,所以这里将迭代结果向整数靠近的误差作为判断迭代是否停止的条件 if sum(x-floor(x).2)Gauss_Seidel A = 10 -2 -1 3 -2 10 -1 15 -1 -2 5 10x= 0 0 0 0.300000000000000 1.560000000000000 2.684000000000000 0.880400000000000 1.944480000000000 2.953872000000000 0.984283200000000 1.992243840000000 2.993754176000000 0.997824185600000 1.998940254720000 2.999140939008000 0.999702144844800 1.999854522869760 2.999882238116864 0.999959128385638 1.999980049488814 2.999983845472654 0.999994394445028 1.999997263436271 2.999997784263514 0.999999231113606 1.999999624649072 2.999999696082350 0.999999894538049 1.999999948515845 2.999999958313948 0.999999985534564 1.999999992938308 2.999999994282236 0.999999998015885 1.999999999031401 2.999999999215737 0.999999999727854 1.999999999867145 2.999999999892429 0.999999999962672 1.999999999981778 2.999999999985246 0.999999999994880 1.999999999997501 2.999999999997976 0.999999999999298 1.999999999999657 2.999999999999722 0.999999999999904 1.999999999999953 2.999999999999962 0.999999999999987 1.999999999999994 2.999999999999995 0.999999999999998 1.999999999999999 2.999999999999999 1.000000000000000 2.000000000000000 3.000000000000000迭代结果:X = 1 2 3 可见对此题Gauss_Seidel法的收敛速度还是很快的2.2 Matlab迭代解法常用函数介绍1. 计算范数命令normnorm(A,option)当A为向量V时:norm(V,1)为A的1-范数.norm(V,2)或norm(A)为A的1-范数.notm(V,P) =sum(abs(V).P)(1/P).norm(V,inf) =max(abs(V)为A的无穷范数.norm(V,-inf) = min(abs(V)为A的最小范数.当A为矩阵A时:norm(A,1)为A的列和范数.norm(A,2)或norm(A)为A的谱范数.norm(A,inf)为A的行和范数.norm(A,fro)时为A的F-范数.2. 计算矩阵的谱半径 r=max(abs(eig(A)2.3 实验作业实验名称：解线性方程组的迭代法实验目的：1掌握解线性方程组的雅可比迭代和高斯-塞德尔迭代算法；2培养编程与上机调试能力实验内容：应用雅可比迭代和高斯-塞德尔迭代算法解线性方程组：实验要求：（1） 用Matlab语言编写雅可比（Jacobi）迭代法程序，并且选择不同的迭代次数，观察输出结果；（2）用Matlab语言编写高斯-塞德尔(Gauss-Seidel)迭代法程序，并且选择不同的迭代次数，观察输出结果第3次 非线性方程求根二分法：算法流程图迭代法：算法流程图牛顿迭代法：算法流程图3.1 例题例3.1 判别方程的实根存在区间，要求区间长度不大于1，并求出最小正根的近似值,精度解：先用画图的方法来粗略估计其根的范围ezplot(x3-3*x+1);%亦可用fplot命令grid on; 观察图可知其根分布区间大概分布在(-2,-1),(0,1)和（1,2）中.编写二分法求解最小正根的近似值程序如下:format long;f=inline(x3-3*x+1);a=0;b=1;Eps=1E-5;for k=1:50A(k)=a;B(k)=b;ya=feval(f,a);yb=feval(f,b);temp=(a+b)/2;X(k)=temp;yt=feval(f,temp);F(k)=yt;if abs(yt)Eps break;endif yt*ya0 a=a;b=temp;elseif yt*yb1E10;%此为x 发散break;end;end;i, x 运行结果如下:f = Inline function: f(x) = x-x3-4*x2+10x= -0.8750 6.7324 -469.7200 1.0275e+008 -1.0849e+024 1.2771e+072i = 6x = 1.2771e+072迭代6次后x的值大得令人吃惊,表明构造的式子并不收敛(2) :求解程序如下:f=inline(10/x-4*x)(1/2)disp(x=);x=feval(f,1.5);disp(x);Eps=1E-5;i=1;while 1 x0=x; i=i+1; x=feval(f,x); disp(x); if isreal(x) % 不是实数不进行迭代 break; end if x1E10; % 发散不进行迭代 break; end; if abs(x-x0)1E10; break; end; if abs(x-x0)x=sym(x)f=sym(x*sin(x)-0.5)求符号函数导数df=diff(f,x)构造Newton迭代公式FX=x-f/df %以下为输入结果x =xf =x*sin(x)-0.5df =sin(x)+x*cos(x)FX =x-(x*sin(x)-.5)/(sin(x)+x*cos(x) 进入迭代计算,计算前10次迭代值：format long;Fx=inline(FX);x0=0.7;for i=1:10disp(x0); x0=feval(Fx,x0);endx0format short; 输出x的中间结果为： 0.700000000000000 0.741579619191183 0.740841172608556 0.740840955095510 0.740840955095491 0.740840955095491 0.740840955095491 0.740840955095491 0.740840955095491 0.740840955095491x0 = 0.740840955095491 可见法迭代5次即趋于稳定.3.2 Matlab非线性方程求根的命令1. 代数方程组的求根rootsr=roots(P)P为代数方程的系数向量,从高次到低次排列,该命令只能求出一个一元多项式的根,返回所有复数和实数根.2. 求零点fzerox=fzero(F,x0,option)F为函数和字符表达式、内联函数或M文件的函数形式,x0为零点的大概位置,option为可选项,可参阅其他资料得到详细说明.使用该命令前,前配合plot、ezplot和fplot先画出函数图形，再猜测x0的大概值.3 求方程组数值解的命令fsolvex=fsolve(fun,x0,options)fun为M函数文件名,一般可用代替单引号进行标识.x0是向量或矩阵,它是探索方程的起点,输出为与x0同维的向量或矩阵,是方程组的数值解.对于更详细的设置option,可用help fsolve获得详细说明.4. 非线性方程组的解析解solvesolve(eqn1,eqn2,eqnN)对N个方程采用默认变量求解.solve(eqn1,eqn2,eqnN,var1,var2,.,varN)对N个方程的var1,var2,.,varN变量求解.S= solve(eqn1,eqn2,eqnN,var1,var2,.,varN)对N个方程的var1,var2,.,varN变量求解.x1,x2,xn= solve(eqn1,eqn2,eqnN,var1,var2,.,varN)将求解分别赋给变量x1,x2,xn.3.3 实验作业实验名称：解非线性方程的迭代法实验目的：掌握解非线性方程的迭代法；实验内容：求方程的根实验要求：用牛顿迭代，同样计算到为止输出迭代初值及各次迭代值和迭代次数49第4次 插值法拉格朗日插值：算法流程图4.1 例题 例4.1 已知，用线性插值法求的近似值解：Matlab中有直接进行线性插值计算的命令interp1,直接使用interp1命令即可 x=4 9; y=2 3; f=interp1(x,y,7,linear) % 选项使用线性插值f = 2.6000故插值计算结果为例4.2 ,用二次Lagrange插值求的近似值解：求解过程描述如下：format long; % 显示15位% 输入初始数据x0=4 9 16;y0=2 3 4;x=7; % 插值点n=length(x0); % 取长度% 初始计算s=0;% 进入公式计算for j=0:(n-1) t=1; for i=0:(n-1) if i=j t=t*(x-x0(i+1)/(x0(j+1)-x0(i+1); end end s=s+t*y0(j+1);ends % 显示输出结果format short;程序输出结果:s =例4.3 设有函数在-5,5上取等距节点上进行Lagrange插值解：采用n次进行插值,为使计算结果可视化,画图显示插值结果,编制如下程序：%输入初始数据clc;x0=-5:5;y0=1./(1+x0.2);x=-5:0.05:5;%插值计算点m=length(x);for k=1:m%分别对每一点进行插值 tx=x(k);%插值点 n=length(x0); s=0; %进入迭代计算过程 for j=0:(n-1) t=1; for i=0:(n-1) if i=j t=t*(tx-x0(i+1)/(x0(j+1)-x0(i+1); end end s=s+t*y0(j+1); end yx(k)=s;end%画图显示结果figure;plot(x,yx,:m,LineWidth,2);hold on;plot(x0,y0,-r,LineWidth,0.5);grid on;title(Lagrange插值时的Runge振荡现象);xlabel(x轴);ylabel(y轴);legend(插值曲线,原曲线);运行程序输出结果如下：由图可知确实在附近发生了振荡,误差很大.例4.4. 设有函数在-5,5上取等距节点上进行分段线性插值并且绘制图像。Matlab程序如下Function y=div_linear(x0, y0, x, n)for j=1: length(x);for i=1:n-1 if (x = x0(i) & (x = x0(i+1) y=(x-x0(i+1)/ (x0(i)-x0(i+1)*y0(i)+ (x-x0(i)/ (x0(i+1)-x0(i)*y0(i+1); else continue;end %end conditional statementendfunction y=f(x) %原来的函数y=5./(x.2+1)；x0=-5:1:5； %插值节点w=length(x0); %插值区间长度n=w-1; %节点个数for x=-5: 0.01:5 %加密了的插值节点 y=div_linear(x0, f(x0), x, n); %分段插值后的点的纵坐标hold on;plot(x, y, r); plot(x, f(x), b);end运行节点画图如下例4.5 设有函数在-5,5上取等距节点上进行三次样条插值，并且绘制图像。注解：我没有理解此程序。上实验课程的时候可以选择使用Matlab已经封装了的函数y=spline(x0,y0,x)直接进行计算并且画图。解：% X为1阶导数横坐标向量% Y为1阶导数纵坐标向量% dx0和dxn导数边界条件function S=csfit(X, Y, dx0, dxn) %3次样条插值函数？N=length(X)-1;H=diff(X);D=diff(Y)./H;A=H(2: N-1);%边界约束B=2*(H(1:N-1)+H(2: N)C= H(2: N);U=6*diff(D);B(1)=B(1)-H(1)/2;U(1)= U(1)-3*(D(1);B(N-1)= B(N-1)-H(N)/2;U(N-1)=U(N-1)-3*(-D(N);for k=2: N-1temp=A(k-1)/B(k-1);B(k)=B(k)-temp*C(k-1);U(k)=U(k)-temp*U(k-1);endM(N)=U(N-1)/B(N-1);for k=N-2: -1: 1M(k+1)=(U(k)-C(k)*M(k+2)/B(k);endM(1)=3*(D(1)-dx0)/H(1)-M(2)/2;M(N+1)=3*(dxn-D(N)/H(N)-M(N)/2;for k=0: N-1S(k+1, 1)=(M(k+2)- M(k+1)/6*H(k+1);S(k+1, 2) = M(k+1)/2;S(k+1, 3) = D(k+1)-H(k+1)*(2*M(k+1)+ M(k+2)/6;S(k+1, 4) = Y(k+1);endclear allclcclfx=-5:1:5;y=5./(x.2+1)；X=-5:1:5;Y=y;dx0=0.07396449704142;dxn=-0.07396449704142;S= csfit(X, Y, dx0, dxn);x1=-5:0.01:-4;y1=polyval(S(1, :), x1-X(1);x2=-4:0.01:-3;y2=polyval(S(2, :), x2-X(2);x3=-3:0.01:-2;y3=polyval(S(3, :), x3-X(3);x4=-2:0.01:-1;y4=polyval(S(4, :), x4-X(4);x5=-1:0.01:0;y5=polyval(S(5, :), x5-X(5)