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 热传导方程动画瞬间状态图
