第四章 数值计算的应用

函数文件dydt.m （2）根据上述方程组，编写 函数文件 ）根据上述方程组，编写M函数文件 function ydot=dydt(t,y) ydot=[y(2);2*(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)) xlabel('t'),title('x(t)') % ode45(@dydt,tspan,y0)
s1=dblquad(fx1,0,1,1,2) s2=dblquad(fx2,0,1,1,2) s3=dblquad(fx3,0,1,1,2) s4=dblquad(@fx4,0,1,1,2) sn = 0.40546626724351
4.1.4 函数极值的数值求解
[x,fval,exitflag,output]=fminbnd(fun,x1,x2,options)
S2=dblquad(fun,xmin,xmax,ymin,ymax,tol) S3=triplequad(fun,xmin,xmax,ymin,ymax,zmin ,zmax,tol)
说明： 说明： 是被积函数， 字符串、 （1）fun是被积函数，可以是字符串、内联对象、匿 ） 是被积函数 可以是字符串 内联对象、 名函数和 函数文件的函数句柄 函数文件的函数句柄。 名函数和M函数文件的函数句柄。 （2）要保证对于向量形式的自变量（一般为字母 ） ）要保证对于向量形式的自变量（一般为字母x） 输入，输出为长度相同的函数值向量。 输入，输出为长度相同的函数值向量。 （3）a，b为一重积分的下限和上限； 为一重积分的下限和上限； ） ， 为一重积分的下限和上限 xmin,xmax,ymin,ymax,zmin,zmax是多重积分由 是多重积分由 内向外的积分限。 内向外的积分限。 是个标量， （4）tol是个标量，用来控制绝对误差。默认时，积 ） 是个标量 用来控制绝对误差。默认时， 分的绝对精度为10-6.（1e-6） 分的绝对精度为 （ ） 注：在编写被积函数时要遵循“数组运算”规则。 在编写被积函数时要遵循“数组运算”规则。
说明： 说明： （1）本例中的指令 ）本例中的指令s_ta=d*trapz(y); 可以用
s_ta=trapz(t,y)替换，即计算由t，y所绘折线下的面积。 替换，即计算由 ， 所绘折线下的面积 所绘折线下的面积。 替换 用作近似积分是错误的。 （2）s=sum(y); s_sa=d*s; 用作近似积分是错误的。阶梯 ） 虚线所占的自变量区间比积分区间多一个采样子区间。 虚线所占的自变量区间比积分区间多一个采样子区间。即 不能把d *sum(y)看作“矩形近似积分”。 看作“ 不能把 看作 矩形近似积分” （3）实际使用中，应该把子区间划分得相当小以便获取较 ）实际使用中， 高精度的近似积分。显然，采样点愈多，积分精度愈高， 高精度的近似积分。显然，采样点愈多，积分精度愈高， 但精度无法定量确定。 但精度无法定量确定。
【例4.1-9】求微分方程 】
d 2x 2 dx − 2(1 − x ) + x = 0 在初始条件 2 dt dt
dx(0) 情况下的解，并图示。 x (0) = 1, = 0 情况下的解，并图示。 dt
（1）把高阶微分方程改写成一阶微分方程 ）
dy1 dt y2 y1 (0) 1 dx , 令 y1 = x, y2 = ，则 dy = = 0 2 2 2(1 − y1 ) y2 − y1 y2 (0) dt dt
分方程组形式。 分方程组形式。
y ' = f ( y, t)
【例】解微分方程 y ' = y − 2t / y, y (0) = 1, 0 < t < 4 odefun=inline('y-2*t/y','t','y'); [t,y]=ode45(odefun,[0,4],1); plot(t,y) % ode45(odefun,[0,4],1)
4.1.3 计算精度可控的数值积分
数值积分有闭型（ 数值积分有闭型（closed－type）算法、开型 － ）算法、 （open－type）算法，其主要区别在于：是否需 － ）算法，其主要区别在于： 要计算积分区间端点处的函数值。 要计算积分区间端点处的函数值。 S1=quad(fun,a,b,tol) S1=quadl(fun,a,b,tol) Simpson法 法 Lobatto法 Lobatto法
format long fx1=‘x.^2’ fx2=inline('x.^y') fx3=@(x,y)x.^y %采用字符串表示被积函数 采用字符串表示被积函数 %采用内联对象表示被积函数 采用内联对象表示被积函数 %采用匿名函数表示被积函数 采用匿名函数表示被积函数 function z=fx4(x,y) z=x.^y
【例】求二元函数在原点附近的最大值： 求二元函数在原点附近的最大值：
f ( x, y ) = 5 − x 4 − y 4 + 4 xy
注：fminsearch被优化函数 中的多元自变量应采用单一变量名 被优化函数fun中的多元自变量应采用单一变量名 被优化函数 的向量形式表达。 的向量形式表达。 ff=@(x)(x(1)^4+x(2)^4-4*x(1)*x(2)-5); x0=[0 0]; [x,g]=fminsearch(ff,x0)
【例4.1-7】已知， 】已知， 间
y = (x + π ) ⋅ e
sin(x +π )在区−π / 2 ≤ x ≤ π / 2
，求函数的极小值。 求函数的极小值。
三种方法： 三种方法： ）（2）优化算法（ （1）教科书方法（自学）（ ）优化算法（fminbnd ）教科书方法（自学）（ 指令）（3）图形法求极小值（ginput指令 指令） 指令）（3）图形法求极小值（ginput指令） ）（
注意： 注意：不要企图借助数值计算求取极限
4.1.2 数值求和与近似数值积分
sx=sum(X) st=trapz(x,y) 沿列方向求和 采样梯形法沿列方向求函数y 采样梯形法沿列方向求函数 关于自变量x的积分 关于自变量 的积分 sct=cumtrapz(x,y) 采样梯形法沿列方向求函数 y 关于自变量x的累计积分 关于自变量 的累计积分
求一元函数在区间（ 求一元函数在区间（x1,x2）中的极小值。 ）中的极小值。
[x,fval,exitflag,output]=fminsearch(fun,x0,options)
单纯形法求多元函数的极值点
fun—字符串、内联对象、匿名函数、函数文件句柄； 字符串、内联对象、匿名函数、函数文件句柄； 字符串 输出量x,fval分别是极值点和相应的目标函数极值； 分别是极值点和相应的目标函数极值； 输出量 分别是极值点和相应的目标函数极值 输出量exitflag若给出大于 的数，说明成功搜索； 若给出大于0的数 说明成功搜索； 输出量 若给出大于 的数， 输出量output给出具体的优化算法和迭代次数。 给出具体的优化算法和迭代次数。 输出量 给出具体的优化算法和迭代次数
[ 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]) title('[0, 2\pi]')
syms x Isym=vpa(int(exp(-x^2),x,0,1)) Isym = .74682413281242702539946743613185 format long d=0.001;x=0:d:1; Itrapz=trapz(x,exp(-x.*x)) Itrapz = 0.74682407149919 fx='exp(-x.^2)'; Ic=quad(fx,0,1,1e-8) Ic = 0.74682413285445
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 数值微积分 4.2 矩阵和代数方程 4.4 多项式运算
4.1数值微积分
4.1.1 近似数值极限及导数
数值计算中， 在MATLAB数值计算中，既没有专门的求极限指令， 数值计算中 既没有专门的求极限指令， 也没有专门的求导指令。 提供了与“ 也没有专门的求导指令。但MATLAB提供了与“求 提供了与 导”概念有关的“求差分”指令。 概念有关的“求差分”指令。