MATLAB语言及应用--第6章数值计算

d b( s ) (3)[q,d]=polyder(b,a)可以求出多项式之比b(s)/a(s)的导数。 [q( s) d ( s)] [ ] 导函数的分子存入q，分母存入d ds a( s)
15
6.3 多项式和卷积
7） 多项式曲线拟合
p=polyfit(x,y,n);
测量数据对 多项式的阶数
内容提要
1、矩阵的相关计算 2、线性方程组的解 3、多项式和卷积 4、拟合与插值 5、matlab的泛函指令
2
6.1 矩阵的相关计算
det(A) inv(A)
求矩阵A的行列式
A=[1,-1,1;5,-4,3;2,1,1];
求矩阵A的逆
B=inv(A); A*B,B*A
如果矩阵A不是方阵，或A是非满秩方阵时，矩 阵A没有逆矩阵，但可以找到一个与A的转置矩 阵A’同型的矩阵B，使得 A*B*A=A, B*A*B=B
有理多项式有两种表示方法： [r,p,k]=residue(b,a);
rn r1 r2 H ( s) ... k (s) s p1 s p2 s pn
b1 b2 s 1 ... bm 1s m H (s) a1 a2 s 1 ... an 1s n
[t,y]=ode45(‘fun’ ,tspan,y0)
28
6.5.1 数值微分
MATLAB中没有直接提供求数值导数的函数， 只有计算向前差分的函数diff DX=diff(X) 计算向量X的向前差分，DX(i)=X(i+1)-X(i) DX=diff(X,n) 计算X的n阶差分，如diff(X,2)=diff(diff(X)) DX=diff(A,n,dim) 计算矩阵A的n阶差分。 dim=1(默认状态）按列计算差分； dim=2按行计算差分。
21
6.4 拟合与插值
插值： 是在认定所给“基准数据”完全正确的情况下，
研 究如何“平滑”的估算出“基准数据”之间其它 %给定数据对 点的函 x0=0:0.1:1; 数值。因此，插值所得曲线一定穿过“基准数 y0=[-.44,1.97,3.11,5.25,5.02,4.66,4.01,4.58,3.45,5.35,9.22]; 据”。 %采用三次多项式进行插值 xi=0:0.02:1; yi=interp1(x0,y0,xi,'cubic'); %绘图 plot(xi,yi,'-b',x0,y0,'.r','MarkerSize',20),xlabel('x')
29
例：设 x由[0,2*pi]间均匀分布的10个点组成，
求sinx的1~3阶差分。 x=linspace(0,2*pi,10);
Y=sin(x);
DY=diff(Y); D2Y=diff(Y,2); D3Y=diff(Y,3);
30
6.5.2
数值积分
例：求以下的定积分。其精确值为0.7468204…
I e
0
1
x2
dx
%（1）使用字符串表示被处理函数 fun='exp(-x.*x)';
%（2）调用积分指令求积分值 E1 = quad (fun,0,1) E2 = quadl (fun,0,1) 区别：求数值积分是采取不同的方法
31
6.5.3
常微分方程求解
常微分方程数值解法的步骤
1、常微分方程化为一阶微分方程组
20
6.4 拟合与插值
p=polyfit(x,y,n);%求得最小二乘拟合多项式的系数
%给定数据对
x0=0:0.1:1; y0=[-.44,1.97,3.11,5.25,5.02,4.66,4.01,4.58,3.45,5.35,9.22];
%求拟合多项式系数 n=3; P=polyfit(x0,y0,n) %图示拟合情况 xx=0:0.01:1; yy=polyval(P,xx); 拟合：寻找一条“平滑”曲线来最好的表现带噪声的“测量数据”。 plot(xx,yy,'-b',x0,y0,'.r','MarkerSize',20),xlabel('x') 但并不要求拟合曲线穿过这些“测量数据”点。
1
3];
y=2*X*X+3*X +eye(3)
12
6.3 多项式和卷积
4）多项式的根
功能：计算多 项式P的根。
R=roots(P)
p1 x p2 x
n
n 1
... pn x pn 1 0
-27 ];
P= [ 1 -6 -72 R=roots(P)
13
6.3 多项式和卷积
5）有理多项式的部分分式展开：
23
6.4 拟合与插值 1）最邻近插值方法（nearest）
插值点的值与其最邻近的点的函数值相等。
*
*
*
*
*
24
6.4 拟合与插值 2）线性插值方法（ linear ）
插值点的值在前后两个数据点所构成 的直线上。
*
*
*
25
*
*
6.4 拟合与插值 3）三次样条插值方法(spline)
利用一系列样条函数获得内插数据点， 从而确定已有数据点之间的函数。
y
(n)
f (t , y , y ,..., y ( n 1) )


y1 y, y2 y,... yn y ( n1)
y 1 y2

令
y 2 y3


yn f (t , y1 , y2 ,... yn )
32
常微分方程数值解法的步骤
确定初值
y1 (0), y2 (0),..., yn (0)
27
6.5 matlab的泛函指令
在matlab中，凡是以“函数”为输入宗量的指令， 都被统称为matlab的泛函指令。 z=fzero ( fun , x0 主要内容： 被处理“函数”fun可以有三种形式： )
1、数值微分 （1）字符串表达式 q=quad(fun,a,b) 2、数值积分 （2）内联函数 （语句函数） （3）“M函数文件”的函数句柄 3、解微分方程 q=quadl(fun,a,b)
y=polyval(p,x); y=polyvalm(p,x);
功能：按矩阵运算规 则计算多项式的值。 其中x必须为方阵。
11
例题开讲
P=[ 2 3 y=2*X.*X 1]; +3*X+1 X=[ 1 2 3 ; 3 2 1 ; 2
y=polyval ( P , X ) y=polyvalm ( P , X )
[b,a]=residue(r,p,k);
14
6.3 多项式和卷积
6） 多项式求导 q=polyder(p);
具体用法:
d q( s ) [ p( s)] ds
(1)q=polyder(p)可以得到多项式p的导数。 d q( s ) [a( s)b( s)] ds (2) q=polyder(a,b)可以求出多项式a,b之积的导数。
4）三次曲线内插方法（cubic）
构造三次曲线函数来拟合已知数据x、y，从 而确定内插点的值。
26
说明: 1、四种插值方法中，x中的数据是单调但不一定是均匀的。 6.4 拟合与插值 2、若已知x为均匀的，则在method前加*，可使执行速度加快。 3、按nearest linear cubic spline的顺序，对内存要求 从小到大，执行速度由快到慢，平滑度由差到好。
17
6.3 多项式和卷积
（2）解卷积：
[q,r]=deconv(c,b);
c( s ) r (s) q( s) b( s ) b( s )
例如： b=[4,5,6]; c=[4,13,28,27,18];
18
6.4 拟合与插值
19
6.4 拟合与插值
在matlab中实现最小二乘法拟合通常可以采用如下两种途径： （1）利用polyfit函数进行多项式拟合。 （2）利用矩阵除法解决复杂型函数的拟合。
A*V=V*D
4
6.2 线性方程组的解
关于线性方程组Ax=b的解法一般可以分为两类： （1）直接法：通过矩阵的变形、消去直接得到方 程的解，适用于解低阶稠密矩阵（阶数大约 <150）方程组。 （2）迭代法：是用某种极限过程去逐渐逼近方程 组精确解的方法，适用于解大型稀疏矩阵方 程组。
5
6.2
线性方程组的解
此时B为A的伪逆。pinv( )
A=[3,1,1,1;1,3,1,1;1,1,3,1];B=pinv(A)
3
6.1 矩阵的相关计算
Matlab的相关指令： rank(A) lamda=eig(A)或 [V,D]=eig(A)
求矩阵的秩 计算A(必须是方阵)的特征值，以向量形式存放 计算A的特征向量阵V和特征值对角阵D
在matlab命令窗中的指令是： a=[0.4096 … 0.2943;…;0.1784 … 0.3927]; b=[0.4043 0.1550 0.4240 -0.2557]’;
X=a\b
6
线性方程组的解
情形一：如果A是n*n的矩阵，b是n维向量， 则表示n个线性方程求解n的未知数。A可逆时， 该线性方程组有唯一解。 情形二：如果A是m*n阵（m<n），方程组个数 小于未知数个数，则该方程组有无穷多组解， maltab可以求出无穷多组解中的一组解。
22
6.4 拟合与插值
说明：如果被插值 函数是一个单变量 函数，则数据插值 称为一维插值， interp1
ys=interp1(x,y,xs,’method’); 其中:(1) x,y是测量数据对； (2) xs是需要内插的点所构成的向量。 (3) method是指所使用的内插方法。 插值算法： ‘nearest’，‘linear’，‘spline’，‘cubic’