ch4数值计算.docx
- 文档编号:6166666
- 上传时间:2023-01-04
- 格式:DOCX
- 页数:46
- 大小:352.55KB
ch4数值计算.docx
《ch4数值计算.docx》由会员分享,可在线阅读,更多相关《ch4数值计算.docx(46页珍藏版)》请在冰豆网上搜索。
ch4数值计算
第4章数值计算
与符号计算相比,数值计算在科研和工程中的应用更为广泛。
MATLAB也正是凭借其卓越的数值计算能力而称雄世界。
随着科研领域、工程实践的数字化进程的深入,具有数字化本质的数值计算就显得愈益重要。
较之十年、二十年前,在计算机软硬件的支持下当今人们所能拥有的计算能力已经得到了巨大的提升。
这就自然地激发了人们从新的计算能力出发去学习、理解概念的欲望,鼓舞了人们用新计算能力试探解决实际问题的雄心。
鉴于当今高校本科教学偏重符号计算和便于手算简单示例的实际,也出于帮助读者克服对数值计算生疏感的考虑,本章在内容安排上仍从“微积分”开始。
这一方面与第2章符号计算相呼应,另方面通过“微积分”说明数值计算离散本质的微观和宏观影响。
为便于读者学习,本章内容的展开脉络基本上沿循高校数学教程,而内容深度力求控制在高校本科水平。
考虑到知识的跳跃和交叉,本书对重要概念、算式、指令进行了尽可能完整地说明。
.1数值微积分
.1.1近似数值极限及导数
〖说明〗
【例4.1-1】设
,
,试用最小正数realmin替代理论0计算极限
,
。
%
x=realmin;
L1=(1-cos(2*x))/(x*sin(x)),%
L2=sin(x)/x,%
L1=
NaN
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
〖说明〗
【例4.1-2】
1)
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】
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-3
〖说明〗
.1.2数值求和与近似数值积分
〖说明〗
【例4.1-4】
clear
d=pi/8;%
t=0:
d:
pi/2;%
y=0.2+sin(t);%
s=sum(y);%
s_sa=d*s;%<6>
s_ta=d*trapz(y);%<7>
disp(['sum求得积分',blanks(3),'trapz求得积分'])
disp([s_sa,s_ta])
t2=[t,t(end)+d];%
y2=[y,nan];%
hs=stairs(t2,y2,':
k','LineWidth',3);%<12>
holdon
ht=plot(t,y,'r','LineWidth',3);%<14>
stem(t,y)%
legend([hs,ht],'sum','trapz','location','best')<16>
axis([0,pi/2+d,0,1.5])%
holdoff
shg
sum求得积分trapz求得积分
1.57621.3013
图4.1-4
〖说明〗
.1.3计算精度可控的数值积分
〖说明〗
表4.2-1积分指令integral的选项名称、取值、及默认值
选项名
Name
允许值
Value
默认值
DefaultValue
含义
'AbsTol'
正实数
'RelTol'
正实数
'ArrayValued'
false
true
'Waypoints'
实数单调序列数组
复数序列数组
【例4.1-5】在
之间,求曲线
与
所夹区域的面积(参见图4.1-5)。
1)
x=linspace(0.01,1.2,50);%
g1=x.^(-0.2);g2=x.^5;%
plot(x,g1,'-r.',x,g2,'-b*')
axis([0,1.2,0,3])
legend('g_1(x)=1/x^{0.2}','g_2(x)=x^5','Location','North')
title('x位于[0,1]间的g_1(x)曲线与g_2(x)曲线所夹的区域')
图4.1-5
2)
formatlong
G1=@(x)x.^-0.2;%<8>
Q1=integral(G1,0,1)%<9>
G2=@(x)x.^5;%<10>
Q2=integral(G2,0,1)%<11>
S12=Q1-Q2%<12>
Q1=
1.250000027856048
Q2=
0.166********6667
S12=
.0833********
3)
G=@(x)x.^[-0.2;5];%<13>
Q=integral(G,0,1,'ArrayValued',true)%<14>
S=[1,-1]*Q%<15>
Q=
1.250000027856048
0.166********6667
S=
1.083333361189381
4)
symst%
Gsym=vpa(int(t.^[-0.2;5],0,1));%
Ssym=Gsym
(1)-Gsym
(2)%
Ssym=
1.0833333333333333333333333333333
〖说明〗
【例4.1-6】验证
。
图4.1-6
F=@(z)(2*z-1)./(z.*(z-1));%<1>
Path=[2+1i,-1+1i,-1-1i,2-1i,2+1i];%<2>
%
sf=integral(F,2+1i,2+1i,'Waypoints',Path)
sf=
0.000000000000000+12.566370614359174i
〖说明〗
.1.4函数极值的数值求解
〖说明〗
【例4.1-7】
1)
x1=-50;x2=5;%
yx=@(x)(sin(x).^2.*exp(-0.1*x)-0.5*sin(x).*(x+0.1));%<2>
%
[xc0,fc0,exitflag]=fminbnd(yx,x1,x2)%<3>
%
xc0=
-8.4867
fc0=
-1.8621
exitflag=
1
2)
xx=-50:
pi/200:
5;%
ezplot(yx,[-50,5])%<5>
xlabel('x'),gridon
图4.1-5
3)
xx=[-23,-20,-18];%
fc=fc0;xc=xc0;%
fork=1:
2
[xw,fw]=fminbnd(yx,xx(k),xx(k+1));%<16>
iffw xc=xw; fc=fw; end end fprintf('函数最小值%6.5f发生在x=%6.5f处',fc,xc) 函数最小值-3.34765发生在x=-19.60721处 〖说明〗 【例4.1-8】求 的极小值点。 1) ff=@(x)(100*(x (2)-x (1)^2)^2+(1-x (1))^2); 2) formatshortg x0=[-5,-2,2,5;-5,-2,2,5];% [sx,sfval,sexit,soutput]=fminsearch(ff,x0) % sx= 0.99998-0.689710.415078.0886 0.99997-1.91684.96437.8004 sfval= 2.4112e-10 sexit= 1 soutput= iterations: 384 funcCount: 615 algorithm: 'Nelder-Meadsimplexdirectsearch' message: [1x194char] 3) formatshorte disp([ff(sx(: 1)),ff(sx(: 2)),ff(sx(: 3)),ff(sx(: 4))]) 2.4112e-105.7525e+022.2967e+033.3211e+05 〖说明〗 .1.5常微分方程的数值解 〖说明〗 【例4.1-9】求微分方程 , ,在初始条件 , 情况下的解,并图示。 1) 2) 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);%<3> plot(tt,yy(: 1)) xlabel('t'),title('x(t)') 图4.1-6 4) plot(yy(: 1),yy(: 2))% xlabel('位移'),ylabel('速度') 图4.1-7 〖说明〗 .2矩阵和代数方程 .2.1矩阵的标量特征参数 表4.2-2 术语 数学含义 MATLAB指令 秩 Rank rank(A) 迹 Trace trace(A) 行列式Determinant det(A) 【例4.2-1】 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-2】 formatshort% rngdefault% 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.7479 tBA= 3.7479 tCD= 3.3399 tDC= 3.3399 % d_A_B=det(A)*det(B) dAB=det(A*B) dBA=det(B*A) d_A_B= -0.0852 dAB= -0.0852 dBA= -0.0852 % % % dCD=det(C*D) dDC=det(D*C)% dCD= -0.0557 dDC= 1.5085e-017 .2.2矩阵的变换和特征值分解 〖说明〗 【例4.2-3】 1) A=magic(4)% [R,ci]=rref(A)% A= 162313 511108 97612 414151 R= 1001 0103 001-3 0000 ci= 123 2) r_A=length(ci) r_A= 3 3) aa=A(: 1: 3)*R(1: 3,4)% err=norm(A(: 4)-aa)% aa= 13 8 12 1 err= 0 【例4.2-4】 A=reshape(1: 15,5,3);% X=null(A)% S=A*X% n=size(A,2);% l=size(X,2);% n-l==rank(A)% X= -0.4082 0.8165 -0.4082 S= 1.0e-014* 0.2665 0.2665 0.3553 0.4441 0.6217 ans= 1 【例4.2-5】 1) 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) [VR,DR]=cdf2rdf(V,D) VR= 0.77460 0.0430-0.6310 DR= 0.83332.4438 -2.44380.8333 3) A1=V*D/V% A1_1=real(A1)% A2=VR*DR/VR err1=norm(A-A1,'fro') err2=norm(A-A2,'fro') A1= 1.0000-3.0000-0.0000i 2.00000.6667 A1_1= 1.0000-3.0000 2.00000.6667 A2= 1.0000-3.0000 2.00000.6667 err1= 4.6613e-016 err2= 4.4409e-016 .2.3线性方程的解 101线性方程解的一般结论 102除法运算解方程 〖说明〗 【例4.2-6】求方程 的解。 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.1374e-014 103矩阵逆 〖说明〗 【例4.2-7】 1) rngdefault A=gallery('randsvd',300,2e13,2);% x=ones(300,1);% b=A*x;% cond(A)% ans= 2.0022e+013 2) tic% xi=inv(A)*b;% ti=toc% eri=norm(x-xi)% rei=norm(A*xi-b)/norm(b)% ti= 0.2034 eri= 0.0812 rei= 0.0047 3) tic; xd=A\b;% td=toc erd=norm(x-xd) red=norm(A*xd-b)/norm(b) td= 0.0125 erd= 0.0274 red= 7.5585e-015 〖说明〗 .2.4一般代数方程的解 解题步骤大致如下: (1) (2) 〖说明〗 【例4.2-8】 1) symst ft=sin(t)^2*exp(-0.1*t)-0.5*abs(t); S=solve(ft,t)%<3> ftS=subs(ft,t,S)% S= 0 ftS= 0 2) % y_C=inline('sin(t).^2.*exp(-0.1*t)-0.5*abs(t)','t'); % % t=-10: 0.01: 10;% Y=y_C(t);% clf, plot(t,Y,'r'); holdon plot(t,zeros(size(t)),'k');% xlabel('t');ylabel('y(t)') holdoff 图4.2-1 % zoomon% [tt,yy]=ginput(5);% zoomoff% 图4.2-2 % tt% tt= -2.0039 -0.5184 -0.0042 0.6052 1.6717 % [t4,y4]=fzero(y_C,0.1)%<18> t4= 0.5993 y4= 1.1102e-016 〖说明〗 .3概率分布和统计分析 .3.1概率函数、分布函数、逆分布函数和随机数的发生 101二项分布(Binomialdistribution) 〖说明〗 【例4.3-1】画出N=100,p=0.5情况下的二项分布概率特性曲线。 N=100;p=0.5;% k=0: N;% pdf=binopdf(k,N,p);% cdf=binocdf(k,N,p);% h=plotyy(k,pdf,k,cdf);%<5> set(get(h (1),'Children'),'Color','b','Marker','.','MarkerSize',13) % set(get(h (1),'Ylabel'),'String','pdf') % set(h (2),'Ycolor',[1,0,0]) % set(get(h (2),'Children'),'Color','r','Marker','+','MarkerSize',4) % set(get(h (2),'Ylabel'),'String','cdf') %<10> xlabel('k')% gridon 图4.3-1 〖说明〗 102正态分布(Normaldistribution) 〖说明〗 【例4.3-2】 mu=3;sigma=0.5;% x=mu+sigma*[-3: -1,1: 3]; yf=normcdf(x,mu,sigma); P=[yf(4)-yf(3),yf(5)-yf (2),yf(6)-yf (1)]; % xd=1: 0.1: 5; yd=normpdf(xd,mu,sigma);% clf fork=1: 3 %------------------------------- xx=x(4-k): sigma/10: x(3+k); yy=normpdf(xx,mu,sigma); %-------------------------------------------------- subplot(3,1,k),plot(xd,yd,'b');% holdon fill([x(4-k),xx,x(3+k)],[0,yy,0],'g');% holdoff ifk<2 text(3.8,0.6,'[{\mu}-{\sigma},{\mu}+{\sigma}]') else kk=int2str(k); text(3.8,0.6,['[{\mu}-',kk,'{\sigma},{\mu}+',kk,'{\sigma}]']) end text(2.8,0.3,num2str(P(k)));shg end xlabel('x');shg 图4.3-2 〖说明〗 103各种概率分布的交互式观察界面 【例4.3-3】概率分布交互界面使用方法介绍。 1) 2) 3) 4) 图4.3-3 〖说明〗 .3.2全局随机流、随机数组和统计分析 101全局随机流的操控 〖说明〗 表4.3-1 generator 算例 可取字符串 含义 'twister' 'v5uniform' 'v5normal' 'v4' 102三个基本随机数组创建指令 〖说明〗 【例4.3-4】 1) rngdefault
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- ch4 数值 计算