MATLAB数值计算与数据分析(4)

2.2.2 查表命令

命令1 table1

功能 一维查表

格式 Y = table1(TAB,X0) %返回用表格矩阵TAB中的行线性插值元素，对X0（TAB的第一列查找X0）进行线性插值得到的结果Y。矩阵TAB是第一列包含关键值，而其他列包含数据的矩阵。X0中的每一元素将相应地返回一线性插值行向量。矩阵TAB的第一列必须是单调的。

例2-38

>>tab = [(1:4)' hilb(4)]

>>y = table1(tab,[1 2.3 3.6 4])

查表结果为：

tab =

1.0000 1.0000 0.5000 0.3333 0.2500

2.0000 0.5000 0.3333 0.2500 0.2000

3.0000 0.3333 0.2500 0.2000 0.1667

4.0000 0.2500 0.2000 0.1667 0.1429

Warning: TABLE1 is obsolete and will be removed in future versions. Use INTERP1 or INTERP1Q instead.

> In D:\MATLABR12\toolbox\matlab\polyfun\table1.m at line 31

y =

1.0000 0.5000 0.3333 0.2500

0.4500 0.3083 0.2350 0.1900

0.2833 0.2200 0.1800 0.1524

0.2500 0.2000 0.1667 0.1429

由上面结果可知，table1是一将要废弃的命令。

命令2 table2

功能 二维查表

格式 Z = table1(TAB,X0,Y0) %返回用表格矩阵TAB中的行与列交叉线性线性插值元素，对X0（TAB的第一列查找X0）进行线性插值，对Y0（TAB的第一行查找Y0）进行线性插值，对上述两个数值进行交叉线性插值，得到的结果为Z。矩阵TAB是第一列与第一行列都包含关键值，而其他的元素包含数据的矩阵。TAB(1,1)的关键值将被忽略。[X0,Y0]中的每点将相应地返回一线性插值。矩阵TAB的第一行与第一列必须是单调的。

例2-39

>>tab = [NaN 1:4; (1:4)' magic(4)]

>>y = table2(tab,[2 3 3.7],[1.3 2.3 4])

查表的结果为：

tab =

NaN 1 2 3 4

1 16 2 3 13

2 5 11 10 8

3 9 7 6 12

4 4 14 15 1

Warning: TABLE2 is obsolete and will be removed in future versions. Use INTERP2 instead.

> In D:\MATLABR12\toolbox\matlab\polyfun\table2.m at line 24

Warning: TABLE1 is obsolete and will be removed in future versions. Use INTERP1 or INTERP1Q instead.

> In D:\MATLABR12\toolbox\matlab\polyfun\table1.m at line 31

In D:\MATLABR12\toolbox\matlab\polyfun\table2.m at line 29

Warning: TABLE1 is obsolete and will be removed in future versions. Use INTERP1 or INTERP1Q instead.

> In D:\MATLABR12\toolbox\matlab\polyfun\table1.m at line 31

In D:\MATLABR12\toolbox\matlab\polyfun\table2.m at line 31

y =

6.8000 10.7000 8.0000

8.4000 6.7000 12.0000

7.4200 12.0200 4.3000

由上面的结果可知，table2是将要废弃的命令。

2.3 数值积分

2.3.1 一元函数的数值积分

函数1 quad、quadl、quad8

功能 数值定积分，自适应Simpleson积分法。

格式 q = quad(fun,a,b) %近似地从a到b计算函数fun的数值积分，误差为10-6。若给fun输入向量x，应返回向量y，即fun是一单值函数。

q = quad(fun,a,b,tol) %用指定的绝对误差tol代替缺省误差。tol越大，函数计算的次数越少，速度越快，但结果精度变小。

q = quad(fun,a,b,tol,trace,p1,p2,…) %将可选参数p1,p2,…等传递给函数fun(x,p1,p2,…)，再作数值积分。若tol=[]或trace=[]，则用缺省值进行计算。

[q,n] = quad(fun,a,b,…) %同时返回函数计算的次数n

… = quadl(fun,a,b,…) %用高精度进行计算，效率可能比quad更好。

… = quad8(fun,a,b,…) %该命令是将废弃的命令，用quadl代替。

例2-40

>>fun = inline(‘’3*x.^2./(x.^3-2*x.^2+3)’);

>>Q1 = quad(fun,0,2)

>>Q2 = quadl(fun,0,2)

计算结果为：

Q1 =

3.7224

Q2 =

3.7224

函数2 trapz

功能 梯形法数值积分

格式 T = trapz(Y) %用等距梯形法近似计算Y的积分。若Y是一向量，则trapz(Y)为Y的积分；若Y是一矩阵，则trapz(Y)为Y的每一列的积分；若Y是一多维阵列，则trapz(Y)沿着Y的第一个非单元集的方向进行计算。

T = trapz(X,Y) %用梯形法计算Y在X点上的积分。若X为一列向量，Y为矩阵，且size(Y,1) = length(X)，则trapz(X,Y)通过Y的第一个非单元集方向进行计算。

T = trapz(…,dim) %沿着dim指定的方向对Y进行积分。若参量中包含X，则应有length(X)=size(Y,dim)。

例2-41

>>X = -1:.1:1;

>>Y = 1./(1+25*X.^2);

>>T = trapz(X,Y)

计算结果为：

T =

0.5492

函数3 rat,rats

功能 有理分式近似。虽然所有的浮点数值都是有理数，有时用简单的有理数字（分子与分母都是较小的整数）近似地表示它们是有必要的。函数rat将试图做到这一点。对于有连续出现的小数的数值，将会用有理式近似表示它们。函数rats调用函数rat，且返回字符串。

格式 [N,D] = rat(X) %对于缺省的误差1.e-6*norm(X(:),1)，返回阵列N与D，使N./D近似为X。

[N,D] = rat(X,tol) %在指定的误差tol范围内，返回阵列N与D，使N./D近似为X。

rat(X)、rat(X…) %在没有输出参量时，简单地显示x的连续分数。

S = rats(X,strlen) %返回一包含简单形式的、X中每一元素的有理近似字符串S，若对于分配的空间中不能显示某一元素，则用星号表示。该元素与X中其他元素进行比较而言较小，但并非是可以忽略。参量strlen为函数rats中返回的字符串元素的长度。缺省值为strlen=13，这允许在78个空格中有6个元素。

S = rats(X) %返回与用MATLAB命令format rat显示 X相同的结果给S。

例2-42

>>s = 1-1/2+1/3-1/4+1/5-1/6+1/7

>>format rat

>>S1 = rats(s)

>>S2 = rat(s)

>>[n,d] = rat(s)

>>PI1 = rats(pi)

>>PI2 = rat(pi)

计算结果为：

s =

0.7595

S1 =

319/420

S2 =

1 + 1/(-4 + 1/(-6 + 1/(-3 + 1/(-5))))

n =

319

d =

420

PI1 =

355/113

PI2 =

3 + 1/(7 + 1/(16))

2.3.2 二元函数重积分的数值计算

函数1 dblquad

功能 矩形区域上的二重积分的数值计算

格式 q = dblquad(fun,xmin,xmax,ymin,ymax) %调用函数quad在区域[xmin,xmax, ymin,ymax]上计算二元函数z=f(x,y)的二重积分。输入向量x，标量y，则f(x,y)必须返回一用于积分的向量。

q = dblquad(fun,xmin,xmax,ymin,ymax,tol) %用指定的精度tol代替缺省精度10-6，再进行计算。

q = dblquad(fun,xmin,xmax,ymin,ymax,tol,method) %用指定的算法method代替缺省算法quad。method的取值有@quadl或用户指定的、与命令quad与quadl有相同调用次序的函数句柄。

q = dblquad(fun,xmin,xmax,ymin,ymax,tol,method,p1,p2,…) %将可选参数p1,p2,..等传递给函数fun(x,y,p1,p2,…)。若tol=[]，method=[]，则使用缺省精度和算法quad。

例2-43

>>fun = inline(’y./sin(x)+x.*exp(y)’);

>>Q = dblquad(fun,1,3,5,7)

计算结果为：

Q =

3.8319e+003

函数2 quad2dggen

功能 任意区域上二元函数的数值积分

格式 q = quad2dggen(fun,xlower,xupper,ymin,ymax) %在由[xlower,xupper, ymin,ymax]指定的区域上计算二元函数z=f(x,y)的二重积分。

q = dblquad(fun,xlower,xupper,ymin,ymax,tol) %用指定的精度tol代替缺省精度10-6，再进行计算。

q = dblquad(fun,xmin,xmax,ymin,ymax,tol,method) %用指定的算法method代替缺省算法。method的取值有缺省算法或用户指定的、与缺省命令有相同调用次序的函数句柄。

q=dblquad(fun,xlower,xupper,ymin,ymax,tol,method,p1,p2,…) %将可选参数p1,p2,..等传递给函数fun(x,y,p1,p2,…)。若tol=[]，method=[]，则使用缺省精度和算法。

例2-44

计算单位圆域上的积分：

先把二重积分转化为二次积分的形式：

f = inline(’exp(-x.^2/2).*sin(x.^2+y)’,’x’,’y’);

xlower = inline(’-sqrt(1-y.^2)’,’y’); xupper = inline(’sqrt(1-y.^2)’,’y’);

Q = quad2dggen(fun,xlower,xupper,-1,1,1e-4)

计算结果为：

Q =

0.5368603818

2.4 常微分方程数值解

函数 ode45、ode23、ode113、ode15s、ode23s、ode23t、ode23tb

功能 常微分方程（ODE）组初值问题的数值解

参数说明：

solver为命令ode45、ode23,ode113,ode15s,ode23s,ode23t,ode23tb之一。

Odefun 为显式常微分方程y’=f(t,y)，或为包含一混合矩阵的方程M(t,y)*y’=f(t,y)。命令ode23只能求解常数混合矩阵的问题；命令ode23t与ode15s可以求解奇异矩阵的问题。

Tspan 积分区间（即求解区间）的向量tspan=[t0,tf]。要获得问题在其他指定时间点t0,t1,t2,…上的解，则令tspan=[t0,t1,t2,…,tf]（要求是单调的）。

Y0 包含初始条件的向量。

Options 用命令odeset设置的可选积分参数。

P1,p2,… 传递给函数odefun的可选参数。

格式 [T,Y] = solver(odefun,tspan,y0) %在区间tspan=[t0,tf]上，从t0到tf，用初始条件y0求解显式微分方程y’=f(t,y)。对于标量t与列向量y，函数f=odefun(t,y)必须返回一f(t,y)的列向量f。解矩阵Y中的每一行对应于返回的时间列向量T中的一个时间点。要获得问题在其他指定时间点t0,t1,t2,…上的解，则令tspan=[t0,t1,t2,…,tf]（要求是单调的）。

[T,Y] = solver(odefun,tspan,y0,options) %用参数options（用命令odeset生成）设置的属性（代替了缺省的积分参数），再进行操作。常用的属性包括相对误差值RelTol（缺省值为1e-3）与绝对误差向量AbsTol（缺省值为每一元素为1e-6）。

[T,Y] =solver(odefun,tspan,y0,options,p1,p2…) 将参数p1,p2,p3,..等传递给函数odefun，再进行计算。若没有参数设置，则令options=[]。

1．求解具体ODE的基本过程：

（1）根据问题所属学科中的规律、定律、公式，用微分方程与初始条件进行描述。

F(y,y’,y’’,…,y(n),t) = 0

y(0)=y0,y’(0)=y1,…,y(n-1)(0)=yn-1

而y=[y;y(1);y(2);…,y(m-1)]，n与m可以不等

（2）运用数学中的变量替换：yn=y(n-1),yn-1=y(n-2),…,y2=y1=y，把高阶（大于2阶）的方程（组）写成一阶微分方程组：，

（3）根据（1）与（2）的结果，编写能计算导数的M-函数文件odefile。

（4）将文件odefile与初始条件传递给求解器Solver中的一个，运行后就可得到ODE的、在指定时间区间上的解列向量y（其中包含y及不同阶的导数）。

2．求解器Solver与方程组的关系表见表2-3。

表2-3

函数指令	含 义	函 数	含 义
求解器 Solver	ode23	普通2-3阶法解ODE	odefile	包含ODE的文件
ode23s	低阶法解刚性ODE	选项	odeset	创建、更改Solver选项
ode23t	解适度刚性ODE	odeget	读取Solver的设置值
ode23tb	低阶法解刚性ODE	输出	odeplot	ODE的时间序列图
ode45	普通4-5阶法解ODE	odephas2	ODE的二维相平面图
ode15s	变阶法解刚性ODE	odephas3	ODE的三维相平面图
ode113	普通变阶法解ODE	odeprint	在命令窗口输出结果

3．因为没有一种算法可以有效地解决所有的ODE问题，为此，MATLAB提供了多种求解器Solver，对于不同的ODE问题，采用不同的Solver。

表2-4 不同求解器Solver的特点

求解器Solver	ODE类型	特点	说明
ode45	非刚性	一步算法；4，5阶Runge-Kutta方程；累计截断误差达(△x)3	大部分场合的首选算法
ode23	非刚性	一步算法；2，3阶Runge-Kutta方程；累计截断误差达(△x)3	使用于精度较低的情形
ode113	非刚性	多步法；Adams算法；高低精度均可到10-3～10-6	计算时间比ode45短
ode23t	适度刚性	采用梯形算法	适度刚性情形
ode15s	刚性	多步法；Gear’s反向数值微分；精度中等	若ode45失效时，可尝试使用
ode23s	刚性	一步法；2阶Rosebrock算法；低精度	当精度较低时，计算时间比ode15s短
ode23tb	刚性	梯形算法；低精度	当精度较低时，计算时间比ode15s短

4．在计算过程中，用户可以对求解指令solver中的具体执行参数进行设置（如绝对误差、相对误差、步长等）。

表2-5 Solver中options的属性

属性名	取值	含义
AbsTol	有效值：正实数或向量 缺省值：1e-6	绝对误差对应于解向量中的所有元素；向量则分别对应于解向量中的每一分量
RelTol	有效值：正实数 缺省值：1e-3	相对误差对应于解向量中的所有元素。在每步(第k步)计算过程中，误差估计为： e(k)<=max(RelTol*abs(y(k)),AbsTol(k))
NormControl	有效值：on、off 缺省值：off	为‘on’时，控制解向量范数的相对误差，使每步计算中，满足：norm(e)<=max(RelTol*norm(y),AbsTol)
Events	有效值：on、off	为‘on’时，返回相应的事件记录
OutputFcn	有效值：odeplot、odephas2、odephas3、odeprint 缺省值：odeplot	若无输出参量，则solver将执行下面操作之一： 画出解向量中各元素随时间的变化； 画出解向量中前两个分量构成的相平面图； 画出解向量中前三个分量构成的三维相空间图； 随计算过程，显示解向量
OutputSel	有效值：正整数向量 缺省值：[]	若不使用缺省设置，则OutputFcn所表现的是那些正整数指定的解向量中的分量的曲线或数据。若为缺省值时，则缺省地按上面情形进行操作
Refine	有效值：正整数k>1 缺省值：k = 1	若k>1，则增加每个积分步中的数据点记录，使解曲线更加的光滑
Jacobian	有效值：on、off 缺省值：off	若为‘on’时，返回相应的ode函数的Jacobi矩阵
Jpattern	有效值：on、off 缺省值：off	为‘on’时，返回相应的ode函数的稀疏Jacobi矩阵
Mass	有效值：none、M、 M(t)、M(t,y) 缺省值：none	M：不随时间变化的常数矩阵 M(t)：随时间变化的矩阵 M(t,y)：随时间、地点变化的矩阵
MaxStep	有效值：正实数 缺省值：tspans/10	最大积分步长

例2-45 求解描述振荡器的经典的Ver der Pol微分方程

y(0)=1，y’(0)=0

令x1=y，x2=dy/dx，则

dx1/dt = x2

dx2/dt = μ(1-x2)-x1

编写函数文件verderpol.m：

function xprime = verderpol(t,x)

global MU

xprime = [x(2);MU*(1-x(1)^2)*x(2)-x(1)];

再在命令窗口中执行：

>>global MU

>>MU = 7;

>>Y0=[1;0]

>>[t,x] = ode45(‘verderpol’,0,40,Y0);

>>x1=x(:,1);x2=x(:,2);

>>plot(t,x1,t,x2)

图形结果为图2-20。



图2-20 Ver der Pol微分方程图

2.5 偏微分方程的数值解

MATLAB提供了一个专门用于求解偏微分方程的工具箱—PDE Toolbox (Paticial Difference Equation )。在本章，我们仅提供一些最简单、经典的偏微分方程，如：椭圆型、双曲型、抛物型等少数的偏微分方程，并给出求解方法。用户可以从中了解其解题基本方法，从而解决相类似的问题。

Matlab能解决的偏微分类型

其中u = u(x,y)，

，

c = c(x,y)，，

2.5.1 单的Poission方程

Poission方程是特殊的椭圆型方程：，

即 c = 1，a = 0，f = -1

Poission的解析解为：。在下面计算中，用求得的数值解与精确解进行比较，看误差如何。

方程求解

（1）问题输入：

c = 1；a = 0；f = 1；%方程的输入。给c，a，f赋值即可

g = 'circleg' % 区域G，内部已经定义为 circleg

b = 'circleb1' % u在区域G的边界上的条件，内部已经定义好

（2）对单位圆进行网格化，对求解区域G作剖分，且作的是三角分划：

[p,e,t] = initmesh（g,'hmax'，1）

（3）迭代求解：

error = []；err = 1；

while err > 0.001，

[p,e,t]=refinemesh（'circleg',p,e,t）；

u=assempde（'circleb1',p,e,t,1,0,1）；

exact=-（p（1，：）^2+p（2，：）^2-1）/4；

err=norm（u-exact',inf）；

error=[error,err]；

end

（4）误差显示与区域G内的剖分点数：

Error: 1.292265e-002. Number of nodes: 25

Error: 4.079923e-003. Number of nodes: 81

Error: 1.221020e-003. Number of nodes: 289

Error: 3.547924e-004. Number of nodes: 1089

（5）结果显示：

subplot（2,2,1），pdemesh（p,e,t）%结果显示

title（'数值解'）

subplot（2,2,2），pdesurf（p,t,u）%精确解显示

title（'精确解'）

subplot（2,2,3），pdesurf（p,t,u-exact'）%与精确解的误差

title（'计算误差'）



图2-21 Poission方程图

2.5.2 双曲型偏微分方程

1．Matlab能求解的类型：

其中，，，，，

2．形传递问题：，

即：c =1; a = 0; f = 0; d = 1

3．方程求解

（1）问题的输入：

c = 1; a = 0; f = 0; d = 1; % 输入方程的系数

g = 'squareg' % 输入方形区域G,内部已经定义好

b = 'sqareb3' % 输入边界条件，即初始条件

（2）对单位矩形G进行网格化：

[p,e,t] = initmesh('squareg')

（3）定解条件和求解时间点：

x = p(1,:)'; y = p(2,:)';

u0 = atan(cos(pi/2*x));

ut0 = 3*sin(pi*x).*exp(sin(pi/2*y));

n = 31;

tlist = linspace(0,5,n);

（4）求解：

uu = hyperbolic(uo, ut0,tlist,b,p,e,t,c,a,f,d);

结果显示：计算过程中的时间点和信息

Time：0.166667

Time：0.333333

Time：4.33333

……

Time：4.66667

Time：4.83333

Time：5

428 successful steps

62 failed attempts

982 function evaluations

1 partial derivatives

142 LU decompositions

981 solutions of linear systems

（5）动画显示：

delta=-1:0.1:1;

[uxy,tn,a2,a3]=tri2grid(p,t,uu(:,1),delta,delta);

gp=[tn;a2;a3];

umax=max(max(uu));

umin=min(min(uu));

newplot;M=moviein(n);

for i=1:n,

pdeplot(p,e,t,'xydata',uu(:,i),'zdata',uu(:,i),…

'mesh','off','xygrid','on','gridparam',gp,…

'colorbar','off','zstyle','continuous');

axis([-1 1 -1 1 umin umax]);

caxis([umin umax]);

M(:,i)=getframe;

end

movie(M,5)

图2-22是动画过程中的一个状态。



图2-22 波动方程动画中的一个状态

2.5.3 抛物型偏微分方程

1．Matlab能求解的类型：

其中，，，，，

2．热传导方程：，

即：c = 1; a = 0; f = 1; d = 1;

3．问题计算

（1）问题的输入：

c = 1; a = 0; f = 1; d = 1; % 输入方程的系数

g = 'squareg'; % 输入方形区域G

b = 'squareb1'; % 输入边界条件

（2）对单位矩形的网格化：

[p,e,t] = initmesh(g);

（3）定解条件和求解的时间点：

u0 = zeros(size(p, 2), 1);

ix = find(sqrt(p(1, :).^2+p(2, :).^2) < 0.4);

u0(ix) = ones(size(ix));

nframes = 20;

tlist=linspace(0,0.1,nframes) % 在时间[0, 0.1]内20个点上计算，生成20帧

（4）求解方程：

u1 = parabolic(u0, tlist, b, p, e, t, c, a, f, d)

计算结果如下：

Time: 0.00526316

Time: 0.0105263

……

Time: 0.0947368

Time: 0.1

75 successful steps

1 failed attempts

154 function evaluations

1 partial derivatives

17 LU decompositions

153 solutions of linear systems

（5）动画显示：

x = linspace(-1,1,31); y = x;

newplot;

Mv = moviein(nframes);

umax=max(max(u1));

umin=min(min(u1));

for j=1:nframes

u=tri2grid(p,t,u1(:,j),tn,a2,a3);i=find(isnan(u));u(i)=zeros(size(i));…

surf(x,y,u);caxis([umin umax]);colormap(cool),…

axis([-1 1 -1 1 0 1]);…

Mv(:,j) = getframe;…

end

movie(Mv,10)

图2-23是动画过程中的瞬间状态。



图2-23 热传导方程动画瞬间状态图