第三章数值计算3.1LU分解和恰定方程组的解Matlab矩阵的分解形式主要有：三角分解、正交化、特征值分解。
3.1.1LU三角分解、行列式和逆1、LU分解是利用Gaussian列主元消去法进行的。
为了保证主元消去策略的实施，一般来说，必须对被分解矩阵实施行置换因此有：LU=PA 式中L为主对角元为1的下三角矩阵，U是上三角矩阵，P是由0或1组成的行置换矩阵：［L UP］=lu(A)2、A-1=U-1L-1P Matlab根据以上算法编制了相应的指令：det(a) 求矩阵a的行列式inv(a) 求矩阵a的逆矩阵3.1.2恰定方程组的解在求解方程式，尽量不要使用inv(a)*b指令，而应采用a\b，后者计算速度比前者快，精确度高。
【例3.1.2-1】“求逆”法和“左除”法解恰定方程的性能对比（1）randn('state',0);A=gallery('randsvd',100,2e13,2); ％产生条件数为2e13的100阶随机矩阵x=ones(100,1);b=A*x;（2）ticxi=inv(A)*b;ti=toc（3）tic;xd=A\b;td=toc;3.1.3矩阵特征值和矩阵函数矩阵A与向量x相乘,即表示矩阵对向量的变换(transformation)，一般说来，向量在变换的作用下将发生旋转（rotation），反射(reflection)和放大缩小。
但对于任何一个矩阵来说，中存在那么一些特殊的向量，再对其变化的作用下，向量的方向不变，而仅长短发生变化。
这种向量就是所谓的特征向量（eigenvector），它满足方程 Ax=λx3.1.3.1特征值和特征向量的求取d=eig(a) 仅仅计算a的特征值[v d]=eig(a) 计算矩阵a的特征向量阵v和特征对角阵d，使av=vd成立。
【例3.1.3.1-1】简单实阵的特征值问题。
A=[1,-3;2,2/3];[V,D]=eig(A)V =0.7746 0.77460.0430 - 0.6310i 0.0430 + 0.6310iD =0.8333 + 2.4438i 00 0.8333 - 2.4438i3.2数据分析3.2.1 基本统计函数指令note: 1、median(x) 当x为向量时，先把x元素有小到大排列，然后在新排成数组中；n为奇取（n+1）/2个元素；n为偶数时取n/2与(n/2+1)元素的平均值，作为总位数。
2、输入为向量，则运算对整个向量进行的。
若输入是数组，那么指令运算时按列进行的。
结果为一行向量。
3、Max(a,b)取a b矩阵中加大的数.clear all; a=reshape(1:9,3,3);b=reshape(1:16,4,4);median(a),median(b) ans =2 5 8ans =2.5000 6.5000 10.5000 14.5000clear all;a=reshape(1:16,4,4);b=max(a),c=max(max(a)),b =4 8 12 16c =16[x y]=max(a)x =4 8 12 16y =4 4 4 4【例】hist(histogram)指令的使用示例。
randn('state',1),rand('state',31)x=randn(1000,1);y=rand(1000,1);%图 3.2-13.3函数的数值导数(1) 数值差分和导数、偏导数dx=diff(x) ％求X相邻行元素间的一阶差分。
dx=diff(x,n) ％求X相邻元素间的n阶差分。
dx=diff(x,n,dim) ％在dim指定的维上，求X相邻元素间的n阶差分。
note：1、diff是基于前向差分概念设计的，即dx/dt=(x(t+h)-x(t))/h2、数值导数的求取应尽量避免。
3、一般需先通过多项式拟合，或通过样条拟合，然后再从元数据拟合函数求导。
4、利用多项式求导。
【例3.3-1】y=sin(x)的导数为cos(x)，利用数值求导及多项式拟合求导，并进行比较。
x=0:pi/10:2*pi; y=sin(x);plot(x,y);hold on;y1=diff(y);plot(x(2:end),y1);plot(x,cos(x),':r');p=polyfit(x,y,5);p1=polyder(p);y2=polyval(p1,x);p（2）数值梯度[fx fy]=gradient(F,h)note:1、数值梯度使用的场合与数值导数相同2、MATLAB约定：数组的“行”数据点沿x轴取得，“列”数据点沿y轴取的。
3、h为步长。
（3）方向导数的可视化quiver(x,y,u,v,scale) 在（x,y）二维平面点上，画（u,v）方向的箭头。
scale 表示箭头的长度，默认值为1。
【例3.3-2】用一个简单矩阵表现diff和gradient指令计算方式。
F=[1,2,3;4,5,6;7,8,9]Dx=diff(F)Dx_2=diff(F,1,2)[FX,FY]=gradient(F)[FX_2,FY_2]=gradient(F,0.5)【例 3.3-3】研究偶极子的电势和电场强度。
设在),(b a 处有电荷q +，在),(b a --处有电荷q -。
那么在电荷所在平面上任何一点的电势和场强分别为)11(4),(0-+-=r r q y x V πε，V E -∇= 。
其中2222)()(,)()(b y a x r b y a x r +++=-+-=-+。
9010941⋅==πεk 。
又设电荷6102-⋅=q ，5.1=a ，5.1-=b 。
clear;clf;q=2e-6;k=9e9;a=1.5;b=-1.5;x=-6:0.6:6;y=x;[X,Y]=meshgrid(x,y);rp=sqrt((X-a).^2+(Y-b).^2);rn=sqrt((X+a).^2+(Y+b).^2); V=q*k*(1./rp-1./rn); [ex,ey]=gradient(-V);ae=sqrt(ex.^2+ey.^2);ex=ex./ae;ey=ey./ae; cv=linspace(min(min(V)),max(max(V)),49); contourf(X,Y,V,cv,'k-') axis('square')title('\fontname{隶书}\fontsize{22}偶极子的场'),hold on quiver(X,Y,ex,ey,0.7) plot(a,b,'wo',a,b,'w+')plot(-a,-b,'wo',-a,-b,'w-')xlabel('x');ylabel('y'),hold off3.4 函数的零点和极点很多情况下，需要求函数数组子矩阵的零点和极点3.4.1 多项式的根出于计算考虑，Matlab 不对多项式直接求根，而是通过求他的伴随矩阵的特征值进行，对于阶数不超过20的多项式，这种处理方法被认为是最适当可靠的计算方法。
指令：roots(p) 求多项式P 的根，其中p 为多项式系数的行向量。
【例 3.4.1-1】利用多项式计算y=sin(x)，在零到2＊pi 区间内5阶多项式拟合，并计算此多项式的的零点。
clear all;x=0:pi/10:2*pi;y=sin(x); p=polyfit(x,y,5); roots(p)clear all;p=[1 -5 6]; poly2str(p,'x'),roots(p)一元函数的零点与多项式不同，任意函数f(x)=0可能有零点，也可能没有零点，可能有一个零点，也可能有多个零点，因此，很难说出一个通用解法。
一般来说，零点的数值计算过程是，先猜测一个初始零点或该令点所在的区间，然后通过一些计算，是猜测值不断精确化，或是猜测区间不断收缩，直到达到预先指定的精度，终止计算。
指令: z=fzero(fun,xo)note:1、该函数只能求去一元连续函数穿越横轴的零点，而不会确定连续函数曲线，那种子触及横轴而不穿越横轴的零点。
例如abs(sin(x))。
2、fun 为字符串或内联函数（inline function ）4、x0时表示零点初始猜测，此点应该靠近零点5、无法判断x0时，可利用plot指令。
Example 1Calculate π by finding the zero of the sine function near 3.x = fzero(@sin,3)x =3.1416Example 2To find the zero of cosine between 1 and 2x = fzero(@cos,[1 2])x =1.5708Note that cos(1) and cos(2) differ in sign.Example 3To find a zero of the function f(x) = x3–2x–5, write an anonymous function f:f = @(x)x.^3-2*x-5;Then find the zero near 2:z = fzero(f,2)z =2.0946Because this function is a polynomial, the statement roots([1 0 -2 -5]) finds the same real zero, and a complex conjugate pair of zeros.2.0946-1.0473 + 1.1359i-1.0473 - 1.1359iIf fun is parameterized, you can use anonymous functions to capture the problem-dependent parameters. For example, suppose you want to minimize the objective function myfun defined by the following function file:function f = myfun(x,a)f = cos(a*x);Example 4求sin(x)靠近pi的零点z=fzero('sin(x)',3)Example 5求sin(x)/x的主瓣上的零点clear all;x=-2*pi:pi/10:2*pi;y=sin(x+eps)./(x+eps);figuregrid ;holdplot(x,y);hold off[a b]=ginput(2);fs=inline('sin(x)/x','x');z1=fzero(fs,a(1));z2=fzero(fs,a(2));3.4.2函数极值点许多科学研究和工程计算问题口可以归结为一个极值问题,如能量最小,时间最短,最佳拟合等等。