湖南大学matlab教案04.docx
《湖南大学matlab教案04.docx》由会员分享，可在线阅读，更多相关《湖南大学matlab教案04.docx（48页珍藏版）》请在冰点文库上搜索。
湖南大学matlab教案04
第4章数值计算
与符号计算相比，数值计算在科研和工程中的应用更为广泛。
MATLAB也正是凭借其卓越的数值计算能力而称雄世界。
随着科研领域、工程实践的数字化进程的深入，具有数字化本质的数值计算就显得愈益重要。
.1数值微积分
.1.1近似数值极限及导数
Matlab数值计算中，没有提供专门的求极限和求导的指令，但提供了与求导数相关的求差分指令。
1、差分
对于向量x（n），一阶差分为:
[x
（2）-x
（1）,x（3）-x
（2）,…,x（n）-x（n-1）],其是相邻元素相减。
二阶差分是相邻行向量相减。
指令：
DX=diff（X）一阶差分
DX=diff（X,n）n阶差分
DX=diff（x,n,dim）在dim指定的维上，
求n阶差分
2、数值梯度
FX=gradient（F,h）
求一元函数的行元素梯度。
F为向量时，
FX
（1）=[F
（2）-F
（1）]/h
FX（end）=[F（end）-F（end-1）]/h,
FX（2:
end-1）=[（F（3:
end）-F（1:
end-2））/2]/h
h是沿坐标取点的步长，默认步长为1。
[FX,FY]=gradient（F，h）
求二元函数的梯度，FX,FY与F的大小相同，分
别表示行、列元素间的“梯度”，h是沿坐标取
点的步长，默认步长为1。
【例4.1-1】设
，
，试用机器零阈值eps替代理论0计算极限
，
。
x=eps;
L1=（1-cos（2*x））/（x*sin（x））,
L2=sin（x）/x,
L1=
0
L2=
1
symst
f1=（1-cos（2*t））/（t*sin（t））;
f2=sin（t）/t;
Ls1=limit（f1,t,0）
Ls2=limit（f2,t,0）
Ls1=
2
Ls2=
1
注意：
尽量不用数值计算求极限
【例】用一个简单矩阵表现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）
F=
123
456
789
Dx=
333
333
Dx_2=
11
11
11
FX=
111
111
111
FY=
333
333
333
FX_2=
222
222
222
FY_2=
666
666
666
【例4.1-2】已知
，求该函数在区间
中的近似导函数。
d=pi/100;
t=0:
d:
2*pi;
x=sin（t）;
dt=5*eps;
x_eps=sin（t+dt）;
dxdt_eps=（x_eps-x）/dt;
plot（t,x,'LineWidth',5）
holdon
plot（t,dxdt_eps）
holdoff
legend（'x（t）','dx/dt'）
xlabel（'t'）
图4.1-1增量过小引起有效数字严重丢失后的毛刺曲线
（2）增量取适当
x_d=sin（t+d）;
dxdt_d=（x_d-x）/d;
plot（t,x,'LineWidth',5）
holdon
plot（t,dxdt_d）
holdoff
legend（'x（t）','dx/dt'）
xlabel（'t'）
图4.1-2增量适当所得导函数比较光滑
注意：
数值导数的使用应十分谨慎，自变量增量的取值应适当
【例4.1-3】已知
，采用diff和gradient计算该函数在区间
中的近似导函数。
。
clf
d=pi/100;
t=0:
d:
2*pi;
x=sin（t）;
dxdt_diff=diff（x）/d;
dxdt_grad=gradient（x）/d;
subplot（1,2,1）
plot（t,x,'b'）
holdon
plot（t,dxdt_grad,'m','LineWidth',8）
plot（t（1:
end-1）,dxdt_diff,'.k','MarkerSize',8）
axis（[0,2*pi,-1.1,1.1]）
title（'[0,2\pi]'）
legend（'x（t）','dxdt_{grad}','dxdt_{diff}','Location','North'）
xlabel（'t'）,boxoff
holdoff
subplot（1,2,2）
kk=（length（t）-10）:
length（t）;
holdon
plot（t（kk）,dxdt_grad（kk）,'om','MarkerSize',8）
plot（t（kk-1）,dxdt_diff（kk-1）,'.k','MarkerSize',8）
title（'[end-10,end]'）
legend（'dxdt_{grad}','dxdt_{diff}','Location','SouthEast'）
xlabel（'t'）,boxoff
holdoff
图4.1-3diff和gradient求数值近似导数的异同比较
.1.2数值求和与近似数值积分
指令：
1：
SX=sum（X）
按列求和，生成一各行向量，各元素是X各列元素之和。
2：
Scs=cumsum（X）
Scs大小与X的相同，Sxs（I,k）是X的第k列的前i个元素之和。
3：
St=trapz（x,y）
给出采样点（x,y）所连折线下的面积，即函数y在自变量区间x上的近似积分
4：
Sct=cumtrapz（x,y）
求累计积分
[例]
a=[123;456;786],b=sum（a）,c=cumsum（a）
a=
123
456
786
b=
121515
c=
123
579
121515
【例4.1-4】求积分
，其中
。
clear
d=pi/8;
t=0:
d:
pi/2;
y=0.2+sin（t）;
s=sum（y）;
s_sa=d*s;%小矩形面积之和
s_ta=d*trapz（y）;%折线下面积
disp（['sum求得积分',blanks（3）,'trapz求得积分']）
disp（[s_sa,s_ta]）
t2=[t,t（end）+d];
y2=[y,nan];
stairs（t2,y2,':
k'）
holdon
plot（t,y,'r','LineWidth',3）
h=stem（t,y,'LineWidth',2）;
set（h
（1）,'MarkerSize',10）
axis（[0,pi/2+d,0,1.5]）
holdoff
shg
sum求得积分trapz求得积分
1.57621.3013
图4.1-4sum和trapz求积模式示意
注意：
以上求积分办法精度难控制，不常用。
.1.3计算精度可控的数值积分
S1=quad（fun,a,b,tol,）Simpson法
s1=quad1（fun,a,b,tol,）lobatto法
功能：
a,b为上下限，tol为绝对误差，默认为10-6
S2=dblquad（f,xmin,xmax,ymin,ymax,tol,）
该函数求f（x,y）在[xmin,xmax]×[ymin,ymax]区域上的二重定积分。
参数tol的用法与函数quad完全相同
S3=triplequad（f,xmin,xmax,ymin,ymax,tol,）
该函数求f（x,y）在[xmin,xmax]×[ymin,ymax]区域上的三重定积分。
参数tol的用法与函数quad完全相同
【例4.1-5】求
。
方法1：
symsx
Isym=vpa（int（exp（-x^2）,x,0,1））
Isym=
.74682413281242702539946743613185
方法2：
formatlong
d=0.001;x=0:
d:
1;
Itrapz=d*trapz（exp（-x.*x））
Itrapz=
0.746824071499185
方法3：
fx='exp（-x.^2）';
Ic=quad（fx,0,1,1e-8）
Ic=
0.746824132854452
【例4.1-6】求
。
symsxy
s=vpa（int（int（x^y,x,0,1）,y,1,2））
s=
.40546510810816438197801311546432
formatlong
s_n=dblquad（@（x,y）x.^y,0,1,1,2）或s_n=dblquad（‘x.^y’,0,1,1,2）
s_n=
0.405466267243508
.1.4函数极值的数值求解
调用格式
1:
[x,fval,exitflag,output]=fminbnd（fun,x1,x2,option）
其中fminbnd函数用于求一元函数的最小值点。
fun是目标函数名，x1和x2限定自变量的取值范围。
四个输出参数分别表示：
极小值点；函数值；判断是否成功搜到极值点（是否大于0）；算法。
2:
[x,fval,exitflag,output]=fminsearch（fun,x0,option）
求多元函数fun从x0开始的极小值x,该点的函数值为fval.
【例4.1-7】已知
，在
区间，求函数的极小值。
（2）优化算法
x1=-pi/2;x2=pi/2;
yx=@（x）（x+pi）*exp（abs（sin（x+pi）））;
[xn0,fval,exitflag,output]=fminbnd（yx,x1,x2）
xn0=
-1.2999e-005
fval=
3.1416
exitflag=
1
output=
iterations:
21
funcCount:
22
algorithm:
[1x46char]
message:
[1x112char]
（3）图形法
xx=-pi/2:
pi/200:
pi/2;
yxx=（xx+pi）.*exp（abs（sin（xx+pi）））;
plot（xx,yxx）
xlabel（'x'）,gridon
图4.1-5在[-pi/2,pi/2]区间中的函数曲线
[xx,yy]=ginput
（1）
%执行后吧十字点对准极值点点击
xx=
1.5054e-008
yy=
3.1416
图4.1-6函数极值点附近的局部放大和交互式取值
【例4.1-8】求
的极小值点。
它即是著名的Rosenbrock's"Banana"测试函数，它的理论极小值是
。
（1）
ff=@（x）（100*（x
（2）-x
（1）.^2）^2+（1-x
（1））^2）;
（2）
x0=[-5,-2,2,5;-5,-2,2,5];%提供4个搜索起点
[sx,sfval,sexit,soutput]=fminsearch（ff,x0）
sx=
1.0000-0.68970.41518.0886
1.0000-1.91684.96437.8004
sfval=
2.4112e-010
sexit=
1
soutput=
iterations:
384
funcCount:
615
algorithm:
[1x33char]
message:
[1x196char]
（3）检查目标函数值
formatshorte
disp（[ff（sx（:
1））,ff（sx（:
2））,ff（sx（:
3））,ff（sx（:
4））]）
Columns1through3
2.4112e-0105.7525e+0022.2967e+003
Column4
3.3211e+005
.1.5常微分方程的数值解
一阶常微分方程的初值问题为求解下式：
对于高阶微分方程，需要运用变量变换转换为一阶微分方程，初始条件也作相应的变换。
设有高阶微分方程：
作变量替换，
上式为：
其中Y0、Y’和Y都为列向量。
Matlab指令格式：
[t,Y]=solver（‘F’,tspan,y0）
Solver的典型代表是ode45
以上括号中：
F是函数文件名；tspan是求数值解的时间区间，形如[t0,t1,…tf],必须单调排列；y0是初始条件列向量；t是所求数值解的自变量列向量,而Y（:
k）就是yk分量的解。
【例4.1-9】求微分方程
，
，在初始条件
情况下的解，并图示。
（1）高阶变一阶
令y1=x,y2=dx/dt,则原二阶方程可变为一阶方程组
（2）编写M函数文件
[DyDt.m]
functionydot=DyDt（t,y）
mu=2;
ydot=[y
（2）;mu*（1-y
（1）^2）*y
（2）-y
（1）];
（3）解算微分方程
tspan=[0,30];
y0=[1;0];
[tt,yy]=ode45（@DyDt,tspan,y0）;
plot（tt,yy（:
1））%y1的解
xlabel（'t'）,title（'x（t）'）
图4.1-7微分方程解
（4）画相平面图
plot（yy（:
1）,yy（:
2））
xlabel（'位移'）,ylabel（'速度'）
图4.1-8平面相轨迹
【例】求方程
初始值y（0）=-1,y’（0）=1,y’’（0）=0，时间区间[0，12]内的数值解
解答：
1、阶变换低阶：
令y1=y,y2=y’,y3=y’’,初始条件为y1（0）=-1,y2（0）=1,y3（0）=0
2、编写方程函数名为odet3.m
functiondy=odet3（t,y）
dy=zeros（3,1）;
dy
（1）=y
（2）;
dy
（2）=y（3）;
dy（3）=y
（1）.^（1/2）-2*y（3）.^2*（y
（2）-4）;
或dy=[y
（2）;y（3）;y
（1）.^（1/2）-2*y（3）.^2*（y
（2）-4）]
tspan=[012];
%在同一目录下，计算方程数值解。
给出时间区间和初始值
y0=[-110];
[t,Y]=ode45（'odet3',tspan,y0）;
%采用ode45算法
plot（t,Y（:
1）,'-',t,Y（:
2）,'-.',t,Y（:
3）,'.'）
%绘制计算结果并标注，如图3-15
legend（'Y1','Y2','Y3'）
.2矩阵和代数方程
.2.1矩阵运算和特征参数
10一矩阵运算
按矩阵运算规则计算
【例4.2-2】观察矩阵的转置操作和数组转置操作的差别。
clear;A=zeros（2,3）;
A（:
）=1:
6;
A=A*（1+i）
A_A=A.'
A_M=A'
A=
1.0000+1.0000i3.0000+3.0000i5.0000+5.0000i
2.0000+2.0000i4.0000+4.0000i6.0000+6.0000i
A_A=
1.0000+1.0000i2.0000+2.0000i
3.0000+3.0000i4.0000+4.0000i
5.0000+5.0000i6.0000+6.0000i
A_M=
1.0000-1.0000i2.0000-2.0000i
3.0000-3.0000i4.0000-4.0000i
5.0000-5.0000i6.0000-6.0000i
10二矩阵的标量特征参数
计算矩阵特征参数的指令：
Rank:
矩阵的秩
Trace:
矩阵的迹（主对角元素之和）
Det:
矩阵的行列式
【例4.2-3】矩阵标量特征参数计算示例。
A=reshape（1:
9,3,3）
r=rank（A）
d3=det（A）
d2=det（A（1:
2,1:
2））%子行列式
t=trace（A）
A=
147
258
369
r=
2
d3=
0
d2=
-3
t=
15
【例4.2-4】矩阵标量特征参数的性质。
formatshortg
rand（'state',0）
A=rand（3,3）;
B=rand（3,3）;
C=rand（3,4）;
D=rand（4,3）;
tAB=trace（A*B）
tBA=trace（B*A）
tCD=trace（C*D）
tDC=trace（D*C）
tAB=
3.6697
tBA=
3.6697
tCD=
2.1544
tDC=
2.1544
d_A_B=det（A）*det（B）
dAB=det（A*B）
dBA=det（B*A）
d_A_B=
0.0846
dAB=
0.0846
dBA=
0.0846
dCD=det（C*D）
dDC=det（D*C）%非同阶矩阵相乘交换位置行列式改变
dCD=
-0.012362
dDC=
-2.1145e-018
.2.2矩阵的变换和特征值分解
指令：
[R,ci]=rref（A）：
用初等变换把A编程行阶梯矩阵R,ci是行向量，其元素表示A的第‘几’列是基。
X=null（A）：
求A矩阵零空间的正交基，满足AX=0
eig（A）：
（1）d=eig（A）：
求矩阵A的全部特征值，（构成向量d）。
（2）[V,D]=eig（A）：
求矩阵A的全部特征值，构成对角阵D，并求A的特征向量构成V的列向量，使AV=VD成立.
【例4.2-5】行阶梯阵简化指令rref计算结果的含义。
（1）
A=magic（4）
[R,ci]=rref（A）
A=
162313
511108
97612
414151
R=
1001
0103
001-3
0000
ci=
123
%A的第1，2，3列是基
（2）
r_A=length（ci）%A的秩
r_A=
3
【例4.2-6】矩阵零空间及其含义。
A=reshape（1:
15,5,3）;
X=null（A）%A的零空间
S=A*X
n=size（A,2）;%求A的列数
l=size（X,2）;%求X的列数，即零空间的维数
n-l==rank（A）
X=
0.4082
-0.8165
0.4082
S=
1.0e-014*
-0.2665
-0.1776
-0.0888
-0.0888
-0.0888
ans=1
【例4.2-7】简单实阵的特征值分解。
（1）8
A=[1,-3;2,2/3]
[V,D]=eig（A）
A=
1.0000-3.0000
2.00000.6667
V=
0.77460.7746
0.0430-0.6310i0.0430+0.6310i
D=
0.8333+2.4438i0
00.8333-2.4438i
.2.3线性方程的解
10一线性方程解的一般结论
10二除法运算解方程
解方程
（1）方程AX=b，当A非奇异时，用左除法解方程
指令：
x=A\b
（2）方程XA=b，当A非奇异时，用右除法解方程
指令：
x=b\A
【例4.2-8】求方程
的解。
（1）
A=reshape（1:
12,4,3）;
b=（13:
16）';
（2）
ra=rank（A）
rab=rank（[A,b]）
ra=
2
rab=
2
（3）
xs=A\b;
xg=null（A）;
c=rand
（1）;
ba=A*（xs+c*xg）
norm（ba-b）
Warning:
Rankdeficient,rank=2,tol=1.8757e-014.
ba=
13.0000
14.0000
15.0000
16.0000
ans=
1.0658e-014
10三矩阵逆
指令：
inv（A）
对于方程AX=b的求解方法：
（1）左除法：
x=A\b
（2）逆矩阵法：
x=inv（A）*b
注意：
方法
（1）速度快，且精度比方法
（2）高得多
.2.4一般代数方程的解
对于函数f（x）=0求零点：
方法
（1）:
作图法
先确定一个零点可能存在的自变量空间，通过plot指令作图，观察f（x）与横轴的交点坐标，用zoom指令局部放大，借助ginput指令获得交点坐标值。
方法
（2）指令法
指令：
[x,favl]=fzero（fun,x0）
%求一元函数零点最简格式
[x,favl]x=fsolve（fun,x0）
%非线性方程法求解多元函数零点
注意：
x0为零点初始猜测值，可以通过作图先得到。
例：
求f（x）=x3-2x+5的零点
%第一步，绘制一元曲线图
x=-5:
0.1:
5;
f=x.^3-2*x-5;
xlabel（'x'）;
ylabel（'f（x）'）;
plot（x,f,x,0）%x轴划0刻度
%通过观察零点大致在2附近
编写并存储一元函数表达式，函数名为zero1.m，使f（x）=0
functionF=zero1（x）
F=x.^3-2*x-5
%计算2附近的零点和该点的函数值
[x,f]=fzero（@zero1,2）
x=
2.0946
f=
-8.8818e-016
【例4.2-10】求
的零点
升级会员 