MATLAB---ch4数值计算

[x,favl]=fzero(fun,x0)
>>xzero=fzero(‘ humps ‘ , 1.2) xzero=
1.2995 >>yzero=humps(xzero , 1.2) yzero=
3.5527e-15
% look for a zero near 1.2 % evaluate at xzero
➢sum（X） 按列求和 ➢cumsum （X） 按列求元素累计和 ➢trapz(x , y) 给定数据点x和y，计算y=f(x)下
的梯形面积积分。
➢cumtrapz(x , y) 给定数据点x和y，计算y=f(x)
下的梯形面积累计积分。
6
MATLAB提供最简单的积分函数是梯形法trapz(x,y)， 其中x,y分别代表数目相同 的阵列或矩阵，而y与x的关系可以由 是一函数型态（如y=sin(x)）或是不以函数描述的离散型态 。 我们看一简单积分式
2
而 的 数值微分则为 dy=diff(y)./diff(x)。
例：要计算下列多项式在 [-4, 5] 区间的微分
x=linspace(-4,5); % 产生100个x的离散点 p=[1 -3 -11 27 10 -24]; f=polyval(p,x); plot(x,f) % 将多项式函数绘图 title('Fifth-deg. equation') dfb=diff(f)./diff(x); % 注意要分别计算diff(f)和diff(x) xd=x(2:length(x)); % 注意只有99个df值，对应x2,x3,...,x100的点 plot(xd,dfb ) % 将多项式的微分项绘图 title('Derivative of fifth-deg. equation')
MATLAB没有专门提供求函数最大值的函数，但只要注意到 -f(x)在区间(a,b)上的最小值就是f(x)在(a,b)的最大值，
9
常微分方程的数值求解
基于龙格－库塔法，MATLAB提供了求常微分方程数值解的函数， 一般调用格式为： [t,y]=ode45(odefun,tspan,y0) 其中odefun是定义f(t,y)的函数文件名，该函数文件必须返回一个 列向量。tspan形式为[t0,tf]，表示求解区间。y0是初始状态列向 量。t和y分别给出时间向量和相应的状态向量。
>>A = pascal(3) ％产生3×3的 对称矩阵 A=
111 123 136 >>d = det(A) ％求行列式的值 d= 1 >>C = inv(A) ％求逆矩阵 C= 3 -3 1 -3 5 -2 1 -2 1
14
矩阵特征值
➢特征方程：Ax=λx
▪ 由|A-λI|=0 求特征值λi ▪ 由Ax=λ i x 求特征向量x i
22
X(any(isnan(X)'),:) = [];
➢删除越界的数据点:
load count.dat;%读数据 mu=mean(count);%计算平均值mu = 32.0000 46.5417 65.5833 sigma=std(count);%计算偏差sigma =25.3703 41.4057 68.0281
方程ax=b
1 2 x1 = 8
2 3 x 2 13 a x=b
a=[1 2;2 3];b=[8;13];
x=inv(a)*b x=
x=a\b x=
2.00
2.00
3.00
3.00 19
一般代数方程求解
正如人们对寻找函数的极点感兴趣一样，有时寻找函数 过零或等于其它常数的点也非常重要。一般试图用解析 的方法寻找这类点非常困难，而且很多时候是不可能的。 函数fzero寻找一维函数的零点。
婴幼儿体格生长
MATLAB---ch4数值计算
4.1 数值微积分
有限微分（差分）
与积分相比，数值微分非常困难。一个函数小的变化，容易产生相邻点 的斜率的大的改变。由于微分这个固有的困难，所以尽可能避免数值微分，特 别是对实验获得的数据进行微分。在这种情况下，最好用最小二乘曲线拟合这 种数据，然后对所得到的多项式进行微分。
5
数值求和与数值积分
求解定积分的数值方法多种多样，如简单的梯形法、辛普生 (Simpson)•法、牛顿－柯特斯(Newton-Cotes)法等都是经常采用 的方法。它们的基本思想都是将整个积分区间[a,b]分成n个子区 间[xi,xi+1]，i=1,2,…,n，其中x1=a，xn=b。这样求定积分问题 就分解为求和问题。
• m=n, 恰定方程组，可以精确求解 • m>n, 超定方程组，可以求得最小二乘解 • m<n, 欠定方程组，求得一个基解
17
恰定方程组的解
方程ax＝b(a为非奇异) x=a-1 b 矩阵逆
两种解:
➢x=inv(a)b — 采用求逆运算解方程 ➢x=a\b — 采用左除运算解方程
18
例: x1+2x2=8 2x1+3x2=13
各种概率分布的交互式观察界面 disttool
24
样本分布的频数直方图描述
hist指令的使用示例。 randn('state',1),rand('state',31)%初始化 x=randn(100,1);y=rand(100,1);%生成正态和均匀分布实验样本 %观察正态数据组的频数直方图在不同区间分段数时的变化 subplot(1,2,1),histfit(x,20) %20区间情况 subplot(1,2,2),hist(y,20) %20区间情况（带正态拟合线）
• a = [-2 -1 0 1 2] • isnan(0./a) • Warning: Divide by zero. • ans = • 00100
• 删除NaN
• i = find(~isnan(x));x = x(i) • x = x(find(~isnan(x))) • x = x(~isnan(x)); • x(isnan(x)) = []; • X(any(isnan(X)‘),:) = []; ％删除任意含有NaN的行 • function X = excise(X)
[n,p]=size(count);
outliers=abs(count-mu(ones(n,1),:))>3*sigma(ones(n,1),:)
%判断偏差过大的数据所在位置 nout=sum(outliers); %计算越界数据点数 nout = count(any(outliers‘),:)=[];%删除越界数据
[x,fval,exitflag,output]=fmin(fun,x1,x2,options) [x,fval,exitflag,output]==fmins(fun,x0,options)
这两个函数的调用格式相似。 ✓其中fmin函数用于求单变量函数的最小值点。fun是被最小化 的目标函数名，x1和x2限定自变量的取值范围。 ✓fmins函数用于求多变量函数的最小值点，x0是求解的初始值 向量。
3
➢其他用法：
• DX=diff(X)：计算向量X的向前差分，DX(i)=X(i+1)X(i)，i=1,2,…,n-1。
• DX=diff(X,n)：计算X的n阶向前差分。例如， diff(X,2)=diff(diff(X))。
• DX=diff(A,n,dim)：计算矩阵A的n阶差分，dim=1时 (缺省状态)，按列计算差分；dim=2，按行计算差分。
V=
0.7212 0.4443 0.5315
-0.6863 0.5621 0.4615
-0.0937 -0.6976 0.7103
D=
-0.0166 0
0
0 1.4801 0
0
0 2.5365
16
求解线性方程组
➢矩阵方程AX=b的解
• X=A\b
➢矩阵方程XA=b的解
• X=A/b
➢三种方程组情况的求解（对m×n矩阵A）
以下为 MATLAB 的程式
>> x=0:pi/100:pi; >> y=sin(x); >> k=trapz(x,y) k= 1.9998
7
精度可控的数值积分
S1=quad（fun,a,b,tol) 采用递推自适应Simpson法计算积分。 S1=quadl（fun,a,b,tol) 采用递推自适应Lobatto法求数值积分。 S2=dblqual（fun,xmin,xmax,ymin,ymax,tol) 二重(闭型)数值 积分指令。 S3=triplequal （fun,xmin,xmax,ymin,ymax,zmin,zmax,tol)三重 (闭型)数值积分指令。
10
单个高阶常微分方程处理方法
11
4.2 矩阵与线性代数-矩阵运算
➢加(+)、减(-) ➢乘：(*）
• 矩阵之间的乘、向量与矩阵相乘、标量与矩阵相乘
➢除：右除（/）、左除（\）
• 数学上没有矩阵除法的定义。
➢幂 ^ ➢转置 ‘ ➢其他：
• 指数运算 Y = expm(A) 还有:expm1、expm2、
2．矩阵的迹
>> trace(A)
矩阵的迹等于矩阵的对 ans =
角线元素之和，也等于 5
矩阵的特征值之和。在
MATLAB中，求矩阵的
迹的函数是trace(A)。 13
行列式和逆
➢行列式：d = det(A)
• 计算方阵的行列式值
➢矩阵的逆：C= inv(A)
• AC=CA=I ，则：C=A-1 • 其中A是非奇异方阵 • 实际数值计算中较少使用
➢eig求解Ax=λx特征值问题
• d = eig(A) 只返回特征值向量d
• [V,D] = eig(A) 返回特征向量阵V和特征值对角阵 D
• [V,D] = eig(A,‘nobalance’)矩阵中有与截断误差 （eps)相当的元素时,计算精度更高。
• [V,D] = eig(A,B) 计算广义特征向量阵V和广义特 征值阵D，使AV=BVD成立。
例: 求定积分。
(1) 建立被积函数文件fesin.m。
function f=fesin(x)
f=exp(-0.5*x).*sin(x+pi/6);
(2) 调用数值积分函数quad求定积分。
S=quad('fesin',0,3*pi)
S=
0.9008
8
函数极值的数值求解
MATLAB提供了基于单纯形算法求解函数极值的函数 fminbnd和fminsearch， 它们分别用于求单变量函数和多变量函 数的最小值，其调用格式为：
在MATLAB中，没有直接提供求数值导数的函数，只有计算向前
差分的函数diff，其调用格式为：Y = diff(X)
计算向量中连续两个元素的差Y=[X(2)-X(1) X(3)-X(2) ... X(n)-X(n-1)]
例如：A = [9 -2 3 0 1 5 4];
diff(A)
ans =
-11 5 -3 1 4 -1
➢几个应用：
• diff(x)==0 相邻元素是否相等 • all(diff(x)>0) 是否单调递增 • all(diff(diff(x))==0) 是否均匀分布
4
梯度
➢FX = gradient（F） 求一元（函数）梯度 ➢[FX,FY]= gradient（F） 求二元（函数）梯度
计算梯度，绘制引力线图：
100
23
概率分布函数
二项分布(Binomial distribution) Pk=binopdf(k,N,p) Fk=binocdf(k,N,p) R=binordf(N,p,m,n)
正态分布（Normal distribution） Px=normpdf(x,Mu,Sigma) Fx=normcdf(x,Mu,Sigma) R=normrnd(Mu,Sigma,m,n)
expm3 • 开方 Y = sqrtm(A) • 对数 Y = logm(A)
12
矩阵的秩与迹
1．矩阵的秩 矩阵线性无关的行数与 列数称为矩阵的秩。在 MATLAB中，求矩阵秩 的函数是rank(A)。
>> A=[2 1; 4 3];
>> rank(A)
ans =
2 % 表示A秩数为2且 等于矩阵的列数
20
4.3 概率分布和统计数据分析 —数据集
➢单变量的统计数据存储在1×n或
n×1的向量里。
➢多变量的数据用矩阵来表示，不同
的变量放在不同的列。
21
数据预处理
➢MATLAB用NaN表示非数 ➢用NaN处理缺少的数据为数据分析带来方便
• TF = isnaБайду номын сангаас(A) 返回与A同维的矩阵TF，A中为NaN的得到 逻辑1，否则为逻辑0
15
【例】求下列对称矩阵A的特征值/特征向量
>>A =[1.0000 1.0000 0.5000
1.0000 1.0000 0.2500
0.5000 0.2500 2.0000 ]
>>eig(A)%返回特征值 ans = -0.0166
1.4801 2.5365
>> [V,D] = eig(A) %返回特征值及特征向量