第 4 章 数值计算
与符号计算相比，数值计算在科研和工程中的应用更为广泛。MATLAB也正是凭借其卓越的数值计算能力而称雄世界。随着科研领域、工程实践的数字化进程的深入，具有数字化本质的数值计算就显得愈益重要。 较之十年、二十年前，在计算机软硬件的支持下当今人们所能拥有的计算能力已经得到了巨大的提升。这就自然地激发了人们从新的计算能力出发去学习、理解概念的欲望，鼓舞了人们用新计算能力试探解决实际问题的雄心。 鉴于当今高校本科教学偏重符号计算和便于手算简单示例的实际，也出于帮助读者克服对数值计算生疏感的考虑，本章在内容安排上仍从“微积分”开始。这一方面与第2章符号计算相呼应，另方面通过“微积分”说明数值计算离散本质的微观和宏观影响。 为便于读者学习，本章内容的展开脉络基本上沿循高校数学教程，而内容深度力求控制在高校本科水平。考虑到知识的跳跃和交叉，本书对重要概念、算式、指令进行了尽可能完整地说明。
4.1
4.1.1
数值微积分
近似数值极限及导数
1?cos2xsinx，f2(x)?，试用机器零阈值eps替代理论0计算
xsinxx极限L1(0)?limf1(x)，L2(0)?limf2(x)。
【例4.1-1】设f1(x)?x?0x?0
x=eps;
L1=(1-cos(2*x))/(x*sin(x)), L2=sin(x)/x, L1 = 0 L2 =
1
syms t
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
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) hold on
【例4.1-2】已知x?sin(t)，求该函数在区间 [0, 2?]中的近似导函数。
1
plot(t,dxdt_eps) hold off
legend('x(t)','dx/dt') xlabel('t')
图 4.1-1 增量过小引起有效数字严重丢失后的毛刺曲线
x_d=sin(t+d);
dxdt_d=(x_d-x)/d; plot(t,x,'LineWidth',5) hold on
plot(t,dxdt_d) hold off
legend('x(t)','dx/dt') xlabel('t')
图 4.1-2 增量适当所得导函数比较光滑
【例4.1-3】已知x?sin(t)，采用diff和gradient计算该函数在区间 [0, 2?]中的近似导函数。。
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') hold on
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])
2
title('[0, 2\\pi]')
legend('x(t)','dxdt_{grad}','dxdt_{diff}','Location','North') xlabel('t'),box off hold off
subplot(1,2,2)
kk=(length(t)-10):length(t); hold on
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'),box off hold off
图 4.1-3 diff和gradient求数值近似导数的异同比较
4.1.2 数值求和与近似数值积分
【例 4.1-4】求积分s(x)???/20y(t)dt，其中y?0.2?sin(t)。
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') hold on
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]) hold off shg
sum求得积分 trapz求得积分 1.5762 1.3013
3
图 4.1-4 sum 和trapz求积模式示意
4.1.3 计算精度可控的数值积分
【例 4.1-5】求I??e01?x2dx 。
syms x
Isym=vpa(int(exp(-x^2),x,0,1)) Isym =
.74682413281242702539946743613185
format long d=0.001;x=0:d:1;
Itrapz=d*trapz(exp(-x.*x)) Itrapz =
0.74682407149919
fx='exp(-x.^2)'; Ic=quad(fx,0,1,1e-8) Ic =
0.74682413285445
【例 4.1-6】求s??? 1 2 1 0xydxdy。
syms x y
s=vpa(int(int(x^y,x,0,1),y,1,2)) s =
.40546510810816438197801311546432
format long
s_n=dblquad(@(x,y)x.^y,0,1,1,2) s_n =
0.40546626724351
4.1.4
函数极值的数值求解
【例4.1-7】已知
y?(x??)?esin(x??)，在??/2?x??/2区间，求函数的极小值。
4
（1）
syms x
y=(x+pi)*exp(abs(sin(x+pi))); yd=diff(y,x); xs0=solve(yd) yd_xs0=vpa(subs(yd,x,xs0),6) y_xs0=vpa(subs(y,x,xs0),6) y_m_pi=vpa(subs(y,x,-pi/2),6) y_p_pi=vpa(subs(y,x,pi/2),6) xs0 =
-1.0676598444985783372948670854801 yd_xs0 = .1e-4 y_xs0 = 4.98043 y_m_pi = 4.26987 y_p_pi = 12.8096
（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: 24
algorithm: 'golden section search, parabolic interpolation' message: [1x112 char]
（3）
xx=-pi/2:pi/200:pi/2; yxx=(xx+pi).*exp(abs(sin(xx+pi))); plot(xx,yxx)
xlabel('x'),grid on
图 4.1-5 在[-pi/2,pi/2]区间中的函数曲线
5