matlab第五课.docx
《matlab第五课.docx》由会员分享，可在线阅读，更多相关《matlab第五课.docx（33页珍藏版）》请在冰点文库上搜索。
matlab第五课
Matlab_5数值方法
多项式
代数方程（优化）
数值微积分
数据统计、插值和拟合
常微分方程
偏微分方程
多项式
Matlab提供了一组函数用于处理多项式运算。
（不适合处理大于10的高阶多项式）
常用的函数：
roots,polyval,polyfit,…
%%多项式表达：
%用一个行向量表示各阶系数，阶次降序排列。
p1=[1-32]%x2-3x+2
p2=[10-23]%x3-2x+3
%%多项式求值：
polyval（p1,2）
polyval（p2,1）
%绘制多项式x3-2x+3的图形
x=-2:
0.1:
2;
y=polyval（p2,x）;
plot（x,y,'.-'）
%%多项式求根：
roots（p1）
roots（p2）
roots（[1100]）
%%构建多项式：
p=1:
5
r=roots（p）
pp=poly（r）
pp-p
%由于截断误差，函数poly生成的系数有微小的偏差
%有时候结果出现复数，可以用real函数提取实部，消除虚部的影响。
%多项式加减运算：
多项式加减就是系数向量的加减。
（Matlab没有提供多项式加减的函数）
%多项式乘法：
（系数向量的卷积运算）
a=[1357]
p=conv（a,1:
3）
%多项式除法：
[q,r]=deconv（p,a）
%q商,r余数
[q,r]=deconv（2:
6,a）
conv（a,q）+r
%多项式微分：
a=1:
3,b=[11111]
polyder（a）
polyder（b）
polyder（a,b）%多项式乘积的导数。
polyder（conv（a,b））
%有理多项式的导数，返回分子和分母多项式。
%注意Matlab函数的输入输出参数
[bb,aa]=polyder（b,a）
[bb,aa]=polyder（[10],[1-1]）
%多项式积分：
p=polyint（[111]）
polyder（p）
polyint（[111],10）%指定积分常数为10
多项式拟合
%调用格式：
[p,s,mu]=polyfit（x,y,n）
%%例：
x=0:
0.1:
1;
y=[-0.4471.9783.286.167.087.347.669.569.489.311.2];
plot（x,y,'or'）
n=2;
p=polyfit（x,y,n）%二阶拟合
x2=linspace（0,1）;
y2=polyval（p,x2）;
holdon,plot（x2,y2,'.-'）,holdoff
p10=polyfit（x,y,10）;%十阶拟合
p10.'
y10=polyval（p10,x2）;
holdon,plot（x2,y10,'-m'）,holdoff
%%查看拟合效果
[p,s]=polyfit（x,y,2）
yp=polyval（p,x）;
sqrt（sum（（yp-y）.^2））%s.normr
R2=1-sum（（yp-y）.^2）/sum（y.^2）
%%尺度变换的调用格式
[p,s,mu]=polyfit（x,y,2）
mean（x）%mu
（1）
std（x）%mu
（2）
xp=（x-mu
（1））/mu
（2）;
yp=polyval（p,xp）;
plot（xp,y,'o',xp,yp,'.-'）;
优化
函数y=f（x）在x取何值时，能够得到特定的y值或极值，Matlab称为优化。
对于实际问题直接利用反函数f-1（y）确定x是很困难的，通常需要迭代求解。
这个迭代过程实际上是求解y–f（x）=0的过程。
%humps函数
x=linspace（-0.5,1.5）;
y=humps（x）;
plot（x,y,'.-'）;
gridon;
一元函数寻零---fzero
%指定点附近搜索
fzero（@humps,1）%在x=1附近开始搜索
[x,r,exitflag,output]=fzero（@humps,-1）%返回值r为偏差
%指定搜索区域
fzero（@humps,[12]）%区域两段函数值符号相反
fzero（@humps,[01]）
fzero（@humps,[-33]）
%待处理的函数需在搜索区域连续
fzero（@tan,0.1）
fzero（@tan,[0.1,3]）
[x,r]=fzero（@tan,[0.1,3]）
求解选项的设置和获取---optimset和optimget函数
helpoptimset%查看帮助文件
%
options=optimset（'Display','iter','TolX',1e-3）
[x,value]=fzero（@humps,[-12],options）
一维最小值---fminbnd
调用格式：
[xmin,value,exitflag,output]=fminbnd（fun,x1,x2,options）
%指定在[x1,x2]区间内搜索最小值。
%在x=0.5~0.8区域搜索humps函数最小值
options=optimset（'Display','iter'）;
[xmin,value,flag,out]=fminbnd（@humps,0.5,0.8,options）
%恰当选择搜索区域
[xmin,value,flag]=fminbnd（@humps,0.5,1.4,options）
%最大值问题
%等同于寻找函数–f（x）的最小值。
h=@（x）-1*humps（x）;
[xmax,value,flag]=fminbnd（h,0.5,1.4,options）
多维最小值---fminsearch函数
%Rosenbrock函数和图形
[x,y]=meshgrid（-1.5:
0.1:
1.5,0:
0.1:
2.0）;
z=100*（y-x.^2）.^2+（1-x）.^2;
surf（x,y,z）;
%函数fminsearch输入函数的自变量写成向量形式
%对于较复杂的待求解函数，最好写成M函数
h=@（x）（100*（x
（2）-x
（1）.^2）.^2+（1-x
（1））.^2）;
[xymin,value,flag,output]=fminsearch（h,[-1,2]）%在（-1,2）附近搜索
求解非线性代数方程组---fsolve
docfsolve%查看函数fsolve用法
%调用格式：
%[x,fval,exitflag,output,jacobian]=fsolve（fun,x0,options）
%待求解问题必须写成标准形式F（x）=0
%例：
求解方程组
2x1–x2=e-x1
-x1+2x2=e-x2
%---------------------------------------------------------
functionF=funfsolve（x）
F=[2*x
（1）-x
（2）-exp（-x
（1））
-x
（1）+2*x
（2）-exp（-x
（2））];
%可以写成匿名函数的形式：
%h=@（x）[2-1;-12]*x-[exp（-x
（1））;exp（-x
（2））]
%---------------------------------------------------------
x0=[33];%testforx0=[00]&[-100100]
options=optimset（'Display','iter'）;
%[x,value]=fsolve（@funfsolve,x0,options）
[x,value,exitflag,output,jacobian]=fsolve（@funfsolve,x0,options）
%例：
求解满足X*X*X=[12;34]的矩阵X
h=@（x）x*x*x-[12;34];
[x,value,flag]=fsolve（h,[11;11]）
%例：
求解方程x3–2x2–x+2=0
%
h=@（x）x.^3-2*x.^2-x+2;
x=-2:
0.1:
3;
plot（x,h（x）,'-r'）;
gridon;
r1=fsolve（h,0）
r2=fzero（h,0）
r3=roots（[1-2-12]）
函数fzero/fminbnd/fminsearch/fsolve…小结：
1）在待求解域内，函数是连续的。
求解时总是假设问题有解，但有些问题不收敛或收敛速度极慢，计算会中断。
2）只能得到一个解。
如问题有多解，只能尝试改变不同的设定初值。
3）恰当设定初值，否则问题可能无解。
4）对多项式、线性方程组等问题不要用此类求解器。
5）更多的优化函数，参阅OptimizationToolbox。
6）可以直接使用交互式优化工具optimtool。
数值微积分
%%数值积分
%函数trapz按给定的间隔抽样、梯形分割近似计算。
%trapz可以直接处理离散数据。
%函数trapz的计算方法
x=[0134];
y=[0112];
plot（x,y,'-or'）
trapz（x,y）
%计算humps函数在[-1,2]区间的积分值
x=linspace（-1,2,200）;y=humps（x）;plot（x,y,'.-'）;
formatlong
x=linspace（-1,2,20）;%更多的抽样点n=20,50…
y=humps（x）;
area=trapz（x,y）
%函数cumtrapz计算
x=linspace（-1,2）;
y=humps（x）;
z=cumtrapz（x,y）;
plotyy（x,y,x,z）
%函数quad和quadl采用任意间隔抽样和高阶估计。
formatlong
area=quad（@humps,-1,2）
area=quadl（@humps,-1,2）
%例：
空间曲线的长度
%x（t）=sin（2t）,y（t）=cos（t）,z=t;t=[0,3*pi]
t=0:
0.1:
3*pi;
plot3（sin（2*t）,cos（t）,t）
h=@（t）sqrt（4*cos（2*t）.^2+sin（t）.^2+1）;
curvelength=quadl（h,0,3*pi）
广义积分---quadgk
docquadgk%查看几种积分函数的比较说明
%无界点
formatlong
quadgk（@（x）log（x）,0,1）
quadgk（@（x）abs（x）.^-0.5,-1,1）
%无穷区间
quadgk（@（x）exp（-x）,0,inf）
quadgk（@（x）1./（1+x.^2）,-inf,inf）
%（特殊点）
quadgk（@（x）sin（x）./x,-1,1）%2*sinint
（1）
quadgk（@（x）sin（1./x）,0,1）%sin
（1）-cosint
（1）
多重积分---dblquad/triquad/quad2d
%计算函数sin（x）cos（y）+1在[0,pi]&[-pi,pi]区域的二重积分
h=@（x,y）sin（x）.*cos（y）+1;
[x,y]=meshgrid（0:
pi/20:
pi,-pi:
pi/20:
pi）;
z=h（x,y）;
surf（x,y,z）
volume=dblquad（h,0,pi,-pi,pi）
%精确解
exact_solution=2*pi^2;
err=abs（volume-exact_solution）/exact_solution
%函数triplequad计算三重积分。
%如二维积分区域不是矩形可以用函数quad2d
%计算函数[（x+y）1/2（1+x+y）2]-1在[0,1]&[0,1-x]的二重积分值
[x,y]=meshgrid（0:
0.1:
1）;
z=1./（sqrt（x+y）.*（1+x+y）.^2）;
surf（x,y,z）;
h=@（x,y）1./（sqrt（x+y）.*（1+x+y）.^2）;
ymax=@（x）1-x;
volume=quad2d（h,0,1,0,ymax）
%精确解
exact_solution=pi/4-1/2
%一维差分diff
x=0:
0.1:
1;
y=[-0.4471.9783.286.167.087.347.669.569.489.311.2];
diff（x）
diff（y）
y（2:
end）-y（1:
end-1）
dydx=diff（y）./diff（x）
plot（x,y,'.-',x（1:
end-1）,dydx,'-or'）;
legend（'y（x）','dydx'）;
%微分描述函数局部特性，对函数的波形变化非常敏感。
%处理离散数据时，最好先对数据进行拟合，再进行微分。
p=polyfit（x,y,2）
dp=polyder（p）
dydx2=polyval（dp,x）
holdon
plot（x,dydx2,'--'）
legend（'y（x）','dydx','dydx2'）
holdoff
%一阶微分的几种格式
向前差分：
（yi+1–yi）/（xi+1–xi）
向后差分：
（yi–yi-1）/（xi–xi-1）
中心差分：
（yi+1–yi-1）/（xi+1–xi-1）
%比较三种差分格式
step=1;%trystep=0.5,0.25,0.1…
x=0:
step:
5;y=sin（x）;
x1=x（1:
end-1）;dydx1=diff（y）/step;
x2=x（2:
end）;dydx2=diff（y）/step;
x3=x（2:
end-1）;dydx3=（y（3:
end）-y（1:
end-2））/2/step;
x0=0:
0.01:
5;y0=cos（x0）;
plot（x0,y0,x1,dydx1,'.-r',x2,dydx2,'.-g',x3,dydx3,'-ok'）;
legend（'精确值','向前差分','向后差分','中心差分'）
%梯度函数gradient
%在绘图中的应用：
流线图
[x,y,z]=peaks（20）;
dx=x（1,2）-x（1,1）;dy=y（2,1）-y（1,1）;
[dzdx,dzdy]=gradient（z,dx,dy）;
c=contour（x,y,z,10）;clabel（c）;
holdon
quiver（x,y,dzdx,dzdy）;
holdoff
%另一个例子
[x,y]=meshgrid（-2:
0.2:
2）;
z=x.*exp（-x.^2-y.^2）;
[u,v]=gradient（z,0.2,0.2）;
c=contour（x,y,z,10）;clabel（c）;
holdon,quiver（x,y,u,v）;holdoff
%离散Laplace算子del2
%调用格式：
%L=del2（f,hx,hy,…）
%L=2f/2n，n~维数
%曲面的表面曲率
[x,y,z]=peaks;
dx=x（1,2）-x（1,1）;dy=y（2,1）-y（1,1）;
L=del2（z,dx,dy）;
surf（x,y,z,abs（L））
shadinginterp
数据统计
常用函数
✓sum：
求和
✓mean：
平均值（数学期望）
✓median：
中数
✓mode：
出现最多的数值
✓max/min：
最大值/最小值
✓var：
方差（无偏估计量）
✓std：
标准差（无偏估计量）
✓cov：
协方差（矩阵）
✓corrcoef：
协方差系数（矩阵）
%对多维数组，按列优先原则处理。
%方差函数中，输入参数1为计算方差，无偏估计缺省为0。
%协助方差函数只能按列处理，列数就是随机变量数目。
%%例：
单变量的方差和标准差
a=1:
10;
var（a）%缺省：
无偏估计
var（a,0）%无偏估计/n-1
var（a,1）%方差/n
v0=sum（（a-mean（a））.^2）/（length（a）-1）
v1=sum（（a-mean（a））.^2）/length（a）
std（a）,std（a,0）,std（a,1）
std（a）^2,std（a,1）^2%标准差的平方=方差
%%多变量的计算
a=magic（3）
var（a）
cov（a）%协方差矩阵
corrcoef（a）
a1=a（:
1）;a2=a（:
2）;
cov12=sum（（a1-mean（a1））.*（a2-mean（a2）））/（length（a1）-1）
%协方差如需按行计算，先做转置运算
cov（a'）
%例：
三个城市某月的日最高气温数据
t=[12818
15922
12519
14823
12622
11919
15915
81020
19718
12718
141019
11817
9723
8819
15818
8920
10717
12722
9819
12821
12820
10917
131218
91020
10622
14721
12522
13718
151023
131124
121222];
scales=size（t）
plot（t）,legend（'City-A','City-B','City-C'）;xlabel（'Days'）;ylabel（'Temperature'）;
avg=mean（t）%平均气温
dev=std（t）%气温波动
%每日气温与月平均值的差异
ft=t-avg（ones（length（t）,1）,:
）%...
figure,plot（ft）,legend（'City-A','City-B','City-C'）;
%出现频率最高的气温值
[n,f,c]=mode（t）%n：
频次最高的（最小）数值，f：
最高频次
%c：
频次最高的（全部）数值
%各城市温度的相关性
corrcoef（t）
%连续两天的温度变化量
dt=diff（t）;
plot（dt）,legend（'City-A','City-B','City-C'）;
mean（dt）,std（dt）,corrcoef（dt）
%残缺数表的处理
%读取数据，如data=xlsread（'*.xlsx'）
%生成残缺数表
t（5,1）=NaN;t（10,2）=NaN;t（15,2）=NaN;
plot（t）;
mean（t）,std（t）%nan运算结果不确定
%%编程处理
%对sum,mean等函数可以按列分别处理
t2=t（:
2）;%提取一列
t2（find（isnan（t2）））=[];%剔除nan
size（t2）
sum（t2）,var（t2）
%%Matlab提供的nan-簇函数，可以直接处理含缺数据
nansum（t）,nanvar（t）
%nanvar处理剔除nan后的一列数据。
%nancov剔除含nan的整行数据，保持矩阵结构。
%验证nanvar和nancov的处理方法：
a=magic（4）;b=a;
a（1,1）=NaN;c=a（2:
4,:
）;
nanvar（a）,var（b）,var（c）
nancov（a）,cov（b）,cov（c）
数据插值
三次样条插值方法：
给定数据系列（xi,yi）,i=1,2,…,N，用（N-1）个三阶多项式，逐段逼近。
插值函数f的约束条件：
yi=f（xi）
f（xi+）=f（xi-）
f（xi+）'=f（xi-）'
f（xi+）''=f（xi-）''
边界区域f'''相等
三次样条插值函数---spline
格式：
yy=spline（x,y,xx）
pp=spline（x,y）
%（x,y）是离散的观测数据（x即间断点），xx待插值点
%yy为插值，pp为插值的三次多项式
%y可以是矩阵或多维数组
x=linspace（0,2*pi,6）,y=sin（x）
xx=linspace（0,2*pi）;%
升级会员 