ch4 数值计算Word下载.docx
《ch4 数值计算Word下载.docx》由会员分享，可在线阅读，更多相关《ch4 数值计算Word下载.docx（14页珍藏版）》请在冰豆网上搜索。
B=diag（diag（A））
C=A-B
5．先运行指令x=-3*pi:
pi/15:
3*pi;
y=x;
[X,Y]=meshgrid（x,y）;
warningoff;
Z=sin（X）.*sin（Y）./X./Y;
产生矩阵Z。
（1）请问矩阵Z中有多少个“非数”数据？
（2）用指令surf（X,Y,Z）;
shadinginterp观察所绘的图形。
（3）请写出绘制相应的“无裂缝”图形的全部指令。
x=-3*pi:
A=sum（sum（isnan（Z）））
surf（X,Y,Z）;
shadinginterp
xx=x+（x==0）*eps;
yy=y+（y==0）*eps;
[X,Y]=meshgrid（xx,yy）;
Z=sin（X）.*sin（Y）./X./Y;
shadinginterp;
6．下面有一段程序，企图用来解决如下计算任务：
有矩阵
，当
依次取10,9,8,7,6,5,4,3,2,1时，计算矩阵
“各列元素的和”，并把此求和结果存放为矩阵Sa的第k行。
例如
时，A阵为
，此时它各列元素的和是一个
行数组
，并把它保存为Sa的第3行。
问题：
该段程序的计算结果对吗？
假如计算结果不正确，请指出错误发生的根源，并改正之。
fork=10:
-1:
1
A=reshape（1:
10*k,k,10）;
Sa（k,:
）=sum（A）;
end
Sa
第4章数值计算
.1数值微积分
.1.1近似数值极限及导数
MATLAB数值计算中，没有专门的求极限和导数的指令。
原因：
数值精度有限，不能表示无穷小量，不能准确描述一个数的邻域。
【例4.1-1】设
，
，试用机器零阈值eps替代理论0计算极限
。
x=eps;
L1=（1-cos（2*x））/（x*sin（x））,
L2=sin（x）/x,
symst
f1=（1-cos（2*t））/（t*sin（t））;
f2=sin（t）/t;
Ls1=limit（f1,t,0）
Ls2=limit（f2,t,0）
理论分析表明：
，
●借助符号计算所求的极限与理论值一致；
●用数值法近似计算的极限与理论不一致。
提醒：
除非数值近似法求的极限经过理论验证，否则绝不要借助数值法求取极限。
函数
在点
处的导数
【例4.1-2】已知
，求该函数在区间
中的近似导函数。
d=pi/100;
t=0:
d:
2*pi;
x=sin（t）;
dt=5*eps;
x_eps=sin（t+dt）;
%自变量的增量为5eps。
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增量过小引起有效数字严重丢失后的毛刺曲线
x_d=sin（t+d）;
dxdt_d=（x_d-x）/d;
%以d=pi/100为增量算得的数值导数
plot（t,dxdt_d）
）
图4.1-2增量适当所得导函数比较光滑
结论：
自变量增量的选取一定要避免太小。
MATLAB提供的与“求导”概念有关的指令：
dx=diff（X）求差分
当X是向量时，dx=X（2:
n）–X（1:
n-1）；
当X是矩阵时，dx=X（2:
n,:
）–X（1:
n-1,:
）；
dx的长度比X的“长度”短一个元素（或一行）。
x=[12345];
y=diff（x）
y=
1111
diff（X,2）isthesameasdiff（diff（X））
z=diff（x,2）
z=
000
FX=gradient（F）求梯度
[FX,FY]=gradient（F）二元函数的梯度
当F是向量时，FX
（1）=F
（2）–F
（1）,
FX（end）=F（end）–F（end-1）,
FX（2:
end-1）=（F（3:
end）–F（1:
end-2））/2,
FX的“长度”与F相同。
当F是矩阵时，FX和FY是与F同样大小的矩阵。
FX的每行给出F相应行元素间的“梯度”；
FY的每列给出F相应列元素间的“梯度”。
【例4.1-3】已知
，采用diff和gradient计算该函数在区间
clf
dxdt_diff=diff（x）/d;
dxdt_grad=gradient（x）/d;
subplot（1,2,1）%宏观上看，diff和gradient所求得近似导数大体相同
b'
plot（t,dxdt_grad,'
m'
8）
plot（t（1:
end-1）,dxdt_diff,'
.k'
MarkerSize'
axis（[0,2*pi,-1.1,1.1]）
title（'
[0,2\pi]'
dxdt_{grad}'
dxdt_{diff}'
Location'
North'
）,boxoff
subplot（1,2,2）
kk=（length（t）-10）:
length（t）;
%最后11个数据的下标
plot（t（kk）,dxdt_grad（kk）,'
om'
8）%微观上看，不仅数值上有差异，
plot（t（kk-1）,dxdt_diff（kk-1）,'
8）%而且diff没有给出最后一点的导数。
[end-10,end]'
SouthEast'
holdoff
图4.1-3diff和gradient求数值近似导数的异同比较
.1.2数值求和与近似数值积分
Sx=sum（X）沿列方向求和
Scs=cumsum（X）沿列方向求累计和
St=trapz（x,y）采用梯形法沿列方向求函数y关于自变量x的积分
Sct=cumtrapz（x,y）采用梯形法沿列方向求函数y关于自变量x的累计积分
采样点数愈多，精度越高，但无法定量控制积分精度。
【例4.1-4】求积分
，其中
clear
d=pi/8;
pi/2;
y=0.2+sin（t）;
s=sum（y）;
s_sa=d*s;
s_ta=d*trapz（y）;
%可用s_ta=trapz（t,y）替换，此通用格式适用于“非等间隔采样”的情况。
disp（['
sum求得积分'
blanks（3）,'
trapz求得积分'
]）
disp（[s_sa,s_ta]）
t2=[t,t（end）+d];
y2=[y,nan];
stairs（t2,y2,'
:
k'
）
plot（t,y,'
r'
3）
h=stem（t,y,'
2）;
set（h,'
10）
axis（[0,pi/2+d,0,1.5]）%由图可见：
阶梯虚线所占的自变量区间比积分区间多一个采样子区间。
shg%showgraphwindow
图4.1-4sum和trapz求积模式示意
采样点数愈多，精度越高。
存在的问题：
无法预先设置欲求积分的精度。
.1.3计算精度可控的数值积分
数值积分：
闭型算法：
需要计算积分区间端点处的函数值
开型算法：
不需要计算积分区间端点处的函数值
闭型数值积分指令：
S1=quad（fun,a,b,tol）采用递推自适应Simpson法计算积分
S1=quadl（fun,a,b,tol）采用递推自适应Lobatto法计算积分
S2=dblquad（fun,xmin,xmax,ymin,ymax,tol）二重数值积分
S3=triplequad（fun,xmin,xmax,ymin,ymax,zmin,zmax,tol）三重数值积分
被积函数fun：
字符串、匿名函数、M函数文件的函数句柄等
被积函数的自变量一般采用字母x。
tol：
用来控制绝对误差，默认积分绝对精度为10-6。
【例4.1-5】求
symsx
Isym=vpa（int（exp（-x^2）,x,0,1））%给出32位精度的积分值
formatlong%15位数字显示
d=0.001;
x=0:
1;
Itrapz=d*trapz（exp（-x.*x））%梯形法
fx='
exp（-x.^2）'
;
%字符串，见附录A.1
Ic=quad（fx,0,1,1e-8）%控制绝对精度达1e-10
【例4.1-6】求
symsxy
s=vpa（int（int（x^y,x,0,1）,y,1,2））
%符号计算结果
formatlong
s_n=dblquad（'
x.^y'
0,1,1,2）%被积函数一定要写成“数组运算”格式
%数值积分结果
匿名函数
sqr=@（x）x.^2;
升级会员 