ch4 数值计算Word下载.docx
- 文档编号:15898259
- 上传时间:2022-11-16
- 格式:DOCX
- 页数:14
- 大小:73.82KB
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;
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- ch4 数值计算 数值 计算