第二讲-matlab数学建模之数值计算.ppt

* A(:) 将矩阵A的每一列元素叠加起来,成为一个列向量. * prod(prod(A,2)) prod(prod(A,1)) %求A的全部元素的乘积(or prod(A(:))) * 标准方差就是方差的平方根：一组数据中的每一个数与这组数据的平均数的差的平方的和再除以数据的个数，取平方根既是。 S1 无偏 * 线性,最近 * * p=polyder(P)：求多项式P的导函数 P是多项式的向量表示. sinx * diff(X)/h 导函数(h为步长) * 差商近似代替sinx导数 在0处的值约等于0.9207(cos0=1, 步长h取的不好, 太大) * Trapz 梯形公式 * Trapz 梯形公式 Lobatto洛巴托 trace不起作用 * f=inline(‘exp(-x.^2/2).*sin(x.^2+y)’,’x’,’y’) 其中’x’,’y’可删除, 结果一样. * % or b=[13;-9;6;0]; * (Hessenberg)人名；(德)黑森贝格 * error ??? Error using == error Not enough input arguments. Y相当于x(k+1) x0相当于x(k) * Y相当于x(k+1) x0相当于x(k) * 满秩方阵: n阶矩阵非奇异 * Compute the rational basis for the null space of the same matrix A. ZR = null(A,r) * eps ans = 2.220446049250313e-016 * F(-4)0 f(-6)0 * F(-4)0 f(-6)0 * 注意 积分是点运算, 非线性方程(组)不用点运算也可 * Fx书上不一致 * Fx书上不一致 * Inline 及‘’ 不好用 * t和y分别给出时间向量和相应的状态向量。 * t‘ y‘ y1 * 2．Gauss-Seidel迭代法 Jacobi迭代公式Dx(k+1)=(L+U)x(k)+b可以改进为Dx(k+1)=Lx(k+1)+Ux(k)+b，于是得到： x(k+1)=(D-L)-1Ux(k)+(D-L)-1b 该式即为Gauss-Seidel迭代公式。 和Jacobi迭代相比，Gauss-Seidel迭代用新分量代替旧分量，精度会高些。 * Gauss-Seidel迭代法的MATLAB函数文件gauseidel.m如下： function [y,n]=gauseidel(A,b,x0,eps) if nargin==3 eps=1.0e-6; elseif nargin3 error return end D=diag(diag(A)); %求A的对角矩阵 L=-tril(A,-1); %求A的下三角阵 U=-triu(A,1); %求A的上三角阵 G=(D-L)\U; f=(D-L)\b; y=G*x0+f; n=1; %迭代次数 while norm(y-x0)=eps x0=y; y=G*x0+f; n=n+1; end * 例 用Gauss-Seidel迭代法求解下列线性方程组。设迭代初值为0，迭代精度为10-6。 在命令中调用函数文件gauseidel.m，命令如下： A=[10,-1,0;-1,10,-2;0,-2,10]; b=[9,7,6]; [x,n]=gauseidel(A,b,[0,0,0],1.0e-6) % or [x,n]=gauseidel(A,b,[0,0,0]) % or x=gauseidel(A,b,[0,0,0]) x = 0.9958 0.9579 0.7916 n = 7 * 例 分别用Jacobi迭代和Gauss-Seidel迭代法求解下列线性方程组，看是否收敛。 命令如下： a=[1,2,-2;1,1,1;2,2,1]; b=[9;7;6]; [x,n]=jacobi(a,b,[0;0;0]) [x,n]=gauseidel(a,b,[0;0;0]) x = -27 26 8 n = 4 x = NaN NaN NaN